00001 #include "Validation/RecoTrack/interface/MultiTrackValidator.h"
00002 #include "Validation/Tools/interface/FitSlicesYTool.h"
00003
00004 #include "FWCore/Framework/interface/MakerMacros.h"
00005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00006
00007 #include "DataFormats/TrackReco/interface/Track.h"
00008 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00009 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
00010 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
00011 #include "SimTracker/TrackAssociation/interface/TrackAssociatorByChi2.h"
00012 #include "SimTracker/TrackAssociation/interface/TrackAssociatorByHits.h"
00013 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
00014 #include "SimTracker/Records/interface/TrackAssociatorRecord.h"
00015 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
00016 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
00017 #include "TrackingTools/PatternTools/interface/TrajectoryStateClosestToBeamLineBuilder.h"
00018
00019
00020 #include "TMath.h"
00021 #include <TF1.h>
00022
00023 using namespace std;
00024 using namespace edm;
00025
00026 void MultiTrackValidator::beginRun(Run const&, EventSetup const& setup) {
00027
00028
00029
00030 int j=0;
00031 for (unsigned int ww=0;ww<associators.size();ww++){
00032 for (unsigned int www=0;www<label.size();www++){
00033
00034 dbe_->cd();
00035 InputTag algo = label[www];
00036 string dirName=dirName_;
00037 if (algo.process()!="")
00038 dirName+=algo.process()+"_";
00039 if(algo.label()!="")
00040 dirName+=algo.label()+"_";
00041 if(algo.instance()!="")
00042 dirName+=algo.instance()+"_";
00043 if (dirName.find("Tracks")<dirName.length()){
00044 dirName.replace(dirName.find("Tracks"),6,"");
00045 }
00046 string assoc= associators[ww];
00047 if (assoc.find("Track")<assoc.length()){
00048 assoc.replace(assoc.find("Track"),5,"");
00049 }
00050 dirName+=assoc;
00051 std::replace(dirName.begin(), dirName.end(), ':', '_');
00052 dbe_->setCurrentFolder(dirName.c_str());
00053
00054 setUpVectors();
00055
00056 dbe_->goUp();
00057 string subDirName = dirName + "/simulation";
00058 dbe_->setCurrentFolder(subDirName.c_str());
00059 h_ptSIM.push_back( dbe_->book1D("ptSIM", "generated p_{t}", 5500, 0, 110 ) );
00060 h_etaSIM.push_back( dbe_->book1D("etaSIM", "generated pseudorapidity", 500, -2.5, 2.5 ) );
00061 h_tracksSIM.push_back( dbe_->book1D("tracksSIM","number of simluated tracks",100,-0.5,99.5) );
00062 h_vertposSIM.push_back( dbe_->book1D("vertposSIM","Transverse position of sim vertices",100,0.,120.) );
00063
00064 dbe_->cd();
00065 dbe_->setCurrentFolder(dirName.c_str());
00066 h_tracks.push_back( dbe_->book1D("tracks","number of reconstructed tracks",20,-0.5,19.5) );
00067 h_fakes.push_back( dbe_->book1D("fakes","number of fake reco tracks",20,-0.5,19.5) );
00068 h_charge.push_back( dbe_->book1D("charge","charge",3,-1.5,1.5) );
00069 h_hits.push_back( dbe_->book1D("hits", "number of hits per track", nintHit,minHit,maxHit ) );
00070 h_losthits.push_back( dbe_->book1D("losthits", "number of lost hits per track", nintHit,minHit,maxHit) );
00071 h_nchi2.push_back( dbe_->book1D("chi2", "normalized #chi^{2}", 200, 0, 20 ) );
00072 h_nchi2_prob.push_back( dbe_->book1D("chi2_prob", "normalized #chi^{2} probability",100,0,1));
00073
00074 h_effic.push_back( dbe_->book1D("effic","efficiency vs #eta",nint,min,max) );
00075 h_efficPt.push_back( dbe_->book1D("efficPt","efficiency vs pT",nintpT,minpT,maxpT) );
00076 h_fakerate.push_back( dbe_->book1D("fakerate","fake rate vs #eta",nint,min,max) );
00077 h_fakeratePt.push_back( dbe_->book1D("fakeratePt","fake rate vs pT",nintpT,minpT,maxpT) );
00078 h_effic_vs_hit.push_back( dbe_->book1D("effic_vs_hit","effic vs hit",nintHit,minHit,maxHit) );
00079 h_fake_vs_hit.push_back( dbe_->book1D("fakerate_vs_hit","fake rate vs hit",nintHit,minHit,maxHit) );
00080
00081 h_recoeta.push_back( dbe_->book1D("num_reco_eta","N of reco track vs eta",nint,min,max) );
00082 h_assoceta.push_back( dbe_->book1D("num_assoc(simToReco)_eta","N of associated tracks (simToReco) vs eta",nint,min,max) );
00083 h_assoc2eta.push_back( dbe_->book1D("num_assoc(recoToSim)_eta","N of associated (recoToSim) tracks vs eta",nint,min,max) );
00084 h_simuleta.push_back( dbe_->book1D("num_simul_eta","N of simulated tracks vs eta",nint,min,max) );
00085 h_recopT.push_back( dbe_->book1D("num_reco_pT","N of reco track vs pT",nintpT,minpT,maxpT) );
00086 h_assocpT.push_back( dbe_->book1D("num_assoc(simToReco)_pT","N of associated tracks (simToReco) vs pT",nintpT,minpT,maxpT) );
00087 h_assoc2pT.push_back( dbe_->book1D("num_assoc(recoToSim)_pT","N of associated (recoToSim) tracks vs pT",nintpT,minpT,maxpT) );
00088 h_simulpT.push_back( dbe_->book1D("num_simul_pT","N of simulated tracks vs pT",nintpT,minpT,maxpT) );
00089 h_recohit.push_back( dbe_->book1D("num_reco_hit","N of reco track vs hit",nintHit,minHit,maxHit) );
00090 h_assochit.push_back( dbe_->book1D("num_assoc(simToReco)_hit","N of associated tracks (simToReco) vs hit",nintHit,minHit,maxHit) );
00091 h_assoc2hit.push_back( dbe_->book1D("num_assoc(recoToSim)_hit","N of associated (recoToSim) tracks vs hit",nintHit,minHit,maxHit) );
00092 h_simulhit.push_back( dbe_->book1D("num_simul_hit","N of simulated tracks vs hit",nintHit,minHit,maxHit) );
00093
00094 h_eta.push_back( dbe_->book1D("eta", "pseudorapidity residue", 1000, -0.1, 0.1 ) );
00095 h_pt.push_back( dbe_->book1D("pullPt", "pull of p_{t}", 100, -10, 10 ) );
00096 h_pullTheta.push_back( dbe_->book1D("pullTheta","pull of #theta parameter",250,-25,25) );
00097 h_pullPhi.push_back( dbe_->book1D("pullPhi","pull of #phi parameter",250,-25,25) );
00098 h_pullDxy.push_back( dbe_->book1D("pullDxy","pull of dxy parameter",250,-25,25) );
00099 h_pullDz.push_back( dbe_->book1D("pullDz","pull of dz parameter",250,-25,25) );
00100 h_pullQoverp.push_back( dbe_->book1D("pullQoverp","pull of qoverp parameter",250,-25,25) );
00101
00102 if (associators[ww]=="TrackAssociatorByChi2"){
00103 h_assochi2.push_back( dbe_->book1D("assocChi2","track association #chi^{2}",1000000,0,100000) );
00104 h_assochi2_prob.push_back(dbe_->book1D("assocChi2_prob","probability of association #chi^{2}",100,0,1));
00105 } else if (associators[ww]=="TrackAssociatorByHits"){
00106 h_assocFraction.push_back( dbe_->book1D("assocFraction","fraction of shared hits",200,0,2) );
00107 h_assocSharedHit.push_back(dbe_->book1D("assocSharedHit","number of shared hits",20,0,20));
00108 }
00109
00110 chi2_vs_nhits.push_back( dbe_->book2D("chi2_vs_nhits","#chi^{2} vs nhits",25,0,25,100,0,10) );
00111 h_chi2meanhitsh.push_back( dbe_->book1D("chi2mean_vs_nhits","mean #chi^{2} vs nhits",25,0,25) );
00112
00113 etares_vs_eta.push_back( dbe_->book2D("etares_vs_eta","etaresidue vs eta",nint,min,max,200,-0.1,0.1) );
00114 nrec_vs_nsim.push_back( dbe_->book2D("nrec_vs_nsim","nrec vs nsim",20,-0.5,19.5,20,-0.5,19.5) );
00115
00116 chi2_vs_eta.push_back( dbe_->book2D("chi2_vs_eta","chi2_vs_eta",nint,min,max, 200, 0, 20 ));
00117 h_chi2meanh.push_back( dbe_->book1D("chi2mean","mean #chi^{2} vs #eta",nint,min,max) );
00118 chi2_vs_phi.push_back( dbe_->book2D("chi2_vs_phi","#chi^{2} vs #phi",nintPhi,minPhi,maxPhi, 200, 0, 20 ) );
00119 h_chi2mean_vs_phi.push_back( dbe_->book1D("chi2mean_vs_phi","mean of #chi^{2} vs #phi",nintPhi,minPhi,maxPhi) );
00120
00121 nhits_vs_eta.push_back( dbe_->book2D("nhits_vs_eta","nhits vs eta",nint,min,max,25,0,25) );
00122 h_hits_eta.push_back( dbe_->book1D("hits_eta","mean #hits vs eta",nint,min,max) );
00123 nhits_vs_phi.push_back( dbe_->book2D("nhits_vs_phi","#hits vs #phi",nintPhi,minPhi,maxPhi,25,0,25) );
00124 h_hits_phi.push_back( dbe_->book1D("hits_phi","mean #hits vs #phi",nintPhi,minPhi,maxPhi) );
00125
00126 nlosthits_vs_eta.push_back( dbe_->book2D("nlosthits_vs_eta","nlosthits vs eta",nint,min,max,25,0,25) );
00127 h_losthits_eta.push_back( dbe_->book1D("losthits_eta","losthits_eta",nint,min,max) );
00128
00129
00130
00131
00132
00133
00134
00135 ptres_vs_eta.push_back(dbe_->book2D("ptres_vs_eta","ptres_vs_eta",nint,min,max, 100, -0.1, 0.1));
00136 h_ptresmean_vs_eta.push_back( dbe_->book1D("ptresmean_vs_eta","mean of p_{t} resolution vs #eta",nint,min,max));
00137 h_ptrmsh.push_back( dbe_->book1D("sigmapt","#sigma(#deltap_{t}/p_{t}) vs #eta",nint,min,max) );
00138
00139 ptres_vs_phi.push_back( dbe_->book2D("ptres_vs_phi","p_{t} res vs #phi",nintPhi,minPhi,maxPhi, 100, -0.1, 0.1));
00140 h_ptresmean_vs_phi.push_back( dbe_->book1D("ptresmean_vs_phi","mean of p_{t} resolution vs #phi",nintPhi,minPhi,maxPhi));
00141 h_ptrmshPhi.push_back( dbe_->book1D("sigmaptPhi","#sigma(#deltap_{t}/p_{t}) vs #phi",nintPhi,minPhi,maxPhi) );
00142
00143 ptres_vs_pt.push_back(dbe_->book2D("ptres_vs_pt","ptres_vs_pt",nintpT,minpT,maxpT, 100, -0.1, 0.1));
00144 h_ptrmshPt.push_back( dbe_->book1D("sigmaptPt","#sigma(#deltap_{t}/p_{t}) vs pT",nintpT,minpT,maxpT) );
00145
00146 cotThetares_vs_eta.push_back(dbe_->book2D("cotThetares_vs_eta","cotThetares_vs_eta",nint,min,max, 120, -0.01, 0.01));
00147 h_cotThetarmsh.push_back( dbe_->book1D("sigmacotTheta","#sigma(#deltacot(#theta)) vs #eta",nint,min,max) );
00148
00149 cotThetares_vs_pt.push_back(dbe_->book2D("cotThetares_vs_pt","cotThetares_vs_pt",nintpT,minpT,maxpT, 120, -0.01, 0.01));
00150 h_cotThetarmshPt.push_back( dbe_->book1D("sigmacotThetaPt","#sigma(#deltacot(#theta)) vs pT",nintpT,minpT,maxpT) );
00151
00152 phires_vs_eta.push_back(dbe_->book2D("phires_vs_eta","phires_vs_eta",nint,min,max, 100, -0.003, 0.003));
00153 h_phiresmean_vs_eta.push_back(dbe_->book1D("phiresmean_vs_eta","mean of #phi res vs #eta",nint,min,max));
00154 h_phirmsh.push_back( dbe_->book1D("sigmaphi","#sigma(#delta#phi) vs #eta",nint,min,max) );
00155
00156 phires_vs_pt.push_back(dbe_->book2D("phires_vs_pt","phires_vs_pt",nintpT,minpT,maxpT, 100, -0.003, 0.003));
00157 h_phirmshPt.push_back( dbe_->book1D("sigmaphiPt","#sigma(#delta#phi) vs pT",nintpT,minpT,maxpT) );
00158
00159 phires_vs_phi.push_back(dbe_->book2D("phires_vs_phi","#phi res vs #phi",nintPhi,minPhi,maxPhi, 100, -0.003, 0.003));
00160 h_phiresmean_vs_phi.push_back(dbe_->book1D("phiresmean_vs_phi","mean of #phi res vs #phi",nintPhi,minPhi,maxPhi));
00161 h_phirmshPhi.push_back( dbe_->book1D("sigmaphiPhi","#sigma(#delta#phi) vs #phi",nintPhi,minPhi,maxPhi) );
00162
00163 dxyres_vs_eta.push_back(dbe_->book2D("dxyres_vs_eta","dxyres_vs_eta",nint,min,max, 100, -0.01, 0.01));
00164 h_dxyrmsh.push_back( dbe_->book1D("sigmadxy","#sigma(#deltadxy) vs #eta",nint,min,max) );
00165
00166 dxyres_vs_pt.push_back( dbe_->book2D("dxyres_vs_pt","dxyres_vs_pt",nintpT,minpT,maxpT, 100, -0.01, 0.01));
00167 h_dxyrmshPt.push_back( dbe_->book1D("sigmadxyPt","#sigmadxy vs pT",nintpT,minpT,maxpT) );
00168
00169 dzres_vs_eta.push_back(dbe_->book2D("dzres_vs_eta","dzres_vs_eta",nint,min,max, 150, -0.05, 0.05));
00170 h_dzrmsh.push_back( dbe_->book1D("sigmadz","#sigma(#deltadz) vs #eta",nint,min,max) );
00171
00172 dzres_vs_pt.push_back(dbe_->book2D("dzres_vs_pt","dzres_vs_pt",nintpT,minpT,maxpT, 150, -0.05, 0.05));
00173 h_dzrmshPt.push_back( dbe_->book1D("sigmadzPt","#sigma(#deltadz vs pT",nintpT,minpT,maxpT) );
00174
00175 ptmean_vs_eta_phi.push_back(dbe_->bookProfile2D("ptmean_vs_eta_phi","mean p_{t} vs #eta and #phi",nintPhi,minPhi,maxPhi,nint,min,max,nintpT,minpT,maxpT));
00176 phimean_vs_eta_phi.push_back(dbe_->bookProfile2D("phimean_vs_eta_phi","mean #phi vs #eta and #phi",nintPhi,minPhi,maxPhi,nint,min,max,nintPhi,minPhi,maxPhi));
00177
00178
00179 dxypull_vs_eta.push_back(dbe_->book2D("dxypull_vs_eta","dxypull_vs_eta",nint,min,max,100,-10,10));
00180 ptpull_vs_eta.push_back(dbe_->book2D("ptpull_vs_eta","ptpull_vs_eta",nint,min,max,100,-10,10));
00181 dzpull_vs_eta.push_back(dbe_->book2D("dzpull_vs_eta","dzpull_vs_eta",nint,min,max,100,-10,10));
00182 phipull_vs_eta.push_back(dbe_->book2D("phipull_vs_eta","phipull_vs_eta",nint,min,max,100,-10,10));
00183 thetapull_vs_eta.push_back(dbe_->book2D("thetapull_vs_eta","thetapull_vs_eta",nint,min,max,100,-10,10));
00184 h_dxypulleta.push_back( dbe_->book1D("h_dxypulleta","#sigma of dxy pull vs #eta",nint,min,max) );
00185 h_ptpulleta.push_back( dbe_->book1D("h_ptpulleta","#sigma of p_{t} pull vs #eta",nint,min,max) );
00186 h_dzpulleta.push_back( dbe_->book1D("h_dzpulleta","#sigma of dz pull vs #eta",nint,min,max) );
00187 h_phipulleta.push_back( dbe_->book1D("h_phipulleta","#sigma of #phi pull vs #eta",nint,min,max) );
00188 h_thetapulleta.push_back( dbe_->book1D("h_thetapulleta","#sigma of #theta pull vs #eta",nint,min,max) );
00189 h_ptshifteta.push_back( dbe_->book1D("h_ptshifteta","<#deltapT/pT>[%] vs #eta",nint,min,max) );
00190
00191
00192 ptpull_vs_phi.push_back(dbe_->book2D("ptpull_vs_phi","p_{t} pull vs #phi",nintPhi,minPhi,maxPhi,100,-10,10));
00193 phipull_vs_phi.push_back(dbe_->book2D("phipull_vs_phi","#phi pull vs #phi",nintPhi,minPhi,maxPhi,100,-10,10));
00194 thetapull_vs_phi.push_back(dbe_->book2D("thetapull_vs_phi","#theta pull vs #phi",nintPhi,minPhi,maxPhi,100,-10,10));
00195 h_ptpullphi.push_back( dbe_->book1D("h_ptpullphi","#sigma of p_{t} pull vs #phi",nintPhi,minPhi,maxPhi) );
00196 h_phipullphi.push_back( dbe_->book1D("h_phipullphi","#sigma of #phi pull vs #phi",nintPhi,minPhi,maxPhi) );
00197 h_thetapullphi.push_back( dbe_->book1D("h_thetapullphi","#sigma of #theta pull vs #phi",nintPhi,minPhi,maxPhi) );
00198
00199 j++;
00200 }
00201 }
00202 if (UseAssociators) {
00203 edm::ESHandle<TrackAssociatorBase> theAssociator;
00204 for (unsigned int w=0;w<associators.size();w++) {
00205 setup.get<TrackAssociatorRecord>().get(associators[w],theAssociator);
00206 associator.push_back( theAssociator.product() );
00207 }
00208 }
00209 }
00210
00211 void MultiTrackValidator::analyze(const edm::Event& event, const edm::EventSetup& setup){
00212 using namespace reco;
00213
00214 edm::LogInfo("TrackValidator") << "\n====================================================" << "\n"
00215 << "Analyzing new event" << "\n"
00216 << "====================================================\n" << "\n";
00217
00218 edm::Handle<TrackingParticleCollection> TPCollectionHeff ;
00219 event.getByLabel(label_tp_effic,TPCollectionHeff);
00220 const TrackingParticleCollection tPCeff = *(TPCollectionHeff.product());
00221
00222 edm::Handle<TrackingParticleCollection> TPCollectionHfake ;
00223 event.getByLabel(label_tp_fake,TPCollectionHfake);
00224 const TrackingParticleCollection tPCfake = *(TPCollectionHfake.product());
00225
00226
00227
00228
00229 edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
00230 event.getByLabel(bsSrc,recoBeamSpotHandle);
00231 reco::BeamSpot bs = *recoBeamSpotHandle;
00232
00233 int w=0;
00234 for (unsigned int ww=0;ww<associators.size();ww++){
00235 for (unsigned int www=0;www<label.size();www++){
00236
00237
00238
00239 edm::Handle<View<Track> > trackCollection;
00240 event.getByLabel(label[www], trackCollection);
00241
00242
00243
00244
00245 reco::RecoToSimCollection recSimColl;
00246 reco::SimToRecoCollection simRecColl;
00247
00248 if(UseAssociators){
00249 edm::LogVerbatim("TrackValidator") << "Analyzing "
00250 << label[www].process()<<":"
00251 << label[www].label()<<":"
00252 << label[www].instance()<<" with "
00253 << associators[ww].c_str() <<"\n";
00254
00255 LogTrace("TrackValidator") << "Calling associateRecoToSim method" << "\n";
00256 recSimColl=associator[ww]->associateRecoToSim(trackCollection,
00257 TPCollectionHfake,
00258 &event);
00259 LogTrace("TrackValidator") << "Calling associateSimToReco method" << "\n";
00260 simRecColl=associator[ww]->associateSimToReco(trackCollection,
00261 TPCollectionHeff,
00262 &event);
00263 }
00264 else{
00265 edm::LogVerbatim("TrackValidator") << "Analyzing "
00266 << label[www].process()<<":"
00267 << label[www].label()<<":"
00268 << label[www].instance()<<" with "
00269 << associatormap.process()<<":"
00270 << associatormap.label()<<":"
00271 << associatormap.instance()<<"\n";
00272
00273 Handle<reco::SimToRecoCollection > simtorecoCollectionH;
00274 event.getByLabel(associatormap,simtorecoCollectionH);
00275 simRecColl= *(simtorecoCollectionH.product());
00276
00277 Handle<reco::RecoToSimCollection > recotosimCollectionH;
00278 event.getByLabel(associatormap,recotosimCollectionH);
00279 recSimColl= *(recotosimCollectionH.product());
00280 }
00281
00282
00283
00284
00285 edm::LogVerbatim("TrackValidator") << "\n# of TrackingParticles: " << tPCeff.size() << "\n";
00286 int ats = 0;
00287 int st=0;
00288 for (TrackingParticleCollection::size_type i=0; i<tPCeff.size(); i++){
00289 TrackingParticleRef tpr(TPCollectionHeff, i);
00290 TrackingParticle* tp=const_cast<TrackingParticle*>(tpr.get());
00291 if( (! tpSelector(*tp))) continue;
00292 st++;
00293 h_ptSIM[w]->Fill(sqrt(tp->momentum().perp2()));
00294 h_etaSIM[w]->Fill(tp->momentum().eta());
00295 h_vertposSIM[w]->Fill(sqrt(tp->vertex().perp2()));
00296
00297 std::vector<std::pair<RefToBase<Track>, double> > rt;
00298 if(simRecColl.find(tpr) != simRecColl.end()){
00299 rt = (std::vector<std::pair<RefToBase<Track>, double> >) simRecColl[tpr];
00300 if (rt.size()!=0) {
00301 ats++;
00302 edm::LogVerbatim("TrackValidator") << "TrackingParticle #" << st
00303 << " with pt=" << sqrt(tp->momentum().perp2())
00304 << " associated with quality:" << rt.begin()->second <<"\n";
00305 }
00306 }else{
00307 edm::LogVerbatim("TrackValidator") << "TrackingParticle #" << st
00308 << " with pt=" << sqrt(tp->momentum().perp2())
00309 << " NOT associated to any reco::Track" << "\n";
00310 }
00311
00312 for (unsigned int f=0; f<etaintervals[w].size()-1; f++){
00313 if (getEta(tp->momentum().eta())>etaintervals[w][f]&&
00314 getEta(tp->momentum().eta())<etaintervals[w][f+1]) {
00315 totSIMeta[w][f]++;
00316 if (rt.size()!=0) {
00317 totASSeta[w][f]++;
00318 }
00319 }
00320 }
00321
00322 for (unsigned int f=0; f<pTintervals[w].size()-1; f++){
00323 if (getPt(sqrt(tp->momentum().perp2()))>pTintervals[w][f]&&
00324 getPt(sqrt(tp->momentum().perp2()))<pTintervals[w][f+1]) {
00325 totSIMpT[w][f]++;
00326 if (rt.size()!=0) {
00327 totASSpT[w][f]++;
00328 }
00329 }
00330 }
00331 int tmp = std::min((int)(tp->trackerPSimHit_end()-tp->trackerPSimHit_begin()),int(maxHit-1));
00332 totSIM_hit[w][tmp]++;
00333 if (rt.size()!=0) totASS_hit[w][tmp]++;
00334 }
00335 if (st!=0) h_tracksSIM[w]->Fill(st);
00336
00337
00338
00339
00340
00341 edm::LogVerbatim("TrackValidator") << "\n# of reco::Tracks with "
00342 << label[www].process()<<":"
00343 << label[www].label()<<":"
00344 << label[www].instance()
00345 << ": " << trackCollection->size() << "\n";
00346 int at=0;
00347 int rT=0;
00348 for(View<Track>::size_type i=0; i<trackCollection->size(); ++i){
00349 RefToBase<Track> track(trackCollection, i);
00350 rT++;
00351
00352 std::vector<std::pair<TrackingParticleRef, double> > tp;
00353 if(recSimColl.find(track) != recSimColl.end()){
00354 tp = recSimColl[track];
00355 if (tp.size()!=0) {
00356 at++;
00357 edm::LogVerbatim("TrackValidator") << "reco::Track #" << rT << " with pt=" << track->pt()
00358 << " associated with quality:" << tp.begin()->second <<"\n";
00359 }
00360 } else {
00361 edm::LogVerbatim("TrackValidator") << "reco::Track #" << rT << " with pt=" << track->pt()
00362 << " NOT associated to any TrackingParticle" << "\n";
00363 }
00364
00365
00366 for (unsigned int f=0; f<etaintervals[w].size()-1; f++){
00367 if (getEta(track->momentum().eta())>etaintervals[w][f]&&
00368 getEta(track->momentum().eta())<etaintervals[w][f+1]) {
00369 totRECeta[w][f]++;
00370 if (tp.size()!=0) {
00371 totASS2eta[w][f]++;
00372 }
00373 }
00374 }
00375
00376 for (unsigned int f=0; f<pTintervals[w].size()-1; f++){
00377 if (getPt(sqrt(track->momentum().perp2()))>pTintervals[w][f]&&
00378 getPt(sqrt(track->momentum().perp2()))<pTintervals[w][f+1]) {
00379 totRECpT[w][f]++;
00380 if (tp.size()!=0) {
00381 totASS2pT[w][f]++;
00382 }
00383 }
00384 }
00385 int tmp = std::min((int)track->found(),int(maxHit-1));
00386 totREC_hit[w][tmp]++;
00387 if (tp.size()!=0) totASS2_hit[w][tmp]++;
00388
00389
00390 try{
00391 if (tp.size()==0) continue;
00392
00393 TrackingParticleRef tpr = tp.begin()->first;
00394 const SimTrack * assocTrack = &(*tpr->g4Track_begin());
00395
00396 if (associators[ww]=="TrackAssociatorByChi2"){
00397
00398 double assocChi2 = -tp.begin()->second;
00399 h_assochi2[www]->Fill(assocChi2);
00400 h_assochi2_prob[www]->Fill(TMath::Prob((assocChi2)*5,5));
00401 }
00402 else if (associators[ww]=="TrackAssociatorByHits"){
00403 double fraction = tp.begin()->second;
00404 h_assocFraction[www]->Fill(fraction);
00405 h_assocSharedHit[www]->Fill(fraction*track->numberOfValidHits());
00406 }
00407
00408
00409 h_nchi2[w]->Fill(track->normalizedChi2());
00410 h_nchi2_prob[w]->Fill(TMath::Prob(track->chi2(),(int)track->ndof()));
00411 h_hits[w]->Fill(track->numberOfValidHits());
00412 h_losthits[w]->Fill(track->numberOfLostHits());
00413 chi2_vs_nhits[w]->Fill(track->numberOfValidHits(),track->normalizedChi2());
00414 h_charge[w]->Fill( track->charge() );
00415
00416
00417 edm::ESHandle<MagneticField> theMF;
00418 setup.get<IdealMagneticFieldRecord>().get(theMF);
00419 FreeTrajectoryState
00420 ftsAtProduction(GlobalPoint(tpr->vertex().x(),tpr->vertex().y(),tpr->vertex().z()),
00421 GlobalVector(assocTrack->momentum().x(),assocTrack->momentum().y(),assocTrack->momentum().z()),
00422 TrackCharge(tpr->charge()),
00423 theMF.product());
00424 TrajectoryStateClosestToBeamLineBuilder tscblBuilder;
00425 TrajectoryStateClosestToBeamLine tsAtClosestApproach = tscblBuilder(ftsAtProduction,bs);
00426 GlobalPoint v1 = tsAtClosestApproach.trackStateAtPCA().position();
00427 GlobalVector p = tsAtClosestApproach.trackStateAtPCA().momentum();
00428 GlobalPoint v(v1.x()-bs.x0(),v1.y()-bs.y0(),v1.z()-bs.z0());
00429
00430 double qoverpSim = tsAtClosestApproach.trackStateAtPCA().charge()/p.mag();
00431 double lambdaSim = M_PI/2-p.theta();
00432 double phiSim = p.phi();
00433 double dxySim = (-v.x()*sin(p.phi())+v.y()*cos(p.phi()));
00434 double dzSim = v.z() - (v.x()*p.x()+v.y()*p.y())/p.perp() * p.z()/p.perp();
00435
00436 TrackBase::ParameterVector rParameters = track->parameters();
00437 double qoverpRec = rParameters[0];
00438 double lambdaRec = rParameters[1];
00439 double phiRec = rParameters[2];
00440 double dxyRec = track->dxy(bs.position());
00441 double dzRec = track->dz(bs.position());
00442
00443
00444 double qoverpPull=(qoverpRec-qoverpSim)/track->qoverpError();
00445 double thetaPull=(lambdaRec-lambdaSim)/track->thetaError();
00446 double phiPull=(phiRec-phiSim)/track->phiError();
00447 double dxyPull=(dxyRec-dxySim)/track->dxyError();
00448 double dzPull=(dzRec-dzSim)/track->dzError();
00449
00450 double contrib_Qoverp = ((qoverpRec-qoverpSim)/track->qoverpError())*
00451 ((qoverpRec-qoverpSim)/track->qoverpError())/5;
00452 double contrib_dxy = ((dxyRec-dxySim)/track->dxyError())*((dxyRec-dxySim)/track->dxyError())/5;
00453 double contrib_dz = ((dzRec-dzSim)/track->dzError())*((dzRec-dzSim)/track->dzError())/5;
00454 double contrib_theta = ((lambdaRec-lambdaSim)/track->thetaError())*
00455 ((lambdaRec-lambdaSim)/track->thetaError())/5;
00456 double contrib_phi = ((phiRec-phiSim)/track->phiError())*
00457 ((phiRec-phiSim)/track->phiError())/5;
00458 LogTrace("TrackValidatorTEST") << "assocChi2=" << tp.begin()->second << "\n"
00459 << "" << "\n"
00460 << "ptREC=" << track->pt() << "\n"
00461 << "etaREC=" << track->eta() << "\n"
00462 << "qoverpREC=" << qoverpRec << "\n"
00463 << "dxyREC=" << dxyRec << "\n"
00464 << "dzREC=" << dzRec << "\n"
00465 << "thetaREC=" << track->theta() << "\n"
00466 << "phiREC=" << phiRec << "\n"
00467 << "" << "\n"
00468 << "qoverpError()=" << track->qoverpError() << "\n"
00469 << "dxyError()=" << track->dxyError() << "\n"
00470 << "dzError()=" << track->dzError() << "\n"
00471 << "thetaError()=" << track->thetaError() << "\n"
00472 << "phiError()=" << track->phiError() << "\n"
00473 << "" << "\n"
00474 << "ptSIM=" << sqrt(assocTrack->momentum().perp2()) << "\n"
00475 << "etaSIM=" << assocTrack->momentum().Eta() << "\n"
00476 << "qoverpSIM=" << qoverpSim << "\n"
00477 << "dxySIM=" << dxySim << "\n"
00478 << "dzSIM=" << dzSim << "\n"
00479 << "thetaSIM=" << M_PI/2-lambdaSim << "\n"
00480 << "phiSIM=" << phiSim << "\n"
00481 << "" << "\n"
00482 << "contrib_Qoverp=" << contrib_Qoverp << "\n"
00483 << "contrib_dxy=" << contrib_dxy << "\n"
00484 << "contrib_dz=" << contrib_dz << "\n"
00485 << "contrib_theta=" << contrib_theta << "\n"
00486 << "contrib_phi=" << contrib_phi << "\n"
00487 << "" << "\n"
00488 <<"chi2PULL="<<contrib_Qoverp+contrib_dxy+contrib_dz+contrib_theta+contrib_phi<<"\n";
00489
00490 h_pullQoverp[w]->Fill(qoverpPull);
00491 h_pullTheta[w]->Fill(thetaPull);
00492 h_pullPhi[w]->Fill(phiPull);
00493 h_pullDxy[w]->Fill(dxyPull);
00494 h_pullDz[w]->Fill(dzPull);
00495
00496 double ptres=track->pt()-sqrt(assocTrack->momentum().perp2());
00497 double etares=track->eta()-assocTrack->momentum().Eta();
00498
00499 double ptError = track->ptError();
00500 h_pt[w]->Fill(ptres/ptError);
00501 h_eta[w]->Fill(etares);
00502 etares_vs_eta[w]->Fill(getEta(track->eta()),etares);
00503
00504
00505 chi2_vs_eta[w]->Fill(getEta(track->eta()),track->normalizedChi2());
00506 nhits_vs_eta[w]->Fill(getEta(track->eta()),track->numberOfValidHits());
00507 nlosthits_vs_eta[w]->Fill(getEta(track->eta()),track->numberOfLostHits());
00508
00509
00510 dxyres_vs_eta[w]->Fill(getEta(track->eta()),dxyRec-dxySim);
00511 ptres_vs_eta[w]->Fill(getEta(track->eta()),(track->pt()-sqrt(assocTrack->momentum().perp2()))/track->pt());
00512 dzres_vs_eta[w]->Fill(getEta(track->eta()),dzRec-dzSim);
00513 phires_vs_eta[w]->Fill(getEta(track->eta()),phiRec-phiSim);
00514 cotThetares_vs_eta[w]->Fill(getEta(track->eta()),1/tan(1.570796326794896558-lambdaRec)-1/tan(1.570796326794896558-lambdaSim));
00515
00516
00517 dxyres_vs_pt[w]->Fill(getPt(track->pt()),dxyRec-dxySim);
00518 ptres_vs_pt[w]->Fill(getPt(track->pt()),(track->pt()-sqrt(assocTrack->momentum().perp2()))/track->pt());
00519 dzres_vs_pt[w]->Fill(getPt(track->pt()),dzRec-dzSim);
00520 phires_vs_pt[w]->Fill(getPt(track->pt()),phiRec-phiSim);
00521 cotThetares_vs_pt[w]->Fill(getPt(track->pt()),1/tan(1.570796326794896558-lambdaRec)-1/tan(1.570796326794896558-lambdaSim));
00522
00523
00524 dxypull_vs_eta[w]->Fill(getEta(track->eta()),dxyPull);
00525 ptpull_vs_eta[w]->Fill(getEta(track->eta()),ptres/ptError);
00526 dzpull_vs_eta[w]->Fill(getEta(track->eta()),dzPull);
00527 phipull_vs_eta[w]->Fill(getEta(track->eta()),phiPull);
00528 thetapull_vs_eta[w]->Fill(getEta(track->eta()),thetaPull);
00529
00530
00531 nhits_vs_phi[w]->Fill(track->phi(),track->numberOfValidHits());
00532 chi2_vs_phi[w]->Fill(track->phi(),track->normalizedChi2());
00533 ptmean_vs_eta_phi[w]->Fill(track->phi(),getEta(track->eta()),track->pt());
00534 phimean_vs_eta_phi[w]->Fill(track->phi(),getEta(track->eta()),track->phi());
00535 ptres_vs_phi[w]->Fill(track->phi(),(track->pt()-sqrt(assocTrack->momentum().perp2()))/track->pt());
00536 phires_vs_phi[w]->Fill(track->phi(),phiRec-phiSim);
00537 ptpull_vs_phi[w]->Fill(track->phi(),ptres/ptError);
00538 phipull_vs_phi[w]->Fill(track->phi(),phiPull);
00539 thetapull_vs_phi[w]->Fill(track->phi(),thetaPull);
00540 } catch (cms::Exception e){
00541 LogTrace("TrackValidator") << "exception found: " << e.what() << "\n";
00542 }
00543 }
00544 if (at!=0) h_tracks[w]->Fill(at);
00545 h_fakes[w]->Fill(rT-at);
00546 edm::LogVerbatim("TrackValidator") << "Total Simulated: " << st << "\n"
00547 << "Total Associated (simToReco): " << ats << "\n"
00548 << "Total Reconstructed: " << rT << "\n"
00549 << "Total Associated (recoToSim): " << at << "\n"
00550 << "Total Fakes: " << rT-at << "\n";
00551 nrec_vs_nsim[w]->Fill(rT,st);
00552 w++;
00553 }
00554 }
00555 }
00556
00557 void MultiTrackValidator::endRun(Run const&, EventSetup const&) {
00558
00559 int w=0;
00560 for (unsigned int ww=0;ww<associators.size();ww++){
00561 for (unsigned int www=0;www<label.size();www++){
00562
00563
00564 FitSlicesYTool fsyt_dxy(dxyres_vs_eta[w]);
00565 fsyt_dxy.getFittedSigmaWithError(h_dxyrmsh[w]);
00566 FitSlicesYTool fsyt_dxyPt(dxyres_vs_pt[w]);
00567 fsyt_dxyPt.getFittedSigmaWithError(h_dxyrmshPt[w]);
00568 FitSlicesYTool fsyt_pt(ptres_vs_eta[w]);
00569 fsyt_pt.getFittedSigmaWithError(h_ptrmsh[w]);
00570 fsyt_pt.getFittedMeanWithError(h_ptshifteta[w]);
00571 FitSlicesYTool fsyt_ptPt(ptres_vs_pt[w]);
00572 fsyt_ptPt.getFittedSigmaWithError(h_ptrmshPt[w]);
00573 FitSlicesYTool fsyt_ptPhi(ptres_vs_phi[w]);
00574 fsyt_ptPhi.getFittedSigmaWithError(h_ptrmshPhi[w]);
00575 FitSlicesYTool fsyt_dz(dzres_vs_eta[w]);
00576 fsyt_dz.getFittedSigmaWithError(h_dzrmsh[w]);
00577 FitSlicesYTool fsyt_dzPt(dzres_vs_pt[w]);
00578 fsyt_dzPt.getFittedSigmaWithError(h_dzrmshPt[w]);
00579 FitSlicesYTool fsyt_phi(phires_vs_eta[w]);
00580 fsyt_phi.getFittedSigmaWithError(h_phirmsh[w]);
00581 FitSlicesYTool fsyt_phiPt(phires_vs_pt[w]);
00582 fsyt_phiPt.getFittedSigmaWithError(h_phirmshPt[w]);
00583 FitSlicesYTool fsyt_phiPhi(phires_vs_phi[w]);
00584 fsyt_phiPhi.getFittedSigmaWithError(h_phirmshPhi[w]);
00585 FitSlicesYTool fsyt_cotTheta(cotThetares_vs_eta[w]);
00586 fsyt_cotTheta.getFittedSigmaWithError(h_cotThetarmsh[w]);
00587 FitSlicesYTool fsyt_cotThetaPt(cotThetares_vs_pt[w]);
00588 fsyt_cotThetaPt.getFittedSigmaWithError(h_cotThetarmshPt[w]);
00589
00590
00591 doProfileX(chi2_vs_eta[w],h_chi2meanh[w]);
00592 doProfileX(nhits_vs_eta[w],h_hits_eta[w]);
00593 doProfileX(nlosthits_vs_eta[w],h_losthits_eta[w]);
00594
00595 doProfileX(chi2_vs_nhits[w],h_chi2meanhitsh[w]);
00596 doProfileX(ptres_vs_eta[w],h_ptresmean_vs_eta[w]);
00597 doProfileX(phires_vs_eta[w],h_phiresmean_vs_eta[w]);
00598 doProfileX(chi2_vs_phi[w],h_chi2mean_vs_phi[w]);
00599 doProfileX(nhits_vs_phi[w],h_hits_phi[w]);
00600 doProfileX(ptres_vs_phi[w],h_ptresmean_vs_phi[w]);
00601 doProfileX(phires_vs_phi[w],h_phiresmean_vs_phi[w]);
00602
00603
00604 FitSlicesYTool fsyt_dxyp(dxypull_vs_eta[w]);
00605 fsyt_dxyp.getFittedSigmaWithError(h_dxypulleta[w]);
00606 FitSlicesYTool fsyt_ptp(ptpull_vs_eta[w]);
00607 fsyt_ptp.getFittedSigmaWithError(h_ptpulleta[w]);
00608 FitSlicesYTool fsyt_dzp(dzpull_vs_eta[w]);
00609 fsyt_dzp.getFittedSigmaWithError(h_dzpulleta[w]);
00610 FitSlicesYTool fsyt_phip(phipull_vs_eta[w]);
00611 fsyt_phip.getFittedSigmaWithError(h_phipulleta[w]);
00612 FitSlicesYTool fsyt_thetap(thetapull_vs_eta[w]);
00613 fsyt_thetap.getFittedSigmaWithError(h_thetapulleta[w]);
00614
00615 FitSlicesYTool fsyt_ptpPhi(ptpull_vs_phi[w]);
00616 fsyt_ptpPhi.getFittedSigmaWithError(h_ptpullphi[w]);
00617 FitSlicesYTool fsyt_phipPhi(phipull_vs_phi[w]);
00618 fsyt_phipPhi.getFittedSigmaWithError(h_phipullphi[w]);
00619 FitSlicesYTool fsyt_thetapPhi(thetapull_vs_phi[w]);
00620 fsyt_thetapPhi.getFittedSigmaWithError(h_thetapullphi[w]);
00621
00622
00623 fillPlotFromVectors(h_effic[w],totASSeta[w],totSIMeta[w],"effic");
00624 fillPlotFromVectors(h_fakerate[w],totASS2eta[w],totRECeta[w],"fakerate");
00625 fillPlotFromVectors(h_efficPt[w],totASSpT[w],totSIMpT[w],"effic");
00626 fillPlotFromVectors(h_fakeratePt[w],totASS2pT[w],totRECpT[w],"fakerate");
00627 fillPlotFromVectors(h_effic_vs_hit[w],totASS_hit[w],totSIM_hit[w],"effic");
00628 fillPlotFromVectors(h_fake_vs_hit[w],totASS2_hit[w],totREC_hit[w],"fakerate");
00629
00630 fillPlotFromVector(h_recoeta[w],totRECeta[w]);
00631 fillPlotFromVector(h_simuleta[w],totSIMeta[w]);
00632 fillPlotFromVector(h_assoceta[w],totASSeta[w]);
00633 fillPlotFromVector(h_assoc2eta[w],totASS2eta[w]);
00634
00635 fillPlotFromVector(h_recopT[w],totRECpT[w]);
00636 fillPlotFromVector(h_simulpT[w],totSIMpT[w]);
00637 fillPlotFromVector(h_assocpT[w],totASSpT[w]);
00638 fillPlotFromVector(h_assoc2pT[w],totASS2pT[w]);
00639
00640 fillPlotFromVector(h_recohit[w],totREC_hit[w]);
00641 fillPlotFromVector(h_simulhit[w],totSIM_hit[w]);
00642 fillPlotFromVector(h_assochit[w],totASS_hit[w]);
00643 fillPlotFromVector(h_assoc2hit[w],totASS2_hit[w]);
00644 w++;
00645 }
00646 }
00647 if ( out.size() != 0 && dbe_ ) dbe_->save(out);
00648 }
00649
1.5.4