00001 #include "RecoTracker/DebugTools/interface/TestTrackHits.h"
00002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00003 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00004 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00005 #include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
00006 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
00007 #include "SimTracker/Records/interface/TrackAssociatorRecord.h"
00008 #include <TDirectory.h>
00009
00010 #include "DataFormats/SiStripDetId/interface/TIBDetId.h"
00011 #include "DataFormats/SiStripDetId/interface/TOBDetId.h"
00012 #include "DataFormats/SiStripDetId/interface/TIDDetId.h"
00013 #include "DataFormats/SiStripDetId/interface/TECDetId.h"
00014 #include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
00015 #include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
00016 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00017
00018 typedef TrajectoryStateOnSurface TSOS;
00019 typedef TransientTrackingRecHit::ConstRecHitPointer CTTRHp;
00020 using namespace std;
00021 using namespace edm;
00022
00023 TestTrackHits::TestTrackHits(const edm::ParameterSet& iConfig):
00024 conf_(iConfig) {
00025 LogTrace("TestTrackHits") << conf_;
00026 propagatorName = conf_.getParameter<std::string>("Propagator");
00027 builderName = conf_.getParameter<std::string>("TTRHBuilder");
00028 srcName = conf_.getParameter<std::string>("src");
00029 tpName = conf_.getParameter<std::string>("tpname");
00030 updatorName = conf_.getParameter<std::string>("updator");
00031 out = conf_.getParameter<std::string>("out");
00032
00033
00034
00035
00036
00037
00038
00039
00040 }
00041
00042 TestTrackHits::~TestTrackHits(){}
00043
00044 void TestTrackHits::beginRun(edm::Run & run, const edm::EventSetup& iSetup)
00045 {
00046 iSetup.get<TrackerDigiGeometryRecord>().get(theG);
00047 iSetup.get<IdealMagneticFieldRecord>().get(theMF);
00048 iSetup.get<TrackingComponentsRecord>().get(propagatorName,thePropagator);
00049 iSetup.get<TransientRecHitRecord>().get(builderName,theBuilder);
00050 iSetup.get<TrackingComponentsRecord>().get(updatorName,theUpdator);
00051 iSetup.get<TrackAssociatorRecord>().get("TrackAssociatorByHits",trackAssociator);
00052
00053 file = new TFile(out.c_str(),"recreate");
00054 for (int i=0; i!=6; i++)
00055 for (int j=0; j!=9; j++){
00056 if (i==0 && j>2) break;
00057 if (i==1 && j>1) break;
00058 if (i==2 && j>3) break;
00059 if (i==3 && j>2) break;
00060 if (i==4 && j>5) break;
00061 if (i==5 && j>8) break;
00062
00063 title.str("");
00064 title << "Chi2Increment_" << i+1 << "-" << j+1 ;
00065 hChi2Increment[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
00066 title.str("");
00067 title << "Chi2IncrementVsEta_" << i+1 << "-" << j+1 ;
00068 hChi2IncrementVsEta[title.str()] = new TH2F(title.str().c_str(),title.str().c_str(),50,-2.5,2.5,1000,0,100);
00069 title.str("");
00070 title << "Chi2GoodHit_" << i+1 << "-" << j+1 ;
00071 hChi2GoodHit[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
00072 title.str("");
00073 title << "Chi2BadHit_" << i+1 << "-" << j+1 ;
00074 hChi2BadHit[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
00075 title.str("");
00076 title << "Chi2DeltaHit_" << i+1 << "-" << j+1 ;
00077 hChi2DeltaHit[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
00078 title.str("");
00079 title << "Chi2NSharedHit_" << i+1 << "-" << j+1 ;
00080 hChi2NSharedHit[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
00081 title.str("");
00082 title << "Chi2SharedHit_" << i+1 << "-" << j+1 ;
00083 hChi2SharedHit[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
00084
00085 title.str("");
00086 title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts";
00087 hPullGP_X_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00088 title.str("");
00089 title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts";
00090 hPullGP_Y_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00091 title.str("");
00092 title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts";
00093 hPullGP_Z_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00094
00095 title.str("");
00096 title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts";
00097 hPullGM_X_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00098 title.str("");
00099 title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts";
00100 hPullGM_Y_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00101 title.str("");
00102 title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts";
00103 hPullGM_Z_ts[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00104
00105 title.str("");
00106 title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr";
00107 hPullGP_X_tr[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00108 title.str("");
00109 title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr";
00110 hPullGP_Y_tr[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00111 title.str("");
00112 title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr";
00113 hPullGP_Z_tr[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00114
00115 title.str("");
00116 title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs";
00117 hPullGP_X_rs[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00118 title.str("");
00119 title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs";
00120 hPullGP_Y_rs[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00121 title.str("");
00122 title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs";
00123 hPullGP_Z_rs[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00124
00125 if ( ((i==2||i==4)&&(j==0||j==1)) || (i==3||i==5) ){
00126
00127 title.str("");
00128 title << "Chi2Increment_mono_" << i+1 << "-" << j+1 ;
00129 hChi2Increment_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
00130
00131 title.str("");
00132 title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_mono";
00133 hPullGP_X_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00134 title.str("");
00135 title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts_mono";
00136 hPullGP_Y_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00137 title.str("");
00138 title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts_mono";
00139 hPullGP_Z_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00140
00141 title.str("");
00142 title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_mono";
00143 hPullGM_X_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00144 title.str("");
00145 title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_mono";
00146 hPullGM_Y_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00147 title.str("");
00148 title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_mono";
00149 hPullGM_Z_ts_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00150
00151 title.str("");
00152 title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_mono";
00153 hPullGP_X_tr_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00154 title.str("");
00155 title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_mono";
00156 hPullGP_Y_tr_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00157 title.str("");
00158 title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_mono";
00159 hPullGP_Z_tr_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00160
00161 title.str("");
00162 title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_mono";
00163 hPullGP_X_rs_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00164 title.str("");
00165 title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs_mono";
00166 hPullGP_Y_rs_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00167 title.str("");
00168 title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs_mono";
00169 hPullGP_Z_rs_mono[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00170
00171
00172 title.str("");
00173 title << "Chi2Increment_stereo_" << i+1 << "-" << j+1 ;
00174 hChi2Increment_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,0,100);
00175
00176 title.str("");
00177 title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_stereo";
00178 hPullGP_X_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00179 title.str("");
00180 title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
00181 hPullGP_Y_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00182 title.str("");
00183 title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
00184 hPullGP_Z_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00185
00186 title.str("");
00187 title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_stereo";
00188 hPullGM_X_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00189 title.str("");
00190 title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
00191 hPullGM_Y_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00192 title.str("");
00193 title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
00194 hPullGM_Z_ts_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00195
00196 title.str("");
00197 title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_stereo";
00198 hPullGP_X_tr_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00199 title.str("");
00200 title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_stereo";
00201 hPullGP_Y_tr_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00202 title.str("");
00203 title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_stereo";
00204 hPullGP_Z_tr_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00205
00206 title.str("");
00207 title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_stereo";
00208 hPullGP_X_rs_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00209 title.str("");
00210 title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs_stereo";
00211 hPullGP_Y_rs_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00212 title.str("");
00213 title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs_stereo";
00214 hPullGP_Z_rs_stereo[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00215 }
00216 }
00217 hTotChi2Increment = new TH1F("TotChi2Increment","TotChi2Increment",1000,0,100);
00218 hTotChi2GoodHit = new TH1F("TotChi2GoodHit","TotChi2GoodHit",1000,0,100);
00219 hTotChi2BadHit = new TH1F("TotChi2BadHit","TotChi2BadHit",1000,0,100);
00220 hTotChi2DeltaHit = new TH1F("TotChi2DeltaHit","TotChi2DeltaHit",1000,0,100);
00221 hTotChi2NSharedHit = new TH1F("TotChi2NSharedHit","TotChi2NSharedHit",1000,0,100);
00222 hTotChi2SharedHit = new TH1F("TotChi2SharedHit","TotChi2SharedHit",1000,0,100);
00223 hProcess_vs_Chi2 = new TH2F("Process_vs_Chi2","Process_vs_Chi2",1000,0,100,17,-0.5,16.5);
00224 hClsize_vs_Chi2 = new TH2F("Clsize_vs_Chi2","Clsize_vs_Chi2",1000,0,100,17,-0.5,16.5);
00225 hPixClsize_vs_Chi2= new TH2F("PixClsize_vs_Chi2","PixClsize_vs_Chi2",1000,0,100,17,-0.5,16.5);
00226 hPrjClsize_vs_Chi2= new TH2F("PrjClsize_vs_Chi2","PrjClsize_vs_Chi2",1000,0,100,17,-0.5,16.5);
00227 hSt1Clsize_vs_Chi2= new TH2F("St1Clsize_vs_Chi2","St1Clsize_vs_Chi2",1000,0,100,17,-0.5,16.5);
00228 hSt2Clsize_vs_Chi2= new TH2F("St2Clsize_vs_Chi2","St2Clsize_vs_Chi2",1000,0,100,17,-0.5,16.5);
00229 hGoodHit_vs_Chi2 = new TH2F("GoodHit_vs_Chi2","GoodHit_vs_Chi2",10000,0,1000,2,-0.5,1.5);
00230 hClusterSize = new TH1F("ClusterSize","ClusterSize",40,-0.5,39.5);
00231 hPixClusterSize = new TH1F("PixClusterSize","PixClusterSize",40,-0.5,39.5);
00232 hPrjClusterSize = new TH1F("PrjClusterSize","PrjClusterSize",40,-0.5,39.5);
00233 hSt1ClusterSize = new TH1F("St1ClusterSize","St1ClusterSize",40,-0.5,39.5);
00234 hSt2ClusterSize = new TH1F("St2ClusterSize","St2ClusterSize",40,-0.5,39.5);
00235 hSimHitVecSize = new TH1F("hSimHitVecSize","hSimHitVecSize",40,-0.5,39.5);
00236 hPixSimHitVecSize = new TH1F("PixSimHitVecSize","PixSimHitVecSize",40,-0.5,39.5);
00237 hPrjSimHitVecSize = new TH1F("PrjSimHitVecSize","PrjSimHitVecSize",40,-0.5,39.5);
00238 hSt1SimHitVecSize = new TH1F("St1SimHitVecSize","St1SimHitVecSize",40,-0.5,39.5);
00239 hSt2SimHitVecSize = new TH1F("St2SimHitVecSize","St2SimHitVecSize",40,-0.5,39.5);
00240 goodbadmerged = new TH1F("goodbadmerged","goodbadmerged",5,0.5,5.5);
00241 energyLossRatio = new TH1F("energyLossRatio","energyLossRatio",100,0,1);
00242 mergedPull = new TH1F("mergedPull","mergedPull",200,0,20);
00243 probXgood = new TH1F("probXgood","probXgood",110,0,1.1);
00244 probXbad = new TH1F("probXbad","probXbad",110,0,1.1);
00245 probXdelta = new TH1F("probXdelta","probXdelta",110,0,1.1);
00246 probXshared = new TH1F("probXshared","probXshared",110,0,1.1);
00247 probXnoshare = new TH1F("probXnoshare","probXnoshare",110,0,1.1);
00248 probYgood = new TH1F("probYgood","probYgood",110,0,1.1);
00249 probYbad = new TH1F("probYbad","probYbad",110,0,1.1);
00250 probYdelta = new TH1F("probYdelta","probYdelta",110,0,1.1);
00251 probYshared = new TH1F("probYshared","probYshared",110,0,1.1);
00252 probYnoshare = new TH1F("probYnoshare","probYnoshare",110,0,1.1);
00253 }
00254
00255 void TestTrackHits::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00256 {
00257 LogDebug("TestTrackHits") << "new event" ;
00258 iEvent.getByLabel(srcName,trajCollectionHandle);
00259 iEvent.getByLabel(srcName,trackCollectionHandle);
00260 iEvent.getByLabel(srcName,trajTrackAssociationCollectionHandle);
00261 iEvent.getByLabel(tpName,trackingParticleCollectionHandle);
00262 edm::Handle<reco::BeamSpot> beamSpot;
00263 iEvent.getByLabel("offlineBeamSpot",beamSpot);
00264
00265 LogTrace("TestTrackHits") << "Tr collection size=" << trackCollectionHandle->size();
00266 LogTrace("TestTrackHits") << "TP collection size=" << trackingParticleCollectionHandle->size();
00267
00268 hitAssociator = new TrackerHitAssociator(iEvent);
00269
00270 reco::RecoToSimCollection recSimColl=trackAssociator->associateRecoToSim(trackCollectionHandle,
00271 trackingParticleCollectionHandle,
00272 &iEvent);
00273
00274 TrajectoryStateCombiner combiner;
00275
00276 int evtHits = 0;
00277 int i=0;
00278 int yy=0;
00279 int yyy=0;
00280 for(std::vector<Trajectory>::const_iterator it = trajCollectionHandle->begin(); it!=trajCollectionHandle->end();it++){
00281
00282 LogTrace("TestTrackHits") << "\n*****************new trajectory********************" ;
00283 double tchi2 = 0;
00284
00285 std::vector<TrajectoryMeasurement> tmColl = it->measurements();
00286
00287 edm::Ref<std::vector<Trajectory> > traj(trajCollectionHandle, i);
00288 reco::TrackRef tmptrack = (*trajTrackAssociationCollectionHandle.product())[traj];
00289 edm::RefToBase<reco::Track> track(tmptrack);
00290
00291
00292
00293
00294
00295
00296 std::vector<std::pair<TrackingParticleRef, double> > tP;
00297 if(recSimColl.find(track) != recSimColl.end()){
00298 tP = recSimColl[track];
00299 if (tP.size()!=0) {
00300 edm::LogVerbatim("TestTrackHits") << "reco::Track #" << ++yyy << " with pt=" << track->pt()
00301 << " associated with quality:" << tP.begin()->second <<" good track #" << ++yy << " has hits:" << track->numberOfValidHits() << "\n";
00302 }
00303 }else{
00304 edm::LogVerbatim("TestTrackHits") << "reco::Track #" << ++yyy << " with pt=" << track->pt()
00305 << " NOT associated to any TrackingParticle" << "\n";
00306 i++;
00307 continue;
00308 }
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319 TrackingParticleRef tp = tP.begin()->first;
00320 LogTrace("TestTrackHits") << "a tp is associated with fraction=" << tP.begin()->second;
00321
00322 std::vector<unsigned int> tpids;
00323 for (TrackingParticle::g4t_iterator g4T=tp->g4Track_begin(); g4T!=tp->g4Track_end(); ++g4T) {
00324 LogTrace("TestTrackHits") << "tp id=" << g4T->trackId();
00325 tpids.push_back(g4T->trackId());
00326 }
00327
00328
00329 int pp = 0;
00330 for (std::vector<TrajectoryMeasurement>::iterator tm=tmColl.begin();tm!=tmColl.end();++tm){
00331
00332 tchi2+=tm->estimate();
00333
00334 LogTrace("TestTrackHits") << "+++++++++++++++++new hit+++++++++++++++++" ;
00335 CTTRHp rhit = tm->recHit();
00336
00337
00338 TSOS state = combiner(tm->backwardPredictedState(), tm->forwardPredictedState());
00339
00340 if (rhit->isValid()==0 && rhit->det()!=0) continue;
00341 evtHits++;
00342 LogTrace("TestTrackHits") << "valid hit #" << ++pp << "of hits=" << track->numberOfValidHits();
00343
00344 int subdetId = rhit->det()->geographicalId().subdetId();
00345 int layerId = 0;
00346 DetId id = rhit->det()->geographicalId();
00347 if (id.subdetId()==3) layerId = ((TIBDetId)(id)).layer();
00348 if (id.subdetId()==5) layerId = ((TOBDetId)(id)).layer();
00349 if (id.subdetId()==1) layerId = ((PXBDetId)(id)).layer();
00350 if (id.subdetId()==4) layerId = ((TIDDetId)(id)).wheel();
00351 if (id.subdetId()==6) layerId = ((TECDetId)(id)).wheel();
00352 if (id.subdetId()==2) layerId = ((PXFDetId)(id)).disk();
00353 LogTrace("TestTrackHits") << "subdetId=" << subdetId << " layerId=" << layerId ;
00354
00355 const Surface * surf = rhit->surface();
00356 if (surf==0) continue;
00357
00358 double energyLoss_ = 0.;
00359 unsigned int monoId = 0;
00360 std::vector<double> energyLossM;
00361 std::vector<double> energyLossS;
00362 std::vector<PSimHit> assSimHits = hitAssociator->associateHit(*(rhit)->hit());
00363 unsigned int simhitvecsize = assSimHits.size();
00364 if (simhitvecsize==0) continue;
00365 PSimHit shit;
00366 std::vector<unsigned int> trackIds;
00367 energyLossS.clear();
00368 energyLossM.clear();
00369 for(std::vector<PSimHit>::const_iterator m=assSimHits.begin(); m<assSimHits.end(); m++){
00370 unsigned int tId = m->trackId();
00371 if (find(trackIds.begin(),trackIds.end(),tId)==trackIds.end()) trackIds.push_back(tId);
00372 if (m->energyLoss()>energyLoss_) {
00373 shit=*m;
00374 energyLoss_ = m->energyLoss();
00375 }
00376 if (dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())) {
00377 if (monoId==0) monoId = m->detUnitId();
00378 if (monoId == m->detUnitId()){
00379 energyLossM.push_back(m->energyLoss());
00380 }
00381 else {
00382 energyLossS.push_back(m->energyLoss());
00383 }
00384
00385 } else {
00386 energyLossM.push_back(m->energyLoss());
00387 }
00388 }
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403 double chi2increment = tm->estimate();
00404 LogTrace("TestTrackHits") << "tm->estimate()=" << tm->estimate();
00405 title.str("");
00406 title << "Chi2Increment_" << subdetId << "-" << layerId;
00407 hChi2Increment[title.str()]->Fill( chi2increment );
00408 title.str("");
00409 title << "Chi2IncrementVsEta_" << subdetId << "-" << layerId;
00410 hChi2IncrementVsEta[title.str()]->Fill( track->eta(), chi2increment );
00411 hTotChi2Increment->Fill( chi2increment );
00412 hProcess_vs_Chi2->Fill( chi2increment, shit.processType() );
00413
00414 int clustersize = 0;
00415 bool mergedhit = false;
00416 if (dynamic_cast<const SiPixelRecHit*>(rhit->hit())){
00417 clustersize = ((const SiPixelRecHit*)(rhit->hit()))->cluster()->size() ;
00418 hPixClsize_vs_Chi2->Fill(chi2increment, clustersize);
00419 hPixClusterSize->Fill(clustersize);
00420 hPixSimHitVecSize->Fill(simhitvecsize);
00421 if (simhitvecsize>1) mergedhit = true;
00422 }
00423 else if (dynamic_cast<const SiStripRecHit2D*>(rhit->hit())){
00424 clustersize = ((const SiStripRecHit2D*)(rhit->hit()))->cluster()->amplitudes().size() ;
00425 hSt1Clsize_vs_Chi2->Fill(chi2increment, clustersize);
00426 hSt1ClusterSize->Fill(clustersize);
00427 hSt1SimHitVecSize->Fill(simhitvecsize);
00428 if (simhitvecsize>1) mergedhit = true;
00429 }
00430 else if (dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())){
00431 int clsize1 = ((const SiStripMatchedRecHit2D*)(rhit->hit()))->monoCluster().amplitudes().size();
00432 int clsize2 = ((const SiStripMatchedRecHit2D*)(rhit->hit()))->stereoCluster().amplitudes().size();
00433 if (clsize1>clsize2) clustersize = clsize1;
00434 else clustersize = clsize2;
00435 hSt2Clsize_vs_Chi2->Fill(chi2increment, clustersize);
00436 hSt2ClusterSize->Fill(clustersize);
00437 hSt2SimHitVecSize->Fill(simhitvecsize);
00438 if (simhitvecsize>2) mergedhit = true;
00439 }
00440 else if (dynamic_cast<const ProjectedSiStripRecHit2D*>(rhit->hit())){
00441 clustersize = ((const ProjectedSiStripRecHit2D*)(rhit->hit()))->originalHit().cluster()->amplitudes().size();
00442 hPrjClsize_vs_Chi2->Fill(chi2increment, clustersize);
00443 hPrjClusterSize->Fill(clustersize);
00444 hPrjSimHitVecSize->Fill(simhitvecsize);
00445 if (simhitvecsize>1) mergedhit = true;
00446 }
00447 hClsize_vs_Chi2->Fill( chi2increment, clustersize);
00448 hClusterSize->Fill(clustersize);
00449 hSimHitVecSize->Fill(simhitvecsize);
00450
00451
00452
00453
00454
00455
00456 std::vector<SimHitIdpr> simTrackIds = hitAssociator->associateHitId(*(rhit)->hit());
00457 bool goodhit = false;
00458 for(size_t j=0; j<simTrackIds.size(); j++){
00459 LogTrace("TestTrackHits") << "hit id=" << simTrackIds[j].first;
00460 for (size_t jj=0; jj<tpids.size(); jj++){
00461 if (simTrackIds[j].first == tpids[jj]) goodhit = true;
00462 break;
00463 }
00464 }
00465 bool shared = true;
00466 bool ioniOnly = true;
00467 const SiPixelRecHit* pix = dynamic_cast<const SiPixelRecHit*>(rhit->hit());
00468 if (goodhit) {
00469 if (energyLossM.size()>1&&energyLossS.size()<=1) {
00470 sort(energyLossM.begin(),energyLossM.end(),greater<double>());
00471 energyLossRatio->Fill(energyLossM[1]/energyLossM[0]);
00472 }
00473 else if (energyLossS.size()>1&&energyLossM.size()<=1) {
00474 sort(energyLossS.begin(),energyLossS.end(),greater<double>());
00475 energyLossRatio->Fill(energyLossS[1]/energyLossS[0]);
00476 }
00477 else if (energyLossS.size()>1&&energyLossM.size()>1) {
00478 sort(energyLossM.begin(),energyLossM.end(),greater<double>());
00479 sort(energyLossS.begin(),energyLossS.end(),greater<double>());
00480 energyLossM[1]/energyLossM[0] > energyLossS[1]/energyLossS[0]
00481 ? energyLossRatio->Fill(energyLossM[1]/energyLossM[0])
00482 : energyLossRatio->Fill(energyLossS[1]/energyLossS[0]);
00483 }
00484 if ( mergedhit ) {
00485
00486 LogVerbatim("TestTrackHits") << "MERGED HIT" << std::endl;
00487 unsigned int idc = 0;
00488 for (size_t jj=0; jj<tpids.size(); jj++){
00489 idc += std::count(trackIds.begin(),trackIds.end(),tpids[jj]);
00490 }
00491 if (idc==trackIds.size()) {
00492 shared = false;
00493 }
00494 for(std::vector<PSimHit>::const_iterator m=assSimHits.begin()+1; m<assSimHits.end(); m++){
00495 if ((m->processType()!=7&&m->processType()!=8&&m->processType()!=9)&&abs(m->particleType())!=11){
00496 ioniOnly = false;
00497 break;
00498 }
00499 }
00500 if (ioniOnly&&!shared) {
00501 title.str("");
00502 title << "Chi2DeltaHit_" << subdetId << "-" << layerId;
00503 hChi2DeltaHit[title.str()]->Fill( chi2increment );
00504 hTotChi2DeltaHit->Fill( chi2increment );
00505 if (pix) {
00506 probXdelta->Fill(pix->probabilityX());
00507 probYdelta->Fill(pix->probabilityY());
00508 }
00509 } else if(!ioniOnly&&!shared) {
00510 title.str("");
00511 title << "Chi2NSharedHit_" << subdetId << "-" << layerId;
00512 hChi2NSharedHit[title.str()]->Fill( chi2increment );
00513 hTotChi2NSharedHit->Fill( chi2increment );
00514 if (pix) {
00515 probXnoshare->Fill(pix->probabilityX());
00516 probYnoshare->Fill(pix->probabilityY());
00517 }
00518 } else {
00519 title.str("");
00520 title << "Chi2SharedHit_" << subdetId << "-" << layerId;
00521 hChi2SharedHit[title.str()]->Fill( chi2increment );
00522 hTotChi2SharedHit->Fill( chi2increment );
00523 if (pix) {
00524 probXshared->Fill(pix->probabilityX());
00525 probYshared->Fill(pix->probabilityY());
00526 }
00527 }
00528
00529 for(std::vector<PSimHit>::const_iterator m=assSimHits.begin(); m<assSimHits.end(); m++) {
00530 unsigned int tId = m->trackId();
00531 LogVerbatim("TestTrackHits") << "component with id=" << tId << " eLoss=" << m->energyLoss() << " pType=" << m->processType();
00532 if (find(tpids.begin(),tpids.end(),tId)==tpids.end()) continue;
00533 if (m->processType()==2) {
00534 GlobalPoint gpr = rhit->globalPosition();
00535 AlgebraicSymMatrix33 ger = rhit->globalPositionError().matrix();
00536 GlobalPoint gps = surf->toGlobal(m->localPosition());
00537 LogVerbatim("TestTrackHits") << gpr << " " << gps << " " << ger;
00538 ROOT::Math::SVector<double,3> delta;
00539 delta[0]=gpr.x()-gps.x();
00540 delta[1]=gpr.y()-gps.y();
00541 delta[2]=gpr.z()-gps.z();
00542 LogVerbatim("TestTrackHits") << delta << " " << ger ;
00543 double mpull = sqrt(delta[0]*delta[0]/ger[0][0]+delta[1]*delta[1]/ger[1][1]+delta[2]*delta[2]/ger[2][2]);
00544 LogVerbatim("TestTrackHits") << "hit pull=" << mpull;
00545 mergedPull->Fill(mpull);
00546 break;
00547 }
00548 }
00549 } else {
00550 LogVerbatim("TestTrackHits") << "good hit" ;
00551 title.str("");
00552 title << "Chi2GoodHit_" << subdetId << "-" << layerId;
00553 hChi2GoodHit[title.str()]->Fill( chi2increment );
00554 hTotChi2GoodHit->Fill( chi2increment );
00555 if (pix) {
00556 probXgood->Fill(pix->probabilityX());
00557 probYgood->Fill(pix->probabilityY());
00558 }
00559 }
00560 } else {
00561 LogVerbatim("TestTrackHits") << "bad hit" ;
00562 title.str("");
00563 title << "Chi2BadHit_" << subdetId << "-" << layerId;
00564 hChi2BadHit[title.str()]->Fill( chi2increment );
00565 hTotChi2BadHit->Fill( chi2increment );
00566 goodbadmerged->Fill(2);
00567 if (pix) {
00568 probXbad->Fill(pix->probabilityX());
00569 probYbad->Fill(pix->probabilityY());
00570 }
00571 }
00572 hGoodHit_vs_Chi2->Fill(chi2increment,goodhit);
00573
00574 LocalVector shitLMom;
00575 LocalPoint shitLPos;
00576 if (dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())) {
00577 if (simhitvecsize>2 && goodhit) {
00578 if (ioniOnly&&!shared) goodbadmerged->Fill(3);
00579 else if(!ioniOnly&&!shared) goodbadmerged->Fill(4);
00580 else goodbadmerged->Fill(5);
00581 }
00582 else if (goodhit) goodbadmerged->Fill(1);
00583 double rechitmatchedx = rhit->localPosition().x();
00584 double rechitmatchedy = rhit->localPosition().y();
00585 double mindist = 999999;
00586 float distx, disty;
00587 std::pair<LocalPoint,LocalVector> closestPair;
00588 const StripGeomDetUnit* stripDet =(StripGeomDetUnit*) ((const GluedGeomDet *)(rhit)->det())->stereoDet();
00589 const BoundPlane& plane = (rhit)->det()->surface();
00590 for(std::vector<PSimHit>::const_iterator m=assSimHits.begin(); m<assSimHits.end(); m++){
00591
00592 std::pair<LocalPoint,LocalVector> hitPair = projectHit((*m),stripDet,plane);
00593 distx = fabs(rechitmatchedx - hitPair.first.x());
00594 disty = fabs(rechitmatchedy - hitPair.first.y());
00595 double dist = distx*distx+disty*disty;
00596 if(sqrt(dist)<mindist){
00597 mindist = dist;
00598 closestPair = hitPair;
00599 }
00600 }
00601 shitLPos = closestPair.first;
00602 shitLMom = closestPair.second;
00603 } else {
00604 if (simhitvecsize>1 && goodhit) {
00605 if (ioniOnly&&!shared) goodbadmerged->Fill(3);
00606 else if(!ioniOnly&&!shared) goodbadmerged->Fill(4);
00607 else goodbadmerged->Fill(5);
00608 }
00609 else if (goodhit) goodbadmerged->Fill(1);
00610 shitLPos = shit.localPosition();
00611 shitLMom = shit.momentumAtEntry();
00612 }
00613 GlobalVector shitGMom = surf->toGlobal(shitLMom);
00614 GlobalPoint shitGPos = surf->toGlobal(shitLPos);
00615
00616 GlobalVector tsosGMom = state.globalMomentum();
00617 GlobalError tsosGMEr(state.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3,3));
00618 GlobalPoint tsosGPos = state.globalPosition();
00619 GlobalError tsosGPEr = state.cartesianError().position();
00620
00621 GlobalPoint rhitGPos = (rhit)->globalPosition();
00622 GlobalError rhitGPEr = (rhit)->globalPositionError();
00623
00624 LogVerbatim("TestTrackHits") << "assSimHits.size()=" << assSimHits.size() ;
00625 LogVerbatim("TestTrackHits") << "tsos globalPos =" << tsosGPos ;
00626 LogVerbatim("TestTrackHits") << "sim hit globalPos=" << shitGPos ;
00627 LogVerbatim("TestTrackHits") << "rec hit globalPos=" << rhitGPos ;
00628 LogVerbatim("TestTrackHits") << "geographicalId =" << rhit->det()->geographicalId().rawId() ;
00629 LogVerbatim("TestTrackHits") << "surface position =" << surf->position() ;
00630
00631 # if 0
00632 if (rhit->detUnit()) LogTrace("TestTrackHits") << "rhit->detUnit()->geographicalId()="
00633 << rhit->detUnit()->geographicalId().rawId() ;
00634 LogTrace("TestTrackHits") << "rhit->det()->surface().position()="
00635 << rhit->det()->surface().position() ;
00636 if (rhit->detUnit()) LogTrace("TestTrackHits") << "rhit->detUnit()->surface().position()="
00637 << rhit->detUnit()->surface().position() ;
00638 LogTrace("TestTrackHits") << "rhit->det()->position()=" << rhit->det()->position() ;
00639 if (rhit->detUnit()) LogTrace("TestTrackHits") << "rhit->detUnit()->position()="
00640 << rhit->detUnit()->position() ;
00641 LogTrace("TestTrackHits") << "rhit->det()->surface().bounds().length()="
00642 << rhit->det()->surface().bounds().length() ;
00643 if (rhit->detUnit()) LogTrace("TestTrackHits") << "rhit->detUnit()->surface().bounds().length()="
00644 << rhit->detUnit()->surface().bounds().length() ;
00645 LogTrace("TestTrackHits") << "rhit->det()->surface().bounds().width()="
00646 << rhit->det()->surface().bounds().width() ;
00647 if (rhit->detUnit()) LogTrace("TestTrackHits") << "rhit->detUnit()->surface().bounds().width()="
00648 << rhit->detUnit()->surface().bounds().width() ;
00649 LogTrace("TestTrackHits") << "rhit->det()->surface().bounds().thickness()="
00650 << rhit->det()->surface().bounds().thickness() ;
00651 if (rhit->detUnit()) LogTrace("TestTrackHits") << "rhit->detUnit()->surface().bounds().thickness()="
00652 << rhit->detUnit()->surface().bounds().thickness() ;
00653 #endif
00654
00655 double pullGPX_rs = (rhitGPos.x()-shitGPos.x())/sqrt(rhitGPEr.cxx());
00656 double pullGPY_rs = (rhitGPos.y()-shitGPos.y())/sqrt(rhitGPEr.cyy());
00657 double pullGPZ_rs = (rhitGPos.z()-shitGPos.z())/sqrt(rhitGPEr.czz());
00658
00659
00660
00661
00662 LogTrace("TestTrackHits") << "rs" ;
00663
00664 title.str("");
00665 title << "PullGP_X_" << subdetId << "-" << layerId << "_rs";
00666 hPullGP_X_rs[title.str()]->Fill( pullGPX_rs );
00667 title.str("");
00668 title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs";
00669 hPullGP_Y_rs[title.str()]->Fill( pullGPY_rs );
00670 title.str("");
00671 title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs";
00672 hPullGP_Z_rs[title.str()]->Fill( pullGPZ_rs );
00673
00674 double pullGPX_tr = (tsosGPos.x()-rhitGPos.x())/sqrt(tsosGPEr.cxx()+rhitGPEr.cxx());
00675 double pullGPY_tr = (tsosGPos.y()-rhitGPos.y())/sqrt(tsosGPEr.cyy()+rhitGPEr.cyy());
00676 double pullGPZ_tr = (tsosGPos.z()-rhitGPos.z())/sqrt(tsosGPEr.czz()+rhitGPEr.czz());
00677
00678
00679
00680
00681 LogTrace("TestTrackHits") << "tr" ;
00682
00683 title.str("");
00684 title << "PullGP_X_" << subdetId << "-" << layerId << "_tr";
00685 hPullGP_X_tr[title.str()]->Fill( pullGPX_tr );
00686 title.str("");
00687 title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr";
00688 hPullGP_Y_tr[title.str()]->Fill( pullGPY_tr );
00689 title.str("");
00690 title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr";
00691 hPullGP_Z_tr[title.str()]->Fill( pullGPZ_tr );
00692
00693 double pullGPX_ts = (tsosGPos.x()-shitGPos.x())/sqrt(tsosGPEr.cxx());
00694 double pullGPY_ts = (tsosGPos.y()-shitGPos.y())/sqrt(tsosGPEr.cyy());
00695 double pullGPZ_ts = (tsosGPos.z()-shitGPos.z())/sqrt(tsosGPEr.czz());
00696
00697
00698
00699
00700 LogTrace("TestTrackHits") << "ts1" ;
00701
00702 title.str("");
00703 title << "PullGP_X_" << subdetId << "-" << layerId << "_ts";
00704 hPullGP_X_ts[title.str()]->Fill( pullGPX_ts );
00705 title.str("");
00706 title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts";
00707 hPullGP_Y_ts[title.str()]->Fill( pullGPY_ts );
00708 title.str("");
00709 title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts";
00710 hPullGP_Z_ts[title.str()]->Fill( pullGPZ_ts );
00711
00712 double pullGMX_ts = (tsosGMom.x()-shitGMom.x())/sqrt(tsosGMEr.cxx());
00713 double pullGMY_ts = (tsosGMom.y()-shitGMom.y())/sqrt(tsosGMEr.cyy());
00714 double pullGMZ_ts = (tsosGMom.z()-shitGMom.z())/sqrt(tsosGMEr.czz());
00715
00716
00717
00718
00719 LogTrace("TestTrackHits") << "ts2" ;
00720
00721 title.str("");
00722 title << "PullGM_X_" << subdetId << "-" << layerId << "_ts";
00723 hPullGM_X_ts[title.str()]->Fill( pullGMX_ts );
00724 title.str("");
00725 title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts";
00726 hPullGM_Y_ts[title.str()]->Fill( pullGMY_ts );
00727 title.str("");
00728 title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts";
00729 hPullGM_Z_ts[title.str()]->Fill( pullGMZ_ts );
00730 if (dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())) {
00731 Propagator* thePropagatorAnyDir = new PropagatorWithMaterial(anyDirection,0.105,theMF.product(),1.6);
00732
00733 LogTrace("TestTrackHits") << "MONO HIT" ;
00734 auto m = dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())->monoHit();
00735 CTTRHp tMonoHit = theBuilder->build(&m);
00736 if (tMonoHit==0) continue;
00737 vector<PSimHit> assMonoSimHits = hitAssociator->associateHit(*tMonoHit->hit());
00738 if (assMonoSimHits.size()==0) continue;
00739 const PSimHit sMonoHit = *(assMonoSimHits.begin());
00740 const Surface * monoSurf = &( tMonoHit->det()->surface() );
00741 if (monoSurf==0) continue;
00742 TSOS monoState = thePropagatorAnyDir->propagate(state,*monoSurf);
00743 if (monoState.isValid()==0) continue;
00744
00745 LocalVector monoShitLMom = sMonoHit.momentumAtEntry();
00746 GlobalVector monoShitGMom = monoSurf->toGlobal(monoShitLMom);
00747 LocalPoint monoShitLPos = sMonoHit.localPosition();
00748 GlobalPoint monoShitGPos = monoSurf->toGlobal(monoShitLPos);
00749
00750
00751
00752
00753 GlobalVector monoTsosGMom = monoState.globalMomentum();
00754 GlobalError monoTsosGMEr(monoState.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3,3));
00755 GlobalPoint monoTsosGPos = monoState.globalPosition();
00756 GlobalError monoTsosGPEr = monoState.cartesianError().position();
00757
00758 GlobalPoint monoRhitGPos = tMonoHit->globalPosition();
00759 GlobalError monoRhitGPEr = tMonoHit->globalPositionError();
00760
00761 double pullGPX_rs_mono = (monoRhitGPos.x()-monoShitGPos.x())/sqrt(monoRhitGPEr.cxx());
00762 double pullGPY_rs_mono = (monoRhitGPos.y()-monoShitGPos.y())/sqrt(monoRhitGPEr.cyy());
00763 double pullGPZ_rs_mono = (monoRhitGPos.z()-monoShitGPos.z())/sqrt(monoRhitGPEr.czz());
00764
00765
00766
00767
00768 MeasurementExtractor meMo(monoState);
00769 double chi2mono = computeChi2Increment(meMo,tMonoHit);
00770
00771 title.str("");
00772 title << "Chi2Increment_mono_" << subdetId << "-" << layerId ;
00773 hChi2Increment_mono[title.str()]->Fill(chi2mono);
00774
00775 title.str("");
00776 title << "PullGP_X_" << subdetId << "-" << layerId << "_rs_mono";
00777 hPullGP_X_rs_mono[title.str()]->Fill( pullGPX_rs_mono );
00778 title.str("");
00779 title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs_mono";
00780 hPullGP_Y_rs_mono[title.str()]->Fill( pullGPY_rs_mono );
00781 title.str("");
00782 title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs_mono";
00783 hPullGP_Z_rs_mono[title.str()]->Fill( pullGPZ_rs_mono );
00784
00785 double pullGPX_tr_mono = (monoTsosGPos.x()-monoRhitGPos.x())/sqrt(monoTsosGPEr.cxx()+monoRhitGPEr.cxx());
00786 double pullGPY_tr_mono = (monoTsosGPos.y()-monoRhitGPos.y())/sqrt(monoTsosGPEr.cyy()+monoRhitGPEr.cyy());
00787 double pullGPZ_tr_mono = (monoTsosGPos.z()-monoRhitGPos.z())/sqrt(monoTsosGPEr.czz()+monoRhitGPEr.czz());
00788
00789
00790
00791
00792 title.str("");
00793 title << "PullGP_X_" << subdetId << "-" << layerId << "_tr_mono";
00794 hPullGP_X_tr_mono[title.str()]->Fill( pullGPX_tr_mono );
00795 title.str("");
00796 title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr_mono";
00797 hPullGP_Y_tr_mono[title.str()]->Fill( pullGPY_tr_mono );
00798 title.str("");
00799 title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr_mono";
00800 hPullGP_Z_tr_mono[title.str()]->Fill( pullGPZ_tr_mono );
00801
00802 double pullGPX_ts_mono = (monoTsosGPos.x()-monoShitGPos.x())/sqrt(monoTsosGPEr.cxx());
00803 double pullGPY_ts_mono = (monoTsosGPos.y()-monoShitGPos.y())/sqrt(monoTsosGPEr.cyy());
00804 double pullGPZ_ts_mono = (monoTsosGPos.z()-monoShitGPos.z())/sqrt(monoTsosGPEr.czz());
00805
00806
00807
00808
00809 title.str("");
00810 title << "PullGP_X_" << subdetId << "-" << layerId << "_ts_mono";
00811 hPullGP_X_ts_mono[title.str()]->Fill( pullGPX_ts_mono );
00812 title.str("");
00813 title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts_mono";
00814 hPullGP_Y_ts_mono[title.str()]->Fill( pullGPY_ts_mono );
00815 title.str("");
00816 title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts_mono";
00817 hPullGP_Z_ts_mono[title.str()]->Fill( pullGPZ_ts_mono );
00818
00819 double pullGMX_ts_mono = (monoTsosGMom.x()-monoShitGMom.x())/sqrt(monoTsosGMEr.cxx());
00820 double pullGMY_ts_mono = (monoTsosGMom.y()-monoShitGMom.y())/sqrt(monoTsosGMEr.cyy());
00821 double pullGMZ_ts_mono = (monoTsosGMom.z()-monoShitGMom.z())/sqrt(monoTsosGMEr.czz());
00822
00823
00824
00825
00826 title.str("");
00827 title << "PullGM_X_" << subdetId << "-" << layerId << "_ts_mono";
00828 hPullGM_X_ts_mono[title.str()]->Fill( pullGMX_ts_mono );
00829 title.str("");
00830 title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts_mono";
00831 hPullGM_Y_ts_mono[title.str()]->Fill( pullGMY_ts_mono );
00832 title.str("");
00833 title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts_mono";
00834 hPullGM_Z_ts_mono[title.str()]->Fill( pullGMZ_ts_mono );
00835
00836
00837 LogTrace("TestTrackHits") << "STEREO HIT" ;
00838 auto s = dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())->stereoHit();
00839 CTTRHp tStereoHit = theBuilder->build(&s);
00840 if (tStereoHit==0) continue;
00841 vector<PSimHit> assStereoSimHits = hitAssociator->associateHit(*tStereoHit->hit());
00842 if (assStereoSimHits.size()==0) continue;
00843 const PSimHit sStereoHit = *(assStereoSimHits.begin());
00844 const Surface * stereoSurf = &( tStereoHit->det()->surface() );
00845 if (stereoSurf==0) continue;
00846 TSOS stereoState = thePropagatorAnyDir->propagate(state,*stereoSurf);
00847 if (stereoState.isValid()==0) continue;
00848
00849 LocalVector stereoShitLMom = sStereoHit.momentumAtEntry();
00850 GlobalVector stereoShitGMom = stereoSurf->toGlobal(stereoShitLMom);
00851 LocalPoint stereoShitLPos = sStereoHit.localPosition();
00852 GlobalPoint stereoShitGPos = stereoSurf->toGlobal(stereoShitLPos);
00853
00854
00855
00856
00857 GlobalVector stereoTsosGMom = stereoState.globalMomentum();
00858 GlobalError stereoTsosGMEr(stereoState.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3,3));
00859 GlobalPoint stereoTsosGPos = stereoState.globalPosition();
00860 GlobalError stereoTsosGPEr = stereoState.cartesianError().position();
00861
00862 GlobalPoint stereoRhitGPos = tStereoHit->globalPosition();
00863 GlobalError stereoRhitGPEr = tStereoHit->globalPositionError();
00864
00865 MeasurementExtractor meSt(stereoState);
00866 double chi2stereo = computeChi2Increment(meSt,tStereoHit);
00867
00868 title.str("");
00869 title << "Chi2Increment_stereo_" << subdetId << "-" << layerId ;
00870 hChi2Increment_stereo[title.str()]->Fill(chi2stereo);
00871
00872 double pullGPX_rs_stereo = (stereoRhitGPos.x()-stereoShitGPos.x())/sqrt(stereoRhitGPEr.cxx());
00873 double pullGPY_rs_stereo = (stereoRhitGPos.y()-stereoShitGPos.y())/sqrt(stereoRhitGPEr.cyy());
00874 double pullGPZ_rs_stereo = (stereoRhitGPos.z()-stereoShitGPos.z())/sqrt(stereoRhitGPEr.czz());
00875
00876
00877
00878
00879 title.str("");
00880 title << "PullGP_X_" << subdetId << "-" << layerId << "_rs_stereo";
00881 hPullGP_X_rs_stereo[title.str()]->Fill( pullGPX_rs_stereo );
00882 title.str("");
00883 title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs_stereo";
00884 hPullGP_Y_rs_stereo[title.str()]->Fill( pullGPY_rs_stereo );
00885 title.str("");
00886 title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs_stereo";
00887 hPullGP_Z_rs_stereo[title.str()]->Fill( pullGPZ_rs_stereo );
00888
00889 double pullGPX_tr_stereo = (stereoTsosGPos.x()-stereoRhitGPos.x())/sqrt(stereoTsosGPEr.cxx()+stereoRhitGPEr.cxx());
00890 double pullGPY_tr_stereo = (stereoTsosGPos.y()-stereoRhitGPos.y())/sqrt(stereoTsosGPEr.cyy()+stereoRhitGPEr.cyy());
00891 double pullGPZ_tr_stereo = (stereoTsosGPos.z()-stereoRhitGPos.z())/sqrt(stereoTsosGPEr.czz()+stereoRhitGPEr.czz());
00892
00893
00894
00895
00896 title.str("");
00897 title << "PullGP_X_" << subdetId << "-" << layerId << "_tr_stereo";
00898 hPullGP_X_tr_stereo[title.str()]->Fill( pullGPX_tr_stereo );
00899 title.str("");
00900 title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr_stereo";
00901 hPullGP_Y_tr_stereo[title.str()]->Fill( pullGPY_tr_stereo );
00902 title.str("");
00903 title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr_stereo";
00904 hPullGP_Z_tr_stereo[title.str()]->Fill( pullGPZ_tr_stereo );
00905
00906 double pullGPX_ts_stereo = (stereoTsosGPos.x()-stereoShitGPos.x())/sqrt(stereoTsosGPEr.cxx());
00907 double pullGPY_ts_stereo = (stereoTsosGPos.y()-stereoShitGPos.y())/sqrt(stereoTsosGPEr.cyy());
00908 double pullGPZ_ts_stereo = (stereoTsosGPos.z()-stereoShitGPos.z())/sqrt(stereoTsosGPEr.czz());
00909
00910
00911
00912
00913 title.str("");
00914 title << "PullGP_X_" << subdetId << "-" << layerId << "_ts_stereo";
00915 hPullGP_X_ts_stereo[title.str()]->Fill( pullGPX_ts_stereo );
00916 title.str("");
00917 title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts_stereo";
00918 hPullGP_Y_ts_stereo[title.str()]->Fill( pullGPY_ts_stereo );
00919 title.str("");
00920 title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts_stereo";
00921 hPullGP_Z_ts_stereo[title.str()]->Fill( pullGPZ_ts_stereo );
00922
00923 double pullGMX_ts_stereo = (stereoTsosGMom.x()-stereoShitGMom.x())/sqrt(stereoTsosGMEr.cxx());
00924 double pullGMY_ts_stereo = (stereoTsosGMom.y()-stereoShitGMom.y())/sqrt(stereoTsosGMEr.cyy());
00925 double pullGMZ_ts_stereo = (stereoTsosGMom.z()-stereoShitGMom.z())/sqrt(stereoTsosGMEr.czz());
00926
00927
00928
00929
00930 title.str("");
00931 title << "PullGM_X_" << subdetId << "-" << layerId << "_ts_stereo";
00932 hPullGM_X_ts_stereo[title.str()]->Fill( pullGMX_ts_stereo );
00933 title.str("");
00934 title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts_stereo";
00935 hPullGM_Y_ts_stereo[title.str()]->Fill( pullGMY_ts_stereo );
00936 title.str("");
00937 title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts_stereo";
00938 hPullGM_Z_ts_stereo[title.str()]->Fill( pullGMZ_ts_stereo );
00939 }
00940 }
00941 LogTrace("TestTrackHits") << "traj chi2=" << tchi2 ;
00942 LogTrace("TestTrackHits") << "track chi2=" << track->chi2() ;
00943 i++;
00944 }
00945 LogTrace("TestTrackHits") << "end of event: processd hits=" << evtHits ;
00946 delete hitAssociator;
00947 }
00948
00949 void TestTrackHits::endJob() {
00950
00951 TDirectory * chi2d = file->mkdir("Chi2Increment");
00952
00953 TDirectory * gp_ts = file->mkdir("GP_TSOS-SimHit");
00954 TDirectory * gm_ts = file->mkdir("GM_TSOS-SimHit");
00955 TDirectory * gp_tr = file->mkdir("GP_TSOS-RecHit");
00956 TDirectory * gp_rs = file->mkdir("GP_RecHit-SimHit");
00957
00958 TDirectory * gp_tsx = gp_ts->mkdir("X");
00959 TDirectory * gp_tsy = gp_ts->mkdir("Y");
00960 TDirectory * gp_tsz = gp_ts->mkdir("Z");
00961 TDirectory * gm_tsx = gm_ts->mkdir("X");
00962 TDirectory * gm_tsy = gm_ts->mkdir("Y");
00963 TDirectory * gm_tsz = gm_ts->mkdir("Z");
00964 TDirectory * gp_trx = gp_tr->mkdir("X");
00965 TDirectory * gp_try = gp_tr->mkdir("Y");
00966 TDirectory * gp_trz = gp_tr->mkdir("Z");
00967 TDirectory * gp_rsx = gp_rs->mkdir("X");
00968 TDirectory * gp_rsy = gp_rs->mkdir("Y");
00969 TDirectory * gp_rsz = gp_rs->mkdir("Z");
00970
00971 TDirectory * gp_tsx_mono = gp_ts->mkdir("MONOX");
00972 TDirectory * gp_tsy_mono = gp_ts->mkdir("MONOY");
00973 TDirectory * gp_tsz_mono = gp_ts->mkdir("MONOZ");
00974 TDirectory * gm_tsx_mono = gm_ts->mkdir("MONOX");
00975 TDirectory * gm_tsy_mono = gm_ts->mkdir("MONOY");
00976 TDirectory * gm_tsz_mono = gm_ts->mkdir("MONOZ");
00977 TDirectory * gp_trx_mono = gp_tr->mkdir("MONOX");
00978 TDirectory * gp_try_mono = gp_tr->mkdir("MONOY");
00979 TDirectory * gp_trz_mono = gp_tr->mkdir("MONOZ");
00980 TDirectory * gp_rsx_mono = gp_rs->mkdir("MONOX");
00981 TDirectory * gp_rsy_mono = gp_rs->mkdir("MONOY");
00982 TDirectory * gp_rsz_mono = gp_rs->mkdir("MONOZ");
00983
00984 TDirectory * gp_tsx_stereo = gp_ts->mkdir("STEREOX");
00985 TDirectory * gp_tsy_stereo = gp_ts->mkdir("STEREOY");
00986 TDirectory * gp_tsz_stereo = gp_ts->mkdir("STEREOZ");
00987 TDirectory * gm_tsx_stereo = gm_ts->mkdir("STEREOX");
00988 TDirectory * gm_tsy_stereo = gm_ts->mkdir("STEREOY");
00989 TDirectory * gm_tsz_stereo = gm_ts->mkdir("STEREOZ");
00990 TDirectory * gp_trx_stereo = gp_tr->mkdir("STEREOX");
00991 TDirectory * gp_try_stereo = gp_tr->mkdir("STEREOY");
00992 TDirectory * gp_trz_stereo = gp_tr->mkdir("STEREOZ");
00993 TDirectory * gp_rsx_stereo = gp_rs->mkdir("STEREOX");
00994 TDirectory * gp_rsy_stereo = gp_rs->mkdir("STEREOY");
00995 TDirectory * gp_rsz_stereo = gp_rs->mkdir("STEREOZ");
00996
00997 chi2d->cd();
00998 hTotChi2Increment->Write();
00999 hTotChi2GoodHit->Write();
01000 hTotChi2BadHit->Write();
01001 hTotChi2DeltaHit->Write();
01002 hTotChi2NSharedHit->Write();
01003 hTotChi2SharedHit->Write();
01004 hProcess_vs_Chi2->Write();
01005 hClsize_vs_Chi2->Write();
01006 hPixClsize_vs_Chi2->Write();
01007 hPrjClsize_vs_Chi2->Write();
01008 hSt1Clsize_vs_Chi2->Write();
01009 hSt2Clsize_vs_Chi2->Write();
01010 hGoodHit_vs_Chi2->Write();
01011 hClusterSize->Write();
01012 hPixClusterSize->Write();
01013 hPrjClusterSize->Write();
01014 hSt1ClusterSize->Write();
01015 hSt2ClusterSize->Write();
01016 hSimHitVecSize->Write();
01017 hPixSimHitVecSize->Write();
01018 hPrjSimHitVecSize->Write();
01019 hSt1SimHitVecSize->Write();
01020 hSt2SimHitVecSize->Write();
01021 goodbadmerged->Write();
01022 energyLossRatio->Write();
01023 mergedPull->Write();
01024 probXgood->Write();
01025 probXbad->Write();
01026 probXdelta->Write();
01027 probXshared->Write();
01028 probXnoshare->Write();
01029 probYgood->Write();
01030 probYbad->Write();
01031 probYdelta->Write();
01032 probYshared->Write();
01033 probYnoshare->Write();
01034 for (int i=0; i!=6; i++)
01035 for (int j=0; j!=9; j++){
01036 if (i==0 && j>2) break;
01037 if (i==1 && j>1) break;
01038 if (i==2 && j>3) break;
01039 if (i==3 && j>2) break;
01040 if (i==4 && j>5) break;
01041 if (i==5 && j>8) break;
01042 chi2d->cd();
01043 title.str("");
01044 title << "Chi2Increment_" << i+1 << "-" << j+1 ;
01045 hChi2Increment[title.str()]->Write();
01046 title.str("");
01047 title << "Chi2IncrementVsEta_" << i+1 << "-" << j+1 ;
01048 hChi2IncrementVsEta[title.str()]->Write();
01049 title.str("");
01050 title << "Chi2GoodHit_" << i+1 << "-" << j+1 ;
01051 hChi2GoodHit[title.str()]->Write();
01052 title.str("");
01053 title << "Chi2BadHit_" << i+1 << "-" << j+1 ;
01054 hChi2BadHit[title.str()]->Write();
01055 title.str("");
01056 title << "Chi2DeltaHit_" << i+1 << "-" << j+1 ;
01057 hChi2DeltaHit[title.str()]->Write();
01058 title.str("");
01059 title << "Chi2NSharedHit_" << i+1 << "-" << j+1 ;
01060 hChi2NSharedHit[title.str()]->Write();
01061 title.str("");
01062 title << "Chi2SharedHit_" << i+1 << "-" << j+1 ;
01063 hChi2SharedHit[title.str()]->Write();
01064
01065 gp_ts->cd();
01066 gp_tsx->cd();
01067 title.str("");
01068 title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts";
01069 hPullGP_X_ts[title.str()]->Write();
01070 gp_tsy->cd();
01071 title.str("");
01072 title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts";
01073 hPullGP_Y_ts[title.str()]->Write();
01074 gp_tsz->cd();
01075 title.str("");
01076 title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts";
01077 hPullGP_Z_ts[title.str()]->Write();
01078
01079 gm_ts->cd();
01080 gm_tsx->cd();
01081 title.str("");
01082 title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts";
01083 hPullGM_X_ts[title.str()]->Write();
01084 gm_tsy->cd();
01085 title.str("");
01086 title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts";
01087 hPullGM_Y_ts[title.str()]->Write();
01088 gm_tsz->cd();
01089 title.str("");
01090 title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts";
01091 hPullGM_Z_ts[title.str()]->Write();
01092
01093 gp_tr->cd();
01094 gp_trx->cd();
01095 title.str("");
01096 title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr";
01097 hPullGP_X_tr[title.str()]->Write();
01098 gp_try->cd();
01099 title.str("");
01100 title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr";
01101 hPullGP_Y_tr[title.str()]->Write();
01102 gp_trz->cd();
01103 title.str("");
01104 title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr";
01105 hPullGP_Z_tr[title.str()]->Write();
01106
01107 gp_rs->cd();
01108 gp_rsx->cd();
01109 title.str("");
01110 title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs";
01111 hPullGP_X_rs[title.str()]->Write();
01112 gp_rsy->cd();
01113 title.str("");
01114 title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs";
01115 hPullGP_Y_rs[title.str()]->Write();
01116 gp_rsz->cd();
01117 title.str("");
01118 title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs";
01119 hPullGP_Z_rs[title.str()]->Write();
01120
01121 if ( ((i==2||i==4)&&(j==0||j==1)) || (i==3||i==5) ){
01122 chi2d->cd();
01123 title.str("");
01124 title << "Chi2Increment_mono_" << i+1 << "-" << j+1 ;
01125 hChi2Increment_mono[title.str()]->Write();
01126 title.str("");
01127 title << "Chi2Increment_stereo_" << i+1 << "-" << j+1 ;
01128 hChi2Increment_stereo[title.str()]->Write();
01129
01130 gp_ts->cd();
01131 gp_tsx_mono->cd();
01132 title.str("");
01133 title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_mono";
01134 hPullGP_X_ts_mono[title.str()]->Write();
01135 gp_tsy_mono->cd();
01136 title.str("");
01137 title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts_mono";
01138 hPullGP_Y_ts_mono[title.str()]->Write();
01139 gp_tsz_mono->cd();
01140 title.str("");
01141 title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts_mono";
01142 hPullGP_Z_ts_mono[title.str()]->Write();
01143
01144 gm_ts->cd();
01145 gm_tsx_mono->cd();
01146 title.str("");
01147 title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_mono";
01148 hPullGM_X_ts_mono[title.str()]->Write();
01149 gm_tsy_mono->cd();
01150 title.str("");
01151 title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_mono";
01152 hPullGM_Y_ts_mono[title.str()]->Write();
01153 gm_tsz_mono->cd();
01154 title.str("");
01155 title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_mono";
01156 hPullGM_Z_ts_mono[title.str()]->Write();
01157
01158 gp_tr->cd();
01159 gp_trx_mono->cd();
01160 title.str("");
01161 title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_mono";
01162 hPullGP_X_tr_mono[title.str()]->Write();
01163 gp_try_mono->cd();
01164 title.str("");
01165 title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_mono";
01166 hPullGP_Y_tr_mono[title.str()]->Write();
01167 gp_trz_mono->cd();
01168 title.str("");
01169 title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_mono";
01170 hPullGP_Z_tr_mono[title.str()]->Write();
01171
01172 gp_rs->cd();
01173 gp_rsx_mono->cd();
01174 title.str("");
01175 title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_mono";
01176 hPullGP_X_rs_mono[title.str()]->Write();
01177 gp_rsy_mono->cd();
01178 title.str("");
01179 title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs_mono";
01180 hPullGP_Y_rs_mono[title.str()]->Write();
01181 gp_rsz_mono->cd();
01182 title.str("");
01183 title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs_mono";
01184 hPullGP_Z_rs_mono[title.str()]->Write();
01185
01186
01187 gp_ts->cd();
01188 gp_tsx_stereo->cd();
01189 title.str("");
01190 title << "PullGP_X_" << i+1 << "-" << j+1 << "_ts_stereo";
01191 hPullGP_X_ts_stereo[title.str()]->Write();
01192 gp_tsy_stereo->cd();
01193 title.str("");
01194 title << "PullGP_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
01195 hPullGP_Y_ts_stereo[title.str()]->Write();
01196 gp_tsz_stereo->cd();
01197 title.str("");
01198 title << "PullGP_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
01199 hPullGP_Z_ts_stereo[title.str()]->Write();
01200
01201 gm_ts->cd();
01202 gm_tsx_stereo->cd();
01203 title.str("");
01204 title << "PullGM_X_" << i+1 << "-" << j+1 << "_ts_stereo";
01205 hPullGM_X_ts_stereo[title.str()]->Write();
01206 gm_tsy_stereo->cd();
01207 title.str("");
01208 title << "PullGM_Y_" << i+1 << "-" << j+1 << "_ts_stereo";
01209 hPullGM_Y_ts_stereo[title.str()]->Write();
01210 gm_tsz_stereo->cd();
01211 title.str("");
01212 title << "PullGM_Z_" << i+1 << "-" << j+1 << "_ts_stereo";
01213 hPullGM_Z_ts_stereo[title.str()]->Write();
01214
01215 gp_tr->cd();
01216 gp_trx_stereo->cd();
01217 title.str("");
01218 title << "PullGP_X_" << i+1 << "-" << j+1 << "_tr_stereo";
01219 hPullGP_X_tr_stereo[title.str()]->Write();
01220 gp_try_stereo->cd();
01221 title.str("");
01222 title << "PullGP_Y_" << i+1 << "-" << j+1 << "_tr_stereo";
01223 hPullGP_Y_tr_stereo[title.str()]->Write();
01224 gp_trz_stereo->cd();
01225 title.str("");
01226 title << "PullGP_Z_" << i+1 << "-" << j+1 << "_tr_stereo";
01227 hPullGP_Z_tr_stereo[title.str()]->Write();
01228
01229 gp_rs->cd();
01230 gp_rsx_stereo->cd();
01231 title.str("");
01232 title << "PullGP_X_" << i+1 << "-" << j+1 << "_rs_stereo";
01233 hPullGP_X_rs_stereo[title.str()]->Write();
01234 gp_rsy_stereo->cd();
01235 title.str("");
01236 title << "PullGP_Y_" << i+1 << "-" << j+1 << "_rs_stereo";
01237 hPullGP_Y_rs_stereo[title.str()]->Write();
01238 gp_rsz_stereo->cd();
01239 title.str("");
01240 title << "PullGP_Z_" << i+1 << "-" << j+1 << "_rs_stereo";
01241 hPullGP_Z_rs_stereo[title.str()]->Write();
01242 }
01243 }
01244
01245 file->Close();
01246 }
01247
01248
01249
01250
01251 std::pair<LocalPoint,LocalVector>
01252 TestTrackHits::projectHit( const PSimHit& hit, const StripGeomDetUnit* stripDet, const BoundPlane& plane)
01253 {
01254 const StripTopology& topol = stripDet->specificTopology();
01255 GlobalPoint globalpos= stripDet->surface().toGlobal(hit.localPosition());
01256 LocalPoint localHit = plane.toLocal(globalpos);
01257
01258 LocalVector locdir=hit.localDirection();
01259
01260
01261 GlobalVector globaldir= stripDet->surface().toGlobal(locdir);
01262 LocalVector dir=plane.toLocal(globaldir);
01263 float scale = -localHit.z() / dir.z();
01264
01265 LocalPoint projectedPos = localHit + scale*dir;
01266
01267 float selfAngle = topol.stripAngle( topol.strip( hit.localPosition()));
01268
01269 LocalVector stripDir( sin(selfAngle), cos(selfAngle), 0);
01270
01271 LocalVector localStripDir( plane.toLocal(stripDet->surface().toGlobal( stripDir)));
01272
01273 return std::pair<LocalPoint,LocalVector>( projectedPos, localStripDir);
01274 }
01275
01276 template<unsigned int D>
01277 double TestTrackHits::computeChi2Increment(MeasurementExtractor me,
01278 TransientTrackingRecHit::ConstRecHitPointer rhit) {
01279 typedef typename AlgebraicROOTObject<D>::Vector VecD;
01280 typedef typename AlgebraicROOTObject<D,D>::SymMatrix SMatDD;
01281 VecD r = asSVector<D>(rhit->parameters()) - me.measuredParameters<D>(*rhit);
01282
01283 SMatDD R = asSMatrix<D>(rhit->parametersError()) + me.measuredError<D>(*rhit);
01284 R.Invert();
01285 return ROOT::Math::Similarity(r,R) ;
01286 }
01287
01288 #include "FWCore/Framework/interface/ModuleFactory.h"
01289 #include "FWCore/Framework/interface/MakerMacros.h"
01290 DEFINE_FWK_MODULE(TestTrackHits);
1.7.3