|
|
|
#include <HitEff.h>
Public Member Functions | |
| bool | check2DPartner (unsigned int iidd, std::vector< TrajectoryMeasurement > traj) |
| double | checkConsistency (StripClusterParameterEstimator::LocalValues parameters, double xx, double xerr) |
| unsigned int | checkLayer (unsigned int iidd) |
| HitEff (const edm::ParameterSet &conf) | |
| bool | isDoubleSided (unsigned int iidd) const |
| virtual | ~HitEff () |
Private Member Functions | |
| virtual void | analyze (const edm::Event &e, const edm::EventSetup &c) |
| virtual void | beginJob () |
| virtual void | endJob () |
Private Attributes | |
| unsigned int | bunchx |
| float | chi2 |
| float | ClusterLocErrX |
| float | ClusterLocErrY |
| float | ClusterLocX |
| float | ClusterLocY |
| float | ClusterStoN |
| edm::ParameterSet | conf_ |
| bool | DEBUG |
| float | dedx |
| int | dedxNOM |
| unsigned int | event |
| int | events |
| int | EventTrackCKF |
| unsigned int | Id |
| int | istep |
| unsigned int | layers |
| unsigned int | ModIsBad |
| int | nHits |
| int | nLostHits |
| float | p |
| float | pT |
| float | ResX |
| float | ResXSig |
| unsigned int | run |
| unsigned int | SiStripQualBad |
| float | timeDT |
| int | timeDTDOF |
| float | timeDTErr |
| float | timeECAL |
| int | tquality |
| TTree * | traj |
| float | TrajGlbX |
| float | TrajGlbY |
| float | TrajGlbZ |
| unsigned int | trajHitValid |
| float | TrajLocAngleX |
| float | TrajLocAngleY |
| float | TrajLocErrX |
| float | TrajLocErrY |
| float | TrajLocX |
| float | TrajLocY |
| unsigned int | whatlayer |
| bool | withinAcceptance |
| HitEff::HitEff | ( | const edm::ParameterSet & | conf | ) | [explicit] |
Definition at line 74 of file HitEff.cc.
References conf_, DEBUG, edm::ParameterSet::getParameter(), and layers.
: conf_(conf) { layers =conf_.getParameter<int>("Layer"); DEBUG = conf_.getParameter<bool>("Debug"); }
| void HitEff::analyze | ( | const edm::Event & | e, |
| const edm::EventSetup & | c | ||
| ) | [private, virtual] |
Implements edm::EDAnalyzer.
Definition at line 134 of file HitEff.cc.
References abs, anyDirection, edm::EventBase::bunchCrossing(), bunchx, check2DPartner(), checkConsistency(), checkLayer(), chi2, ClusterLocErrX, ClusterLocErrY, ClusterLocX, ClusterLocY, ClusterStoN, conf_, gather_cfg::cout, DEBUG, edm::EventID::event(), events, EventTrackCKF, edm::EventSetup::get(), edm::Event::getByLabel(), edm::ParameterSet::getParameter(), i, edm::EventBase::id(), Id, collect_tpl::input, ires, isDoubleSided(), layers, ClusterParameterEstimator< T >::localParameters(), LayerMeasurements::measurements(), ModIsBad, reco::MuonTime::nDof, nHits, nLostHits, p, Parameters::parameters, edm::Handle< T >::product(), LargeD0_PixelPairStep_cff::propagator, pT, dt_dqm_sourceclient_common_cff::reco, ResX, ResXSig, run, edm::EventID::run(), TECDetId::side(), SiStripClusterInfo::signalOverNoise(), SiStripQualBad, mathSSE::sqrt(), StripSubdetector::TEC, reco::MuonTime::timeAtIpInOut, reco::MuonTime::timeAtIpInOutErr, timeDT, timeDTDOF, timeDTErr, timeECAL, tmp, patCandidatesForDimuonsSequences_cff::tracker, traj, TrajGlbX, TrajGlbY, TrajGlbZ, trajHitValid, TrajLocAngleX, TrajLocAngleY, TrajLocErrX, TrajLocErrY, TrajLocX, TrajLocY, whatlayer, and withinAcceptance.
{
// bool DEBUG = false;
if (DEBUG) cout << "beginning analyze from HitEff" << endl;
using namespace edm;
using namespace reco;
// Step A: Get Inputs
int run_nr = e.id().run();
int ev_nr = e.id().event();
int bunch_nr = e.bunchCrossing();
//CombinatoriaTrack
edm::Handle<reco::TrackCollection> trackCollectionCKF;
edm::InputTag TkTagCKF = conf_.getParameter<edm::InputTag>("combinatorialTracks");
e.getByLabel(TkTagCKF,trackCollectionCKF);
edm::Handle<std::vector<Trajectory> > TrajectoryCollectionCKF;
edm::InputTag TkTrajCKF = conf_.getParameter<edm::InputTag>("trajectories");
e.getByLabel(TkTrajCKF,TrajectoryCollectionCKF);
// Clusters
// get the SiStripClusters from the event
edm::Handle< edmNew::DetSetVector<SiStripCluster> > theClusters;
e.getByLabel("siStripClusters", theClusters);
//get tracker geometry
edm::ESHandle<TrackerGeometry> tracker;
es.get<TrackerDigiGeometryRecord>().get(tracker);
const TrackerGeometry * tkgeom=&(* tracker);
//get Cluster Parameter Estimator
//std::string cpe = conf_.getParameter<std::string>("StripCPE");
edm::ESHandle<StripClusterParameterEstimator> parameterestimator;
es.get<TkStripCPERecord>().get("StripCPEfromTrackAngle", parameterestimator);
const StripClusterParameterEstimator &stripcpe(*parameterestimator);
// get the SiStripQuality records
edm::ESHandle<SiStripQuality> SiStripQuality_;
try {
es.get<SiStripQualityRcd>().get("forCluster",SiStripQuality_);
}
catch (...) {
es.get<SiStripQualityRcd>().get(SiStripQuality_);
}
edm::ESHandle<MagneticField> magFieldHandle;
es.get<IdealMagneticFieldRecord>().get(magFieldHandle);
const MagneticField* magField_ = magFieldHandle.product();
// get the list of module IDs with FED-detected errors
edm::Handle< DetIdCollection > fedErrorIds;
e.getByLabel("siStripDigis", fedErrorIds );
ESHandle<MeasurementTracker> measurementTrackerHandle;
es.get<CkfComponentsRecord>().get(measurementTrackerHandle);
edm::ESHandle<Chi2MeasurementEstimatorBase> est;
es.get<TrackingComponentsRecord>().get("Chi2",est);
edm::ESHandle<Propagator> prop;
es.get<TrackingComponentsRecord>().get("PropagatorWithMaterial",prop);
const Propagator* thePropagator = prop.product();
events++;
// *************** SiStripCluster Collection
const edmNew::DetSetVector<SiStripCluster>& input = *theClusters;
//go through clusters to write out global position of good clusters for the layer understudy for comparison
// Loop through clusters just to print out locations
// Commented out to avoid discussion, should really be deleted.
/*
for (edmNew::DetSetVector<SiStripCluster>::const_iterator DSViter = input.begin(); DSViter != input.end(); DSViter++) {
// DSViter is a vector of SiStripClusters located on a single module
unsigned int ClusterId = DSViter->id();
DetId ClusterDetId(ClusterId);
const StripGeomDetUnit * stripdet=(const StripGeomDetUnit*)tkgeom->idToDetUnit(ClusterDetId);
edmNew::DetSet<SiStripCluster>::const_iterator begin=DSViter->begin();
edmNew::DetSet<SiStripCluster>::const_iterator end =DSViter->end();
for(edmNew::DetSet<SiStripCluster>::const_iterator iter=begin;iter!=end;++iter) {
//iter is a single SiStripCluster
StripClusterParameterEstimator::LocalValues parameters=stripcpe.localParameters(*iter,*stripdet);
const Surface* surface;
surface = &(tracker->idToDet(ClusterDetId)->surface());
LocalPoint lp = parameters.first;
GlobalPoint gp = surface->toGlobal(lp);
unsigned int layer = checkLayer(ClusterId);
if(DEBUG) cout << "Found hit in cluster collection layer = " << layer << " with id = " << ClusterId << " local X position = " << lp.x() << " +- " << sqrt(parameters.second.xx()) << " matched/stereo/rphi = " << ((ClusterId & 0x3)==0) << "/" << ((ClusterId & 0x3)==1) << "/" << ((ClusterId & 0x3)==2) << endl;
}
}
*/
// Tracking
const reco::TrackCollection *tracksCKF=trackCollectionCKF.product();
if (DEBUG) cout << "number ckf tracks found = " << tracksCKF->size() << endl;
//if (tracksCKF->size() == 1 ){
if (tracksCKF->size() > 0 && tracksCKF->size()<100) {
if (DEBUG) cout << "starting checking good event with < 100 tracks" << endl;
EventTrackCKF++;
/*
//get dEdx info if available
Handle<ValueMap<DeDxData> > dEdxUncalibHandle;
if (e.getByLabel("dedxMedianCTF", dEdxUncalibHandle)) {
const ValueMap<DeDxData> dEdxTrackUncalib = *dEdxUncalibHandle.product();
reco::TrackRef itTrack = reco::TrackRef( trackCollectionCKF, 0 );
dedx = dEdxTrackUncalib[itTrack].dEdx();
dedxNOM = dEdxTrackUncalib[itTrack].numberOfMeasurements();
} else {
dedx = -999.0; dedxNOM = -999;
}
*/
//get muon and ecal timing info if available
Handle<MuonCollection> muH;
if(e.getByLabel("muonsWitht0Correction",muH)){
const MuonCollection & muonsT0 = *muH.product();
if(muonsT0.size()!=0) {
MuonTime mt0 = muonsT0[0].time();
timeDT = mt0.timeAtIpInOut;
timeDTErr = mt0.timeAtIpInOutErr;
timeDTDOF = mt0.nDof;
bool hasCaloEnergyInfo = muonsT0[0].isEnergyValid();
if (hasCaloEnergyInfo) timeECAL = muonsT0[0].calEnergy().ecal_time;
}
} else {
timeDT = -999.0; timeDTErr = -999.0; timeDTDOF = -999; timeECAL = -999.0;
}
// actually should do a loop over all the tracks in the event here
for (vector<Trajectory>::const_iterator itraj = TrajectoryCollectionCKF.product()->begin();
itraj != TrajectoryCollectionCKF.product()->end();
itraj++) {
// for each track, fill some variables such as number of hits and momentum
nHits = itraj->foundHits();
nLostHits = itraj->lostHits();
chi2 = (itraj->chiSquared()/itraj->ndof());
pT = sqrt( ( itraj->lastMeasurement().updatedState().globalMomentum().x() *
itraj->lastMeasurement().updatedState().globalMomentum().x()) +
( itraj->lastMeasurement().updatedState().globalMomentum().y() *
itraj->lastMeasurement().updatedState().globalMomentum().y()) );
p = itraj->lastMeasurement().updatedState().globalMomentum().mag();
//Put in code to check track quality
std::vector<TrajectoryMeasurement> TMeas=itraj->measurements();
vector<TrajectoryMeasurement>::iterator itm;
double xloc = 0.;
double yloc = 0.;
double xErr = 0.;
double yErr = 0.;
double angleX = -999.;
double angleY = -999.;
double xglob,yglob,zglob;
for (itm=TMeas.begin();itm!=TMeas.end();itm++){
ConstReferenceCountingPointer<TransientTrackingRecHit> theInHit;
theInHit = (*itm).recHit();
if(DEBUG) cout << "theInHit is valid = " << theInHit->isValid() << endl;
unsigned int iidd = theInHit->geographicalId().rawId();
unsigned int TKlayers = checkLayer(iidd);
if (DEBUG) cout << "TKlayer from trajectory: " << TKlayers << " from module = " << iidd << " matched/stereo/rphi = " << ((iidd & 0x3)==0) << "/" << ((iidd & 0x3)==1) << "/" << ((iidd & 0x3)==2) << endl;
// If Trajectory measurement from TOB 6 or TEC 9, skip it because it's always valid they are filled later
if ( TKlayers == 10 || TKlayers == 22 ) {
if (DEBUG) cout << "skipping original TM for TOB 6 or TEC 9" << endl;
continue;
}
// Make vector of TrajectoryAtInvalidHits to hold the trajectories
std::vector<TrajectoryAtInvalidHit> TMs;
// Make AnalyticalPropagator to use in TAVH constructor
AnalyticalPropagator propagator(magField_,anyDirection);
// for double sided layers check both sensors--if no hit was found on either sensor surface,
// the trajectory measurements only have one invalid hit entry on the matched surface
// so get the TrajectoryAtInvalidHit for both surfaces and include them in the study
if (isDoubleSided(iidd) && ((iidd & 0x3)==0) ) {
// do hit eff check twice--once for each sensor
//add a TM for each surface
TMs.push_back(TrajectoryAtInvalidHit(*itm,tkgeom, propagator, 1));
TMs.push_back(TrajectoryAtInvalidHit(*itm,tkgeom, propagator, 2));
} else if ( isDoubleSided(iidd) && (!check2DPartner(iidd, TMeas)) ) {
// if only one hit was found the trajectory measurement is on that sensor surface, and the other surface from
// the matched layer should be added to the study as well
TMs.push_back(TrajectoryAtInvalidHit(*itm,tkgeom, propagator, 1));
TMs.push_back(TrajectoryAtInvalidHit(*itm,tkgeom, propagator, 2));
if (DEBUG) cout << " found a hit with a missing partner" << endl;
} else {
//only add one TM for the single surface and the other will be added in the next iteration
TMs.push_back(TrajectoryAtInvalidHit(*itm,tkgeom, propagator));
}
//Now check for tracks at TOB6 and TEC9
// to make sure we only propagate on the last TOB5 hit check the next entry isn't also in TOB5
// to avoid bias, make sure the TOB5 hit is valid (an invalid hit on TOB5 could only exist with a valid hit on TOB6)
bool isValid = theInHit->isValid();
bool isLast = (itm==(TMeas.end()-1));
bool isLastTOB5 = true;
if (!isLast) {
if ( checkLayer((++itm)->recHit()->geographicalId().rawId()) == 9 ) isLastTOB5 = false;
else isLastTOB5 = true;
--itm;
}
if ( TKlayers==9 && isValid && isLastTOB5 ) {
// if ( TKlayers==9 && itm==TMeas.rbegin()) {
// if ( TKlayers==9 && (itm==TMeas.back()) ) { // to check for only the last entry in the trajectory for propagation
std::vector< BarrelDetLayer*> barrelTOBLayers = measurementTrackerHandle->geometricSearchTracker()->tobLayers() ;
const DetLayer* tob6 = barrelTOBLayers[barrelTOBLayers.size()-1];
const MeasurementEstimator* estimator = est.product();
const LayerMeasurements* theLayerMeasurements = new LayerMeasurements(&*measurementTrackerHandle);
const TrajectoryStateOnSurface tsosTOB5 = itm->updatedState();
vector<TrajectoryMeasurement> tmp = theLayerMeasurements->measurements(*tob6, tsosTOB5, *thePropagator, *estimator);
if ( !tmp.empty()) {
if (DEBUG) cout << "size of TM from propagation = " << tmp.size() << endl;
// take the last of the TMs, which is always an invalid hit
// if no detId is available, ie detId==0, then no compatible layer was crossed
// otherwise, use that TM for the efficiency measurement
TrajectoryMeasurement tob6TM(tmp.back());
ConstReferenceCountingPointer<TransientTrackingRecHit> tob6Hit;
tob6Hit = tob6TM.recHit();
if (tob6Hit->geographicalId().rawId()!=0) {
if (DEBUG) cout << "tob6 hit actually being added to TM vector" << endl;
TMs.push_back(TrajectoryAtInvalidHit(tob6TM,tkgeom, propagator));
}
}
}
bool isLastTEC8 = true;
if (!isLast) {
if ( checkLayer((++itm)->recHit()->geographicalId().rawId()) == 21 ) isLastTEC8 = false;
else isLastTEC8 = true;
--itm;
}
if ( TKlayers==21 && isValid && isLastTEC8 ) {
std::vector< ForwardDetLayer*> posTecLayers = measurementTrackerHandle->geometricSearchTracker()->posTecLayers() ;
const DetLayer* tec9pos = posTecLayers[posTecLayers.size()-1];
std::vector< ForwardDetLayer*> negTecLayers = measurementTrackerHandle->geometricSearchTracker()->negTecLayers() ;
const DetLayer* tec9neg = negTecLayers[negTecLayers.size()-1];
const MeasurementEstimator* estimator = est.product();
const LayerMeasurements* theLayerMeasurements = new LayerMeasurements(&*measurementTrackerHandle);
const TrajectoryStateOnSurface tsosTEC9 = itm->updatedState();
// check if track on positive or negative z
if (!iidd == StripSubdetector::TEC) cout << "there is a problem with TEC 9 extrapolation" << endl;
TECDetId tecdetid(iidd);
//cout << " tec9 id = " << iidd << " and side = " << tecdetid.side() << endl;
vector<TrajectoryMeasurement> tmp;
if ( tecdetid.side() == 1 ) {
tmp = theLayerMeasurements->measurements(*tec9neg, tsosTEC9, *thePropagator, *estimator);
//cout << "on negative side" << endl;
}
if ( tecdetid.side() == 2 ) {
tmp = theLayerMeasurements->measurements(*tec9pos, tsosTEC9, *thePropagator, *estimator);
//cout << "on positive side" << endl;
}
if ( !tmp.empty()) {
// take the last of the TMs, which is always an invalid hit
// if no detId is available, ie detId==0, then no compatible layer was crossed
// otherwise, use that TM for the efficiency measurement
TrajectoryMeasurement tec9TM(tmp.back());
ConstReferenceCountingPointer<TransientTrackingRecHit> tec9Hit;
tec9Hit = tec9TM.recHit();
unsigned int tec9id = tec9Hit->geographicalId().rawId();
if (DEBUG) cout << "tec9id = " << tec9id << " is Double sided = " << isDoubleSided(tec9id) << " and 0x3 = " << (tec9id & 0x3) << endl;
if (tec9Hit->geographicalId().rawId()!=0) {
if (DEBUG) cout << "tec9 hit actually being added to TM vector" << endl;
// in tec the hit can be single or doubled sided. whenever the invalid hit at the end of vector of TMs is
// double sided it is always on the matched surface, so we need to split it into the true sensor surfaces
if (isDoubleSided(tec9id)) {
TMs.push_back(TrajectoryAtInvalidHit(tec9TM,tkgeom, propagator, 1));
TMs.push_back(TrajectoryAtInvalidHit(tec9TM,tkgeom, propagator, 2));
} else
TMs.push_back(TrajectoryAtInvalidHit(tec9TM,tkgeom, propagator));
}
} //else cout << "tec9 tmp empty" << endl;
}
// Modules Constraints
for(std::vector<TrajectoryAtInvalidHit>::const_iterator TM=TMs.begin();TM!=TMs.end();++TM) {
// --> Get trajectory from combinatedState
iidd = TM->monodet_id();
if (DEBUG) cout << "setting iidd = " << iidd << " before checking efficiency and ";
xloc = TM->localX();
yloc = TM->localY();
xErr = TM->localErrorX();
yErr = TM->localErrorY();
angleX = atan( TM->localDxDz() );
angleY = atan( TM->localDyDz() );
xglob = TM->globalX();
yglob = TM->globalY();
zglob = TM->globalZ();
withinAcceptance = TM->withinAcceptance();
trajHitValid = TM->validHit();
// reget layer from iidd here, to account for TOB 6 and TEC 9 TKlayers being off
TKlayers = checkLayer(iidd);
if ( (layers == TKlayers) || (layers == 0) ) { // Look at the layer not used to reconstruct the track
whatlayer = TKlayers;
if (DEBUG) cout << "Looking at layer under study" << endl;
TrajGlbX = 0.0; TrajGlbY = 0.0; TrajGlbZ = 0.0; ModIsBad = 2; Id = 0; SiStripQualBad = 0;
run = 0; event = 0; TrajLocX = 0.0; TrajLocY = 0.0; TrajLocErrX = 0.0; TrajLocErrY = 0.0;
TrajLocAngleX = -999.0; TrajLocAngleY = -999.0; ResX = 0.0; ResXSig = 0.0;
ClusterLocX = 0.0; ClusterLocY = 0.0; ClusterLocErrX = 0.0; ClusterLocErrY = 0.0; ClusterStoN = 0.0;
bunchx = 0;
// RPhi RecHit Efficiency
if (input.size() > 0 ) {
if (DEBUG) cout << "Checking clusters with size = " << input.size() << endl;
int nClusters = 0;
std::vector< std::vector<float> > VCluster_info; //fill with X residual, X residual pull, local X, sig(X), local Y, sig(Y), StoN
for (edmNew::DetSetVector<SiStripCluster>::const_iterator DSViter = input.begin(); DSViter != input.end(); DSViter++) {
// DSViter is a vector of SiStripClusters located on a single module
//if (DEBUG) cout << "the ID from the DSViter = " << DSViter->id() << endl;
unsigned int ClusterId = DSViter->id();
if (ClusterId == iidd) {
if (DEBUG) cout << "found (ClusterId == iidd) with ClusterId = " << ClusterId << " and iidd = " << iidd << endl;
DetId ClusterDetId(ClusterId);
const StripGeomDetUnit * stripdet=(const StripGeomDetUnit*)tkgeom->idToDetUnit(ClusterDetId);
for(edmNew::DetSet<SiStripCluster>::const_iterator iter=DSViter->begin();iter!=DSViter->end();++iter) {
//iter is a single SiStripCluster
StripClusterParameterEstimator::LocalValues parameters=stripcpe.localParameters(*iter,*stripdet);
float res = (parameters.first.x() - xloc);
float sigma = checkConsistency(parameters , xloc, xErr);
// The consistency is probably more accurately measured with the Chi2MeasurementEstimator. To use it
// you need a TransientTrackingRecHit instead of the cluster
//theEstimator= new Chi2MeasurementEstimator(30);
//const Chi2MeasurementEstimator *theEstimator(100);
//theEstimator->estimate(TM->tsos(), TransientTrackingRecHit);
SiStripClusterInfo clusterInfo = SiStripClusterInfo(*iter, es, ClusterId);
// signal to noise from SiStripClusterInfo not working in 225. I'll fix this after the interface
// redesign in 300 -ku
//float cluster_info[7] = {res, sigma, parameters.first.x(), sqrt(parameters.second.xx()), parameters.first.y(), sqrt(parameters.second.yy()), signal_to_noise};
std::vector< float > cluster_info;
cluster_info.push_back(res);
cluster_info.push_back(sigma);
cluster_info.push_back(parameters.first.x());
cluster_info.push_back(sqrt(parameters.second.xx()));
cluster_info.push_back(parameters.first.y());
cluster_info.push_back(sqrt(parameters.second.yy()));
cout << "before getting signal over noise" << endl;
cluster_info.push_back( clusterInfo.signalOverNoise() );
//cluster_info.push_back( clusterInfo.getSignalOverNoise() );
cout << "after getting signal over noise" << endl;
VCluster_info.push_back(cluster_info);
nClusters++;
if (DEBUG) cout << "Have ID match. residual = " << VCluster_info.back()[0] << " res sigma = " << VCluster_info.back()[1] << endl;
if (DEBUG) cout << "trajectory measurement compatability estimate = " << (*itm).estimate() << endl;
if (DEBUG) cout << "hit position = " << parameters.first.x() << " hit error = " << sqrt(parameters.second.xx()) << " trajectory position = " << xloc << " traj error = " << xErr << endl;
}
}
}
float FinalResSig = 1000.0;
float FinalCluster[7]= {1000.0, 1000.0, 0.0, 0.0, 0.0, 0.0, 0.0};
if (nClusters > 0) {
if (DEBUG) cout << "found clusters > 0" << endl;
if (nClusters > 1) {
//get the smallest one
vector< vector<float> >::iterator ires;
for (ires=VCluster_info.begin(); ires!=VCluster_info.end(); ires++){
if ( abs((*ires)[1]) < abs(FinalResSig)) {
FinalResSig = (*ires)[1];
for (unsigned int i = 0; i<ires->size(); i++) {
if (DEBUG) cout << "filling final cluster. i = " << i << " before fill FinalCluster[i]=" << FinalCluster[i] << " and (*ires)[i] =" << (*ires)[i] << endl;
FinalCluster[i] = (*ires)[i];
if (DEBUG) cout << "filling final cluster. i = " << i << " after fill FinalCluster[i]=" << FinalCluster[i] << " and (*ires)[i] =" << (*ires)[i] << endl;
}
}
if (DEBUG) cout << "iresidual = " << (*ires)[0] << " isigma = " << (*ires)[1] << " and FinalRes = " << FinalCluster[0] << endl;
}
}
else {
FinalResSig = VCluster_info.at(0)[1];
for (unsigned int i = 0; i<VCluster_info.at(0).size(); i++) {
FinalCluster[i] = VCluster_info.at(0)[i];
}
}
nClusters=0;
VCluster_info.clear();
}
if (DEBUG) cout << "Final residual in X = " << FinalCluster[0] << "+-" << (FinalCluster[0]/FinalResSig) << endl;
if (DEBUG) cout << "Checking location of trajectory: abs(yloc) = " << abs(yloc) << " abs(xloc) = " << abs(xloc) << endl;
if (DEBUG) cout << "Checking location of cluster hit: yloc = " << FinalCluster[4] << "+-" << FinalCluster[5] << " xloc = " << FinalCluster[2] << "+-" << FinalCluster[3] << endl;
if (DEBUG) cout << "Final cluster signal to noise = " << FinalCluster[6] << endl;
float exclusionWidth = 0.4;
float TOBexclusion = 0.0;
float TECexRing5 = -0.89;
float TECexRing6 = -0.56;
float TECexRing7 = 0.60;
//Added by Chris Edelmaier to do TEC bonding exclusion
int subdetector = ((iidd>>25) & 0x7);
int ringnumber = ((iidd>>5) & 0x7);
//New TOB and TEC bonding region exclusion zone
if((TKlayers >= 5 && TKlayers < 11)||((subdetector == 6)&&( (ringnumber >= 5)&&(ringnumber <=7) ))) {
//There are only 2 cases that we need to exclude for
float highzone = 0.0;
float lowzone = 0.0;
float higherr = yloc + 5.0*yErr;
float lowerr = yloc - 5.0*yErr;
if(TKlayers >= 5 && TKlayers < 11) {
//TOB zone
highzone = TOBexclusion + exclusionWidth;
lowzone = TOBexclusion - exclusionWidth;
} else if (ringnumber == 5) {
//TEC ring 5
highzone = TECexRing5 + exclusionWidth;
lowzone = TECexRing5 - exclusionWidth;
} else if (ringnumber == 6) {
//TEC ring 6
highzone = TECexRing6 + exclusionWidth;
lowzone = TECexRing6 - exclusionWidth;
} else if (ringnumber == 7) {
//TEC ring 7
highzone = TECexRing7 + exclusionWidth;
lowzone = TECexRing7 - exclusionWidth;
}
//Now that we have our exclusion region, we just have to properly identify it
if((highzone <= higherr)&&(highzone >= lowerr)) withinAcceptance = false;
if((lowzone >= lowerr)&&(lowzone <= higherr)) withinAcceptance = false;
if((higherr <= highzone)&&(higherr >= lowzone)) withinAcceptance = false;
if((lowerr >= lowzone)&&(lowerr <= highzone)) withinAcceptance = false;
}
// fill ntuple varibles
//get global position from module id number iidd
TrajGlbX = xglob;
TrajGlbY = yglob;
TrajGlbZ = zglob;
Id = iidd;
run = run_nr;
event = ev_nr;
bunchx = bunch_nr;
//if ( SiStripQuality_->IsModuleBad(iidd) ) {
if ( SiStripQuality_->getBadApvs(iidd)!=0 ) {
SiStripQualBad = 1;
if(DEBUG) cout << "strip is bad from SiStripQuality" << endl;
} else {
SiStripQualBad = 0;
if(DEBUG) cout << "strip is good from SiStripQuality" << endl;
}
//check for FED-detected errors and include those in SiStripQualBad
for (unsigned int ii=0;ii< fedErrorIds->size();ii++) {
if (iidd == (*fedErrorIds)[ii].rawId() )
SiStripQualBad = 1;
}
TrajLocX = xloc;
TrajLocY = yloc;
TrajLocErrX = xErr;
TrajLocErrY = yErr;
TrajLocAngleX = angleX;
TrajLocAngleY = angleY;
ResX = FinalCluster[0];
ResXSig = FinalResSig;
if (FinalResSig != FinalCluster[1]) if (DEBUG) cout << "Problem with best cluster selection because FinalResSig = " << FinalResSig << " and FinalCluster[1] = " << FinalCluster[1] << endl;
ClusterLocX = FinalCluster[2];
ClusterLocY = FinalCluster[4];
ClusterLocErrX = FinalCluster[3];
ClusterLocErrY = FinalCluster[5];
ClusterStoN = FinalCluster[6];
if (DEBUG) cout << "before check good" << endl;
if ( FinalResSig < 999.0) { //could make requirement on track/hit consistency, but for
//now take anything with a hit on the module
if (DEBUG) cout << "hit being counted as good " << FinalCluster[0] << " FinalRecHit " <<
iidd << " TKlayers " << TKlayers << " xloc " << xloc << " yloc " << yloc << " module " << iidd <<
" matched/stereo/rphi = " << ((iidd & 0x3)==0) << "/" << ((iidd & 0x3)==1) << "/" << ((iidd & 0x3)==2) << endl;
ModIsBad = 0;
traj->Fill();
}
else {
if (DEBUG) cout << "hit being counted as bad ######### Invalid RPhi FinalResX " << FinalCluster[0] << " FinalRecHit " <<
iidd << " TKlayers " << TKlayers << " xloc " << xloc << " yloc " << yloc << " module " << iidd <<
" matched/stereo/rphi = " << ((iidd & 0x3)==0) << "/" << ((iidd & 0x3)==1) << "/" << ((iidd & 0x3)==2) << endl;
ModIsBad = 1;
traj->Fill();
if (DEBUG) cout << " RPhi Error " << sqrt(xErr*xErr + yErr*yErr) << " ErrorX " << xErr << " yErr " << yErr << endl;
} if (DEBUG) cout << "after good location check" << endl;
} if (DEBUG) cout << "after list of clusters" << endl;
} if (DEBUG) cout << "After layers=TKLayers if" << endl;
} if (DEBUG) cout << "After looping over TrajAtValidHit list" << endl;
} if (DEBUG) cout << "end TMeasurement loop" << endl;
} if (DEBUG) cout << "end of trajectories loop" << endl;
}
}
| void HitEff::beginJob | ( | void | ) | [private, virtual] |
Reimplemented from edm::EDAnalyzer.
Definition at line 84 of file HitEff.cc.
References bunchx, chi2, ClusterLocErrX, ClusterLocErrY, ClusterLocX, ClusterLocY, ClusterStoN, dedx, dedxNOM, event, events, EventTrackCKF, Id, istep, ModIsBad, nHits, nLostHits, p, pT, ResX, ResXSig, run, SiStripQualBad, timeDT, timeDTDOF, timeDTErr, timeECAL, tquality, traj, TrajGlbX, TrajGlbY, TrajGlbZ, trajHitValid, TrajLocAngleX, TrajLocAngleY, TrajLocErrX, TrajLocErrY, TrajLocX, TrajLocY, whatlayer, and withinAcceptance.
{
edm::Service<TFileService> fs;
traj = fs->make<TTree>("traj","tree of trajectory positions");
traj->Branch("TrajGlbX",&TrajGlbX,"TrajGlbX/F");
traj->Branch("TrajGlbY",&TrajGlbY,"TrajGlbY/F");
traj->Branch("TrajGlbZ",&TrajGlbZ,"TrajGlbZ/F");
traj->Branch("TrajLocX",&TrajLocX,"TrajLocX/F");
traj->Branch("TrajLocY",&TrajLocY,"TrajLocY/F");
traj->Branch("TrajLocErrX",&TrajLocErrX,"TrajLocErrX/F");
traj->Branch("TrajLocErrY",&TrajLocErrY,"TrajLocErrY/F");
traj->Branch("TrajLocAngleX",&TrajLocAngleX,"TrajLocAngleX/F");
traj->Branch("TrajLocAngleY",&TrajLocAngleY,"TrajLocAngleY/F");
traj->Branch("ClusterLocX",&ClusterLocX,"ClusterLocX/F");
traj->Branch("ClusterLocY",&ClusterLocY,"ClusterLocY/F");
traj->Branch("ClusterLocErrX",&ClusterLocErrX,"ClusterLocErrX/F");
traj->Branch("ClusterLocErrY",&ClusterLocErrY,"ClusterLocErrY/F");
traj->Branch("ClusterStoN",&ClusterStoN,"ClusterStoN/F");
traj->Branch("ResX",&ResX,"ResX/F");
traj->Branch("ResXSig",&ResXSig,"ResXSig/F");
traj->Branch("ModIsBad",&ModIsBad,"ModIsBad/i");
traj->Branch("SiStripQualBad",&SiStripQualBad,"SiStripQualBad/i");
traj->Branch("withinAcceptance",&withinAcceptance,"withinAcceptance/O");
traj->Branch("nHits",&nHits,"nHits/I");
traj->Branch("nLostHits",&nLostHits,"nLostHits/I");
traj->Branch("chi2",&chi2,"chi2/F");
traj->Branch("p",&p,"p/F");
traj->Branch("pT",&pT,"pT/F");
traj->Branch("trajHitValid", &trajHitValid, "trajHitValid/i");
traj->Branch("Id",&Id,"Id/i");
traj->Branch("run",&run,"run/i");
traj->Branch("event",&event,"event/i");
traj->Branch("layer",&whatlayer,"layer/i");
traj->Branch("timeDT",&timeDT,"timeDT/F");
traj->Branch("timeDTErr",&timeDTErr,"timeDTErr/F");
traj->Branch("timeDTDOF",&timeDTDOF,"timeDTDOF/I");
traj->Branch("timeECAL",&timeECAL,"timeECAL/F");
traj->Branch("dedx",&dedx,"dedx/F");
traj->Branch("dedxNOM",&dedxNOM,"dedxNOM/I");
traj->Branch("tquality",&tquality,"tquality/I");
traj->Branch("istep",&istep,"istep/I");
traj->Branch("bunchx",&bunchx,"bunchx/I");
events = 0;
EventTrackCKF = 0;
}
| bool HitEff::check2DPartner | ( | unsigned int | iidd, |
| std::vector< TrajectoryMeasurement > | traj | ||
| ) |
Definition at line 716 of file HitEff.cc.
Referenced by analyze().
{
unsigned int partner_iidd = 0;
bool found2DPartner = false;
// first get the id of the other detector
if ((iidd & 0x3)==1) partner_iidd = iidd+1;
if ((iidd & 0x3)==2) partner_iidd = iidd-1;
// next look in the trajectory measurements for a measurement from that detector
// loop through trajectory measurements to find the partner_iidd
for (std::vector<TrajectoryMeasurement>::const_iterator iTM=traj.begin(); iTM!=traj.end(); ++iTM) {
if (iTM->recHit()->geographicalId().rawId()==partner_iidd) {
found2DPartner = true;
}
}
return found2DPartner;
}
| double HitEff::checkConsistency | ( | StripClusterParameterEstimator::LocalValues | parameters, |
| double | xx, | ||
| double | xerr | ||
| ) |
Definition at line 677 of file HitEff.cc.
References abs, error, and mathSSE::sqrt().
Referenced by analyze().
{
double error = sqrt(parameters.second.xx() + xerr*xerr);
double separation = abs(parameters.first.x() - xx);
double consistency = separation/error;
return consistency;
}
| unsigned int HitEff::checkLayer | ( | unsigned int | iidd | ) |
Definition at line 732 of file HitEff.cc.
References TIBDetId::layer(), TOBDetId::layer(), strip(), DetId::subdetId(), StripSubdetector::TEC, StripSubdetector::TIB, StripSubdetector::TID, StripSubdetector::TOB, TIDDetId::wheel(), and TECDetId::wheel().
Referenced by analyze().
{
StripSubdetector strip=StripSubdetector(iidd);
unsigned int subid=strip.subdetId();
if (subid == StripSubdetector::TIB) {
TIBDetId tibid(iidd);
return tibid.layer();
}
if (subid == StripSubdetector::TOB) {
TOBDetId tobid(iidd);
return tobid.layer() + 4 ;
}
if (subid == StripSubdetector::TID) {
TIDDetId tidid(iidd);
return tidid.wheel() + 10;
}
if (subid == StripSubdetector::TEC) {
TECDetId tecid(iidd);
return tecid.wheel() + 13 ;
}
return 0;
}
| void HitEff::endJob | ( | void | ) | [private, virtual] |
Reimplemented from edm::EDAnalyzer.
Definition at line 669 of file HitEff.cc.
References gather_cfg::cout, events, EventTrackCKF, and traj.
| bool HitEff::isDoubleSided | ( | unsigned int | iidd | ) | const |
Definition at line 684 of file HitEff.cc.
References TIBDetId::layer(), TOBDetId::layer(), TIDDetId::ring(), TECDetId::ring(), strip(), DetId::subdetId(), StripSubdetector::TEC, StripSubdetector::TIB, StripSubdetector::TID, and StripSubdetector::TOB.
Referenced by analyze().
{
StripSubdetector strip=StripSubdetector(iidd);
unsigned int subid=strip.subdetId();
unsigned int layer = 0;
if (subid == StripSubdetector::TIB) {
TIBDetId tibid(iidd);
layer = tibid.layer();
if (layer == 1 || layer == 2) return true;
else return false;
}
else if (subid == StripSubdetector::TOB) {
TOBDetId tobid(iidd);
layer = tobid.layer() + 4 ;
if (layer == 5 || layer == 6) return true;
else return false;
}
else if (subid == StripSubdetector::TID) {
TIDDetId tidid(iidd);
layer = tidid.ring() + 10;
if (layer == 11 || layer == 12) return true;
else return false;
}
else if (subid == StripSubdetector::TEC) {
TECDetId tecid(iidd);
layer = tecid.ring() + 13 ;
if (layer == 14 || layer == 15 || layer == 18) return true;
else return false;
}
else
return false;
}
unsigned int HitEff::bunchx [private] |
Definition at line 77 of file HitEff.h.
Referenced by analyze(), and beginJob().
float HitEff::chi2 [private] |
Definition at line 76 of file HitEff.h.
Referenced by analyze(), and beginJob().
float HitEff::ClusterLocErrX [private] |
Definition at line 72 of file HitEff.h.
Referenced by analyze(), and beginJob().
float HitEff::ClusterLocErrY [private] |
Definition at line 72 of file HitEff.h.
Referenced by analyze(), and beginJob().
float HitEff::ClusterLocX [private] |
Definition at line 72 of file HitEff.h.
Referenced by analyze(), and beginJob().
float HitEff::ClusterLocY [private] |
Definition at line 72 of file HitEff.h.
Referenced by analyze(), and beginJob().
float HitEff::ClusterStoN [private] |
Definition at line 72 of file HitEff.h.
Referenced by analyze(), and beginJob().
edm::ParameterSet HitEff::conf_ [private] |
bool HitEff::DEBUG [private] |
float HitEff::dedx [private] |
Definition at line 80 of file HitEff.h.
Referenced by beginJob().
int HitEff::dedxNOM [private] |
Definition at line 81 of file HitEff.h.
Referenced by beginJob().
unsigned int HitEff::event [private] |
Definition at line 77 of file HitEff.h.
Referenced by beginJob().
int HitEff::events [private] |
Definition at line 62 of file HitEff.h.
Referenced by analyze(), beginJob(), and endJob().
int HitEff::EventTrackCKF [private] |
Definition at line 62 of file HitEff.h.
Referenced by analyze(), beginJob(), and endJob().
unsigned int HitEff::Id [private] |
Definition at line 74 of file HitEff.h.
Referenced by analyze(), and beginJob().
int HitEff::istep [private] |
Definition at line 83 of file HitEff.h.
Referenced by beginJob().
unsigned int HitEff::layers [private] |
unsigned int HitEff::ModIsBad [private] |
Definition at line 74 of file HitEff.h.
Referenced by analyze(), and beginJob().
int HitEff::nHits [private] |
Definition at line 75 of file HitEff.h.
Referenced by analyze(), and beginJob().
int HitEff::nLostHits [private] |
Definition at line 75 of file HitEff.h.
Referenced by analyze(), and beginJob().
float HitEff::p [private] |
Definition at line 76 of file HitEff.h.
Referenced by analyze(), and beginJob().
float HitEff::pT [private] |
Definition at line 76 of file HitEff.h.
Referenced by analyze(), and beginJob().
float HitEff::ResX [private] |
Definition at line 73 of file HitEff.h.
Referenced by analyze(), and beginJob().
float HitEff::ResXSig [private] |
Definition at line 73 of file HitEff.h.
Referenced by analyze(), and beginJob().
unsigned int HitEff::run [private] |
Definition at line 77 of file HitEff.h.
Referenced by analyze(), and beginJob().
unsigned int HitEff::SiStripQualBad [private] |
Definition at line 74 of file HitEff.h.
Referenced by analyze(), and beginJob().
float HitEff::timeDT [private] |
Definition at line 78 of file HitEff.h.
Referenced by analyze(), and beginJob().
int HitEff::timeDTDOF [private] |
Definition at line 79 of file HitEff.h.
Referenced by analyze(), and beginJob().
float HitEff::timeDTErr [private] |
Definition at line 78 of file HitEff.h.
Referenced by analyze(), and beginJob().
float HitEff::timeECAL [private] |
Definition at line 80 of file HitEff.h.
Referenced by analyze(), and beginJob().
int HitEff::tquality [private] |
Definition at line 82 of file HitEff.h.
Referenced by beginJob().
TTree* HitEff::traj [private] |
Definition at line 61 of file HitEff.h.
Referenced by analyze(), beginJob(), and endJob().
float HitEff::TrajGlbX [private] |
Definition at line 70 of file HitEff.h.
Referenced by analyze(), and beginJob().
float HitEff::TrajGlbY [private] |
Definition at line 70 of file HitEff.h.
Referenced by analyze(), and beginJob().
float HitEff::TrajGlbZ [private] |
Definition at line 70 of file HitEff.h.
Referenced by analyze(), and beginJob().
unsigned int HitEff::trajHitValid [private] |
Definition at line 77 of file HitEff.h.
Referenced by analyze(), and beginJob().
float HitEff::TrajLocAngleX [private] |
Definition at line 71 of file HitEff.h.
Referenced by analyze(), and beginJob().
float HitEff::TrajLocAngleY [private] |
Definition at line 71 of file HitEff.h.
Referenced by analyze(), and beginJob().
float HitEff::TrajLocErrX [private] |
Definition at line 71 of file HitEff.h.
Referenced by analyze(), and beginJob().
float HitEff::TrajLocErrY [private] |
Definition at line 71 of file HitEff.h.
Referenced by analyze(), and beginJob().
float HitEff::TrajLocX [private] |
Definition at line 71 of file HitEff.h.
Referenced by analyze(), and beginJob().
float HitEff::TrajLocY [private] |
Definition at line 71 of file HitEff.h.
Referenced by analyze(), and beginJob().
unsigned int HitEff::whatlayer [private] |
Definition at line 66 of file HitEff.h.
Referenced by analyze(), and beginJob().
bool HitEff::withinAcceptance [private] |
Definition at line 74 of file HitEff.h.
Referenced by analyze(), and beginJob().
1.7.3