| 1 |
// -*- C++ -*-
|
| 2 |
//
|
| 3 |
// Package: PVStudy
|
| 4 |
// Class: PVStudy
|
| 5 |
//
|
| 6 |
/**\class PVStudy PVStudy.cc UserCode/PVStudy/plugins/PVStudy.cc
|
| 7 |
|
| 8 |
Description: <one line class summary>
|
| 9 |
|
| 10 |
Implementation:
|
| 11 |
<Notes on implementation>
|
| 12 |
*/
|
| 13 |
//
|
| 14 |
// Original Author: "Geng-yuan Jeng/UC Riverside"
|
| 15 |
// Created: Thu Aug 20 11:55:40 CDT 2009
|
| 16 |
// $Id$
|
| 17 |
//
|
| 18 |
//
|
| 19 |
|
| 20 |
|
| 21 |
// system include files
|
| 22 |
#include <memory>
|
| 23 |
#include <string>
|
| 24 |
#include <vector>
|
| 25 |
#include <iostream>
|
| 26 |
#include <sstream>
|
| 27 |
|
| 28 |
// user include files
|
| 29 |
#include "FWCore/Framework/interface/Frameworkfwd.h"
|
| 30 |
#include "FWCore/Framework/interface/EDAnalyzer.h"
|
| 31 |
|
| 32 |
#include "FWCore/Framework/interface/Event.h"
|
| 33 |
#include "FWCore/Framework/interface/MakerMacros.h"
|
| 34 |
|
| 35 |
#include "FWCore/ParameterSet/interface/ParameterSet.h"
|
| 36 |
#include "FWCore/Utilities/interface/InputTag.h"
|
| 37 |
#include "DataFormats/TrackReco/interface/Track.h"
|
| 38 |
#include "DataFormats/TrackReco/interface/TrackFwd.h"
|
| 39 |
#include "FWCore/ServiceRegistry/interface/Service.h"
|
| 40 |
#include "PhysicsTools/UtilAlgos/interface/TFileService.h"
|
| 41 |
#include "UserCode/PVStudy/interface/PVStudy.h"
|
| 42 |
//
|
| 43 |
#include "DataFormats/VertexReco/interface/Vertex.h"
|
| 44 |
#include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
|
| 45 |
|
| 46 |
// simulated vertices,..., add <use name=SimDataFormats/Vertex> and <../Track>
|
| 47 |
#include <SimDataFormats/Vertex/interface/SimVertex.h>
|
| 48 |
#include <SimDataFormats/Vertex/interface/SimVertexContainer.h>
|
| 49 |
#include <SimDataFormats/Track/interface/SimTrack.h>
|
| 50 |
#include <SimDataFormats/Track/interface/SimTrackContainer.h>
|
| 51 |
|
| 52 |
//root
|
| 53 |
#include "TROOT.h"
|
| 54 |
#include "TF1.h"
|
| 55 |
#include "TString.h"
|
| 56 |
#include "TStyle.h"
|
| 57 |
|
| 58 |
|
| 59 |
using namespace std;
|
| 60 |
|
| 61 |
PVStudy::PVStudy(const edm::ParameterSet& iConfig):
|
| 62 |
nTrkMin_(10),nTrkMax_(150)
|
| 63 |
{
|
| 64 |
simG4_ = iConfig.getParameter<edm::InputTag>( "simG4" );
|
| 65 |
tracksLabel_ = iConfig.getUntrackedParameter<edm::InputTag>("tracks");
|
| 66 |
vtxSample_ = iConfig.getUntrackedParameter<edm::InputTag>("vtxSample");
|
| 67 |
verbose_ = iConfig.getUntrackedParameter<bool>("verbose", false);
|
| 68 |
realData_ = iConfig.getUntrackedParameter<bool>("realData",false);
|
| 69 |
|
| 70 |
//now do what ever initialization is needed
|
| 71 |
pvanainfo.clear();
|
| 72 |
|
| 73 |
edm::Service<TFileService> fs;
|
| 74 |
// TFileDirectory subDir = fs->mkdir( "Summary" );
|
| 75 |
// TFileDirectory subSubDir = subDir.mkdir( "mySubSubDirectory" );
|
| 76 |
|
| 77 |
h_nTrk = fs->make<TH1F>("nTrk","Num of total tracks",300,0,300);
|
| 78 |
h_nTrkPV = fs->make<TH1F>("nTrkPV","Num of Tracks from PV",300,0,300);
|
| 79 |
h_trkPt = fs->make<TH1D>("trkPt","Pt of all tracks",100,0,100);
|
| 80 |
h_trkPtPV = fs->make<TH1D>("trkPtPV","Pt dist of tracks from PV",100,0,100);
|
| 81 |
|
| 82 |
if (!realData_) {
|
| 83 |
h_deltax = fs->make<TH1D>("deltax","x-distance (Rec - Sim)",800,-0.04,0.04);
|
| 84 |
h_deltay = fs->make<TH1D>("deltay","y-distance (Rec - Sim)",800,-0.04,0.04);
|
| 85 |
h_deltaz = fs->make<TH1D>("deltaz","z-distance (Rec - Sim)",800,-0.04,0.04);
|
| 86 |
h_pullx = fs->make<TH1D>("pullx","x-pull (Rec - Sim)",1000,-25.,25.);
|
| 87 |
h_pully = fs->make<TH1D>("pully","y-pull (Rec - Sim)",1000,-25.,25.);
|
| 88 |
h_pullz = fs->make<TH1D>("pullz","z-pull (Rec - Sim)",1000,-25.,25.);
|
| 89 |
// Summary of analysis:
|
| 90 |
h_resx_nTrk = fs->make<TH2D>("resx_ntrk","x-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
|
| 91 |
h_resx_nTrk->SetMarkerStyle(21);
|
| 92 |
h_resx_nTrk->GetYaxis()->SetTitle("#mum");
|
| 93 |
h_resy_nTrk = fs->make<TH2D>("resy_ntrk","y-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
|
| 94 |
h_resy_nTrk->SetMarkerStyle(21);
|
| 95 |
h_resy_nTrk->GetYaxis()->SetTitle("#mum");
|
| 96 |
h_resz_nTrk = fs->make<TH2D>("resz_ntrk","z-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
|
| 97 |
h_resz_nTrk->SetMarkerStyle(21);
|
| 98 |
h_resz_nTrk->GetYaxis()->SetTitle("#mum");
|
| 99 |
h_pullx_nTrk = fs->make<TH2D>("pullx_ntrk","x-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
|
| 100 |
h_pullx_nTrk->SetMarkerStyle(21);
|
| 101 |
h_pullx_nTrk->SetBit(TH1::kCanRebin);
|
| 102 |
h_pullx_nTrk->GetYaxis()->SetTitle("cm");
|
| 103 |
h_pully_nTrk = fs->make<TH2D>("pully_ntrk","y-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
|
| 104 |
h_pully_nTrk->SetMarkerStyle(21);
|
| 105 |
h_pully_nTrk->SetBit(TH1::kCanRebin);
|
| 106 |
h_pully_nTrk->GetYaxis()->SetTitle("cm");
|
| 107 |
h_pullz_nTrk = fs->make<TH2D>("pullz_ntrk","z-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
|
| 108 |
h_pullz_nTrk->SetMarkerStyle(21);
|
| 109 |
h_pullz_nTrk->SetBit(TH1::kCanRebin);
|
| 110 |
h_pullz_nTrk->GetYaxis()->SetTitle("cm");
|
| 111 |
|
| 112 |
for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
|
| 113 |
stringstream ss;
|
| 114 |
ss << ntrk;
|
| 115 |
string suffix = ss.str();
|
| 116 |
h["resx_"+suffix] = fs->make<TH1D> (TString("deltax_"+suffix),"x-distance (Rec - Sim)",400,-0.02,0.02);
|
| 117 |
h["resx_"+suffix]->GetXaxis()->SetTitle("cm");
|
| 118 |
h["resy_"+suffix] = fs->make<TH1D> (TString("deltay_"+suffix),"y-distance (Rec - Sim)",400,-0.02,0.02);
|
| 119 |
h["resy_"+suffix]->GetXaxis()->SetTitle("cm");
|
| 120 |
h["resz_"+suffix] = fs->make<TH1D> (TString("deltaz_"+suffix),"z-distance (Rec - Sim)",400,-0.02,0.02);
|
| 121 |
h["resz_"+suffix]->GetXaxis()->SetTitle("cm");
|
| 122 |
h["pullx_"+suffix] = fs->make<TH1D> (TString("pullx_"+suffix),"x-pull (Rec - Sim)",200,-10.,10.);
|
| 123 |
h["pullx_"+suffix]->GetXaxis()->SetTitle("cm");
|
| 124 |
h["pully_"+suffix] = fs->make<TH1D> (TString("pully_"+suffix),"y-pull (Rec - Sim)",200,-10.,10.);
|
| 125 |
h["pully_"+suffix]->GetXaxis()->SetTitle("cm");
|
| 126 |
h["pullz_"+suffix] = fs->make<TH1D> (TString("pullz_"+suffix),"z-pull (Rec - Sim)",200,-10.,10.);
|
| 127 |
h["pullz_"+suffix]->GetXaxis()->SetTitle("cm");
|
| 128 |
suffix.clear();
|
| 129 |
}
|
| 130 |
}
|
| 131 |
}
|
| 132 |
|
| 133 |
|
| 134 |
PVStudy::~PVStudy()
|
| 135 |
{
|
| 136 |
|
| 137 |
// do anything here that needs to be done at desctruction time
|
| 138 |
// (e.g. close files, deallocate resources etc.)
|
| 139 |
|
| 140 |
}
|
| 141 |
|
| 142 |
|
| 143 |
//
|
| 144 |
// member functions
|
| 145 |
//
|
| 146 |
std::vector<PVStudy::simPrimaryVertex> PVStudy::getSimPVs(const Handle<HepMCProduct> evtMC, std::string suffix="")
|
| 147 |
{
|
| 148 |
std::vector<PVStudy::simPrimaryVertex> simpv;
|
| 149 |
const HepMC::GenEvent* evt=evtMC->GetEvent();
|
| 150 |
if (evt) {
|
| 151 |
if(verbose_) {
|
| 152 |
cout << "process id " << evt->signal_process_id()<<endl;
|
| 153 |
cout << "signal process vertex " << ( evt->signal_process_vertex() ?
|
| 154 |
evt->signal_process_vertex()->barcode() : 0 ) <<endl;
|
| 155 |
cout << "number of vertices " << evt->vertices_size() << endl;
|
| 156 |
}
|
| 157 |
|
| 158 |
int idx=0;
|
| 159 |
for(HepMC::GenEvent::vertex_const_iterator vitr= evt->vertices_begin();
|
| 160 |
vitr != evt->vertices_end(); ++vitr )
|
| 161 |
{ // loop for vertex ...
|
| 162 |
HepMC::FourVector pos = (*vitr)->position();
|
| 163 |
//HepLorentzVector pos = (*vitr)->position();
|
| 164 |
if (pos.t()>0) { continue;} // ?
|
| 165 |
|
| 166 |
bool hasMotherVertex=false;
|
| 167 |
// if(verbose_) cout << "mothers" << endl;
|
| 168 |
for ( HepMC::GenVertex::particle_iterator mother = (*vitr)->particles_begin(HepMC::parents);
|
| 169 |
mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
|
| 170 |
HepMC::GenVertex * mv=(*mother)->production_vertex();
|
| 171 |
if (mv) {hasMotherVertex=true;}
|
| 172 |
// if(verbose_) {
|
| 173 |
// cout << "\t";
|
| 174 |
// (*mother)->print();
|
| 175 |
// }
|
| 176 |
}
|
| 177 |
|
| 178 |
if(hasMotherVertex) {continue;}
|
| 179 |
|
| 180 |
// could be a new vertex, check all primaries found so far to avoid multiple entries
|
| 181 |
const double mm2cm=0.1;
|
| 182 |
simPrimaryVertex sv(pos.x()*mm2cm,pos.y()*mm2cm,pos.z()*mm2cm); // sim unit mm, rec unit cm
|
| 183 |
simPrimaryVertex *vp=NULL; // will become non-NULL if a vertex is found and then point to it
|
| 184 |
for(vector<simPrimaryVertex>::iterator v0=simpv.begin();
|
| 185 |
v0!=simpv.end(); v0++){
|
| 186 |
if( (fabs(sv.x-v0->x)<1e-5) && (fabs(sv.y-v0->y)<1e-5) && (fabs(sv.z-v0->z)<1e-5)){
|
| 187 |
vp=&(*v0);
|
| 188 |
break;
|
| 189 |
}
|
| 190 |
}
|
| 191 |
|
| 192 |
if(!vp){
|
| 193 |
// this is a new vertex
|
| 194 |
if(verbose_) cout << "this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
|
| 195 |
simpv.push_back(sv);
|
| 196 |
vp=&simpv.back();
|
| 197 |
}else{
|
| 198 |
if(verbose_) cout << "this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
|
| 199 |
}
|
| 200 |
vp->genVertex.push_back((*vitr)->barcode());
|
| 201 |
// collect final state descendants
|
| 202 |
for ( HepMC::GenVertex::particle_iterator
|
| 203 |
daughter = (*vitr)->particles_begin(HepMC::descendants);
|
| 204 |
daughter != (*vitr)->particles_end(HepMC::descendants);
|
| 205 |
++daughter ) {
|
| 206 |
if (isFinalstateParticle(*daughter)){
|
| 207 |
if ( find(vp->finalstateParticles.begin(), vp->finalstateParticles.end(),(*daughter)->barcode())
|
| 208 |
== vp->finalstateParticles.end()){
|
| 209 |
vp->finalstateParticles.push_back((*daughter)->barcode());
|
| 210 |
HepMC::FourVector m=(*daughter)->momentum();
|
| 211 |
// the next four lines used to be "vp->ptot+=m;" in the days of CLHEP::HepLorentzVector
|
| 212 |
// but adding FourVectors seems not to be foreseen
|
| 213 |
vp->ptot.setPx(vp->ptot.px()+m.px());
|
| 214 |
vp->ptot.setPy(vp->ptot.py()+m.py());
|
| 215 |
vp->ptot.setPz(vp->ptot.pz()+m.pz());
|
| 216 |
vp->ptot.setE(vp->ptot.e()+m.e());
|
| 217 |
vp->ptsq+=(m.perp())*(m.perp());
|
| 218 |
if ( (m.perp()>0.8) && (fabs(m.pseudoRapidity())<2.5) && isCharged( *daughter ) ){
|
| 219 |
vp->nGenTrk++;
|
| 220 |
}
|
| 221 |
}
|
| 222 |
}
|
| 223 |
}
|
| 224 |
idx++;
|
| 225 |
}
|
| 226 |
}
|
| 227 |
return simpv;
|
| 228 |
}
|
| 229 |
|
| 230 |
std::vector<PVStudy::simPrimaryVertex> PVStudy::getSimPVs(const Handle<HepMCProduct> evtMC,
|
| 231 |
const Handle<SimVertexContainer> simVtxs,
|
| 232 |
const Handle<SimTrackContainer> simTrks)
|
| 233 |
{
|
| 234 |
// simvertices don't have enough information to decide,
|
| 235 |
// genVertices don't have the simulated coordinates ( with VtxSmeared they might)
|
| 236 |
// go through simtracks to get the link between simulated and generated vertices
|
| 237 |
std::vector<PVStudy::simPrimaryVertex> simpv;
|
| 238 |
int idx=0;
|
| 239 |
for(SimTrackContainer::const_iterator t=simTrks->begin();
|
| 240 |
t!=simTrks->end(); ++t){
|
| 241 |
if ( !(t->noVertex()) && !(t->type()==-99) ){
|
| 242 |
double ptsq=0;
|
| 243 |
bool primary=false; // something coming directly from the primary vertex
|
| 244 |
bool resonance=false; // resonance
|
| 245 |
bool track=false; // undecayed, charged particle
|
| 246 |
HepMC::GenParticle* gp=evtMC->GetEvent()->barcode_to_particle( (*t).genpartIndex() );
|
| 247 |
if (gp) {
|
| 248 |
HepMC::GenVertex * gv=gp->production_vertex();
|
| 249 |
if (gv) {
|
| 250 |
for ( HepMC::GenVertex::particle_iterator
|
| 251 |
daughter = gv->particles_begin(HepMC::descendants);
|
| 252 |
daughter != gv->particles_end(HepMC::descendants);
|
| 253 |
++daughter ) {
|
| 254 |
if (isFinalstateParticle(*daughter)){
|
| 255 |
ptsq+=(*daughter)->momentum().perp()*(*daughter)->momentum().perp();
|
| 256 |
}
|
| 257 |
}
|
| 258 |
primary = ( gv->position().t()==0);
|
| 259 |
//resonance= ( gp->mother() && isResonance(gp->mother())); // in CLHEP/HepMC days
|
| 260 |
// no more mother pointer in the improved HepMC GenParticle
|
| 261 |
resonance= ( isResonance(*(gp->production_vertex()->particles_in_const_begin())));
|
| 262 |
if (gp->status()==1){
|
| 263 |
//track=((pdt->particle(gp->pdg_id()))->charge() != 0);
|
| 264 |
track=not isCharged(gp);
|
| 265 |
}
|
| 266 |
}
|
| 267 |
}
|
| 268 |
|
| 269 |
const HepMC::FourVector & v=(*simVtxs)[t->vertIndex()].position();
|
| 270 |
//const HepLorentzVector & v=(*simVtxs)[t->vertIndex()].position();
|
| 271 |
if(primary or resonance){
|
| 272 |
{
|
| 273 |
// check all primaries found so far to avoid multiple entries
|
| 274 |
bool newVertex=true;
|
| 275 |
for(std::vector<PVStudy::simPrimaryVertex>::iterator v0=simpv.begin();
|
| 276 |
v0!=simpv.end(); v0++){
|
| 277 |
if( (fabs(v0->x-v.x())<0.001) && (fabs(v0->y-v.y())<0.001) && (fabs(v0->z-v.z())<0.001) ){
|
| 278 |
if (track) {
|
| 279 |
v0->simTrackIndex.push_back(idx);
|
| 280 |
if (ptsq>(*v0).ptsq){(*v0).ptsq=ptsq;}
|
| 281 |
}
|
| 282 |
newVertex=false;
|
| 283 |
}
|
| 284 |
}
|
| 285 |
if(newVertex && !resonance){
|
| 286 |
PVStudy::simPrimaryVertex anotherVertex(v.x(),v.y(),v.z());
|
| 287 |
if (track) anotherVertex.simTrackIndex.push_back(idx);
|
| 288 |
anotherVertex.ptsq=ptsq;
|
| 289 |
simpv.push_back(anotherVertex);
|
| 290 |
}
|
| 291 |
}//
|
| 292 |
}
|
| 293 |
|
| 294 |
}// simtrack has vertex and valid type
|
| 295 |
idx++;
|
| 296 |
}//simTrack loop
|
| 297 |
return simpv;
|
| 298 |
}
|
| 299 |
|
| 300 |
// ------------ method called to for each event ------------
|
| 301 |
void
|
| 302 |
PVStudy::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
|
| 303 |
{
|
| 304 |
using namespace edm;
|
| 305 |
|
| 306 |
using reco::TrackCollection;
|
| 307 |
|
| 308 |
Handle<TrackCollection> tracks;
|
| 309 |
iEvent.getByLabel(tracksLabel_,tracks);
|
| 310 |
|
| 311 |
Handle<reco::VertexCollection> recVtxs;
|
| 312 |
iEvent.getByLabel(vtxSample_, recVtxs);
|
| 313 |
|
| 314 |
if (!realData_) {
|
| 315 |
Handle<SimVertexContainer> simVtxs;
|
| 316 |
iEvent.getByLabel( simG4_, simVtxs);
|
| 317 |
|
| 318 |
Handle<SimTrackContainer> simTrks;
|
| 319 |
iEvent.getByLabel( simG4_, simTrks);
|
| 320 |
}
|
| 321 |
|
| 322 |
if(verbose_) printRecVtxs(recVtxs);
|
| 323 |
|
| 324 |
try{
|
| 325 |
iSetup.getData(pdt);
|
| 326 |
}catch(const Exception&){
|
| 327 |
cout << "Some problem occurred with the particle data table. This may not work !." <<endl;
|
| 328 |
}
|
| 329 |
|
| 330 |
if (!realData_) {
|
| 331 |
bool MC=false;
|
| 332 |
Handle<HepMCProduct> evtMC;
|
| 333 |
iEvent.getByLabel("generator",evtMC);
|
| 334 |
if (!evtMC.isValid()) {
|
| 335 |
MC=false;
|
| 336 |
if(verbose_){
|
| 337 |
cout << "no HepMCProduct found"<< endl;
|
| 338 |
}
|
| 339 |
} else {
|
| 340 |
if(verbose_){
|
| 341 |
cout << "generator HepMCProduct found"<< endl;
|
| 342 |
}
|
| 343 |
MC=true;
|
| 344 |
}
|
| 345 |
|
| 346 |
if(MC){
|
| 347 |
// make a list of primary vertices:
|
| 348 |
std::vector<simPrimaryVertex> simpv;
|
| 349 |
simpv=getSimPVs(evtMC,"");
|
| 350 |
// simpv=getSimPVs(evtMC, simVtxs, simTrks);
|
| 351 |
|
| 352 |
int nsimtrk=0;
|
| 353 |
int nrectrk=0;
|
| 354 |
for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
|
| 355 |
vsim!=simpv.end(); vsim++){
|
| 356 |
nsimtrk+=vsim->nGenTrk;
|
| 357 |
// look for a matching reconstructed vertex
|
| 358 |
vsim->recVtx=NULL;
|
| 359 |
for(reco::VertexCollection::const_iterator vrec=recVtxs->begin();
|
| 360 |
vrec!=recVtxs->end(); ++vrec){
|
| 361 |
nrectrk=vrec->tracksSize();
|
| 362 |
if(verbose_){
|
| 363 |
cout << "sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z << endl;
|
| 364 |
cout << "Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
|
| 365 |
}
|
| 366 |
if ( matchVertex(*vsim,*vrec) ){
|
| 367 |
// if the matching critera are fulfilled, accept the rec-vertex that is closest in z
|
| 368 |
if( ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
|
| 369 |
|| (!vsim->recVtx) )
|
| 370 |
{
|
| 371 |
vsim->recVtx=&(*vrec);
|
| 372 |
}
|
| 373 |
}
|
| 374 |
}
|
| 375 |
if (vsim->recVtx){
|
| 376 |
if(verbose_) {
|
| 377 |
cout <<"primary matched " << vsim->x << " " << vsim->y << " " << vsim->z << endl;
|
| 378 |
}
|
| 379 |
h_deltax->Fill( vsim->recVtx->x()-vsim->x );
|
| 380 |
h_deltay->Fill( vsim->recVtx->y()-vsim->y );
|
| 381 |
h_deltaz->Fill( vsim->recVtx->z()-vsim->z );
|
| 382 |
h_pullx->Fill( (vsim->recVtx->x()-vsim->x )/vsim->recVtx->xError() );
|
| 383 |
h_pully->Fill( (vsim->recVtx->y()-vsim->y )/vsim->recVtx->yError() );
|
| 384 |
h_pullz->Fill( (vsim->recVtx->z()-vsim->z )/vsim->recVtx->zError() );
|
| 385 |
pvanainfo.push_back(PVStudy::PVAnaInfo( vsim->recVtx->x()-vsim->x,
|
| 386 |
vsim->recVtx->y()-vsim->y,
|
| 387 |
vsim->recVtx->z()-vsim->z,
|
| 388 |
(vsim->recVtx->x()-vsim->x )/vsim->recVtx->xError(),
|
| 389 |
(vsim->recVtx->y()-vsim->y )/vsim->recVtx->yError(),
|
| 390 |
(vsim->recVtx->z()-vsim->z )/vsim->recVtx->zError(),
|
| 391 |
nrectrk )
|
| 392 |
);
|
| 393 |
}
|
| 394 |
else{ // no rec vertex found for this simvertex
|
| 395 |
if(verbose_) {
|
| 396 |
cout <<"primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
|
| 397 |
}
|
| 398 |
}
|
| 399 |
}
|
| 400 |
}
|
| 401 |
}// With MC
|
| 402 |
|
| 403 |
// Loop RecVtxs and find matched SimVtxs
|
| 404 |
for(reco::VertexCollection::const_iterator v=recVtxs->begin();
|
| 405 |
v!=recVtxs->end(); ++v){
|
| 406 |
try {
|
| 407 |
for(reco::Vertex::trackRef_iterator t = v->tracks_begin();
|
| 408 |
t!=v->tracks_end(); t++) {
|
| 409 |
// illegal charge
|
| 410 |
if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
|
| 411 |
// h_trkPtPV->Fill(0.);
|
| 412 |
}
|
| 413 |
else {
|
| 414 |
h_trkPtPV->Fill((**t).pt());
|
| 415 |
}
|
| 416 |
}
|
| 417 |
}
|
| 418 |
catch (...) {
|
| 419 |
// exception thrown when trying to use linked track
|
| 420 |
// h_trkPtPV->Fill(0.);
|
| 421 |
}
|
| 422 |
|
| 423 |
h_nTrkPV->Fill(v->tracksSize());
|
| 424 |
|
| 425 |
}
|
| 426 |
|
| 427 |
h_nTrk->Fill(tracks->size());
|
| 428 |
|
| 429 |
for(TrackCollection::const_iterator itTrack = tracks->begin();
|
| 430 |
itTrack != tracks->end();
|
| 431 |
++itTrack) {
|
| 432 |
int charge = 0;
|
| 433 |
charge = itTrack->charge();
|
| 434 |
h_trkPt->Fill(itTrack->pt());
|
| 435 |
}
|
| 436 |
|
| 437 |
}
|
| 438 |
|
| 439 |
|
| 440 |
// ------------ method called once each job just before starting event loop ------------
|
| 441 |
void
|
| 442 |
PVStudy::beginJob()
|
| 443 |
{
|
| 444 |
}
|
| 445 |
|
| 446 |
// ------------ method called once each job just after ending the event loop ------------
|
| 447 |
void
|
| 448 |
PVStudy::endJob() {
|
| 449 |
double sqt2 = sqrt(2.);
|
| 450 |
if (!realData_) {
|
| 451 |
if (verbose_) cout << "Analyzing PV info" << endl;
|
| 452 |
for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
|
| 453 |
if (verbose_) cout << "Fill point of ntrk = " << ntrk << endl;
|
| 454 |
PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk);
|
| 455 |
if ( ipv.resx > 0 ) h_resx_nTrk->Fill(ntrk,ipv.resx/sqt2, 1.);
|
| 456 |
if ( ipv.resy > 0 ) h_resy_nTrk->Fill(ntrk,ipv.resy/sqt2, 1.);
|
| 457 |
if ( ipv.resz > 0 ) h_resz_nTrk->Fill(ntrk,ipv.resz/sqt2, 1.);
|
| 458 |
if ( ipv.pullx > 0 ) h_pullx_nTrk->Fill(ntrk,ipv.pullx, 1.);
|
| 459 |
if ( ipv.pully > 0 ) h_pully_nTrk->Fill(ntrk,ipv.pully, 1.);
|
| 460 |
if ( ipv.pullz > 0 ) h_pullz_nTrk->Fill(ntrk,ipv.pullz, 1.);
|
| 461 |
}
|
| 462 |
}
|
| 463 |
}
|
| 464 |
|
| 465 |
bool PVStudy::matchVertex(const PVStudy::simPrimaryVertex & vsim,
|
| 466 |
const reco::Vertex & vrec)
|
| 467 |
{
|
| 468 |
return (fabs(vsim.z-vrec.z())<0.0500); // =500um
|
| 469 |
}
|
| 470 |
|
| 471 |
bool PVStudy::isResonance(const HepMC::GenParticle * p){
|
| 472 |
double ctau=(pdt->particle( abs(p->pdg_id()) ))->lifetime();
|
| 473 |
if(verbose_) cout << "isResonance " << p->pdg_id() << " " << ctau << endl;
|
| 474 |
return ctau >0 && ctau <1e-6;
|
| 475 |
}
|
| 476 |
|
| 477 |
bool PVStudy::isFinalstateParticle(const HepMC::GenParticle * p){
|
| 478 |
return ( !p->end_vertex() && p->status()==1 );
|
| 479 |
}
|
| 480 |
|
| 481 |
bool PVStudy::isCharged(const HepMC::GenParticle * p){
|
| 482 |
const ParticleData * part = pdt->particle( p->pdg_id() );
|
| 483 |
if (part){
|
| 484 |
return part->charge()!=0;
|
| 485 |
}else{
|
| 486 |
// the new/improved particle table doesn't know anti-particles
|
| 487 |
return pdt->particle( -p->pdg_id() )!=0;
|
| 488 |
}
|
| 489 |
}
|
| 490 |
|
| 491 |
void PVStudy::printSimVtxs(const Handle<SimVertexContainer> simVtxs){
|
| 492 |
int i=0;
|
| 493 |
for(SimVertexContainer::const_iterator vsim=simVtxs->begin();
|
| 494 |
vsim!=simVtxs->end(); ++vsim){
|
| 495 |
cout << i++ << ")" << scientific
|
| 496 |
<< " evtid=" << vsim->eventId().event()
|
| 497 |
<< " sim x=" << vsim->position().x()
|
| 498 |
<< " sim y=" << vsim->position().y()
|
| 499 |
<< " sim z=" << vsim->position().z()
|
| 500 |
<< " sim t=" << vsim->position().t()
|
| 501 |
<< " parent=" << vsim->parentIndex()
|
| 502 |
<< endl;
|
| 503 |
}
|
| 504 |
}
|
| 505 |
|
| 506 |
void PVStudy::printRecVtxs(const Handle<reco::VertexCollection> recVtxs){
|
| 507 |
int ivtx=0;
|
| 508 |
for(reco::VertexCollection::const_iterator v=recVtxs->begin();
|
| 509 |
v!=recVtxs->end(); ++v){
|
| 510 |
cout << "Recvtx "<< std::setw(3) << std::setfill(' ')<<ivtx++
|
| 511 |
<< "#trk " << std::setw(3) << v->tracksSize()
|
| 512 |
<< " chi2 " << std::setw(4) << v->chi2()
|
| 513 |
<< " ndof " << std::setw(3) << v->ndof() << endl
|
| 514 |
<< " x " << std::setw(8) <<std::fixed << std::setprecision(4) << v->x()
|
| 515 |
<< " dx " << std::setw(8) << v->xError()<< endl
|
| 516 |
<< " y " << std::setw(8) << v->y()
|
| 517 |
<< " dy " << std::setw(8) << v->yError()<< endl
|
| 518 |
<< " z " << std::setw(8) << v->z()
|
| 519 |
<< " dz " << std::setw(8) << v->zError()
|
| 520 |
<< endl;
|
| 521 |
}
|
| 522 |
}
|
| 523 |
|
| 524 |
void PVStudy::printSimTrks(const Handle<SimTrackContainer> simTrks){
|
| 525 |
cout << " simTrks type, (momentum), vertIndex, genpartIndex" << endl;
|
| 526 |
int i=1;
|
| 527 |
for(SimTrackContainer::const_iterator t=simTrks->begin();
|
| 528 |
t!=simTrks->end(); ++t){
|
| 529 |
//HepMC::GenParticle* gp=evtMC->GetEvent()->particle( (*t).genpartIndex() );
|
| 530 |
cout << i++ << ")"
|
| 531 |
<< (*t)
|
| 532 |
<< " index="
|
| 533 |
<< (*t).genpartIndex();
|
| 534 |
//if (gp) {
|
| 535 |
// HepMC::GenVertex *gv=gp->production_vertex();
|
| 536 |
// cout << " genvertex =" << (*gv);
|
| 537 |
//}
|
| 538 |
cout << endl;
|
| 539 |
}
|
| 540 |
}
|
| 541 |
|
| 542 |
PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk) {
|
| 543 |
map<int,double> data;
|
| 544 |
double xmin=0.,xmax=0.;
|
| 545 |
double fpar[2];
|
| 546 |
double cm2um = 10000;
|
| 547 |
stringstream ss;
|
| 548 |
ss << ntk;
|
| 549 |
string suffix = ss.str();
|
| 550 |
|
| 551 |
for ( int i = 0; i < 6; ++i) {
|
| 552 |
switch(i){
|
| 553 |
case 0:
|
| 554 |
case 1:
|
| 555 |
case 2:
|
| 556 |
if (verbose_) cout << "Analyzing Resolution:" << endl;
|
| 557 |
xmin = -0.01; xmax = 0.01;
|
| 558 |
break;
|
| 559 |
default:
|
| 560 |
if (verbose_) cout << "Analyzing Pull:" << endl;
|
| 561 |
xmin = -10.; xmax = 10;
|
| 562 |
break;
|
| 563 |
}
|
| 564 |
TH1D *hdist = new TH1D("hdist","distribution to be fitted", 1000,xmin,xmax);
|
| 565 |
|
| 566 |
for (std::vector<PVStudy::PVAnaInfo>::const_iterator ipvinfo = pvanainfo.begin();
|
| 567 |
ipvinfo != pvanainfo.end(); ++ipvinfo){
|
| 568 |
if ( verbose_ && i==0 ) {
|
| 569 |
cout << "Num of tracks of the PV = " << ipvinfo->nTrk;
|
| 570 |
cout << "; delta x = " << ipvinfo->resx << endl;
|
| 571 |
}
|
| 572 |
if (ipvinfo->nTrk == ntk) {
|
| 573 |
switch(i) {
|
| 574 |
case 0: if (verbose_) cout << "Matched resx" << endl;
|
| 575 |
hdist->Fill(ipvinfo->resx);
|
| 576 |
h["resx_"+suffix]->Fill(ipvinfo->resx);
|
| 577 |
break;
|
| 578 |
case 1: if (verbose_) cout << "Matched resy" << endl;
|
| 579 |
hdist->Fill(ipvinfo->resy);
|
| 580 |
h["resy_"+suffix]->Fill(ipvinfo->resy);
|
| 581 |
break;
|
| 582 |
case 2: if (verbose_) cout << "Matched resz" << endl;
|
| 583 |
hdist->Fill(ipvinfo->resz);
|
| 584 |
h["resz_"+suffix]->Fill(ipvinfo->resz);
|
| 585 |
break;
|
| 586 |
case 3: if (verbose_) cout << "Matched pullx" << endl;
|
| 587 |
hdist->Fill(ipvinfo->pullx);
|
| 588 |
h["pullx_"+suffix]->Fill(ipvinfo->pullx);
|
| 589 |
break;
|
| 590 |
case 4: if (verbose_) cout << "Matched pully" << endl;
|
| 591 |
hdist->Fill(ipvinfo->pully);
|
| 592 |
h["pully_"+suffix]->Fill(ipvinfo->pully);
|
| 593 |
break;
|
| 594 |
case 5: if (verbose_) cout << "Matched pullz" << endl;
|
| 595 |
hdist->Fill(ipvinfo->pullz);
|
| 596 |
h["pullz_"+suffix]->Fill(ipvinfo->pullz);
|
| 597 |
break;
|
| 598 |
}
|
| 599 |
}
|
| 600 |
else {
|
| 601 |
// switch(i) {
|
| 602 |
// case 0: if (verbose_) cout << "bad resx" << endl; break;
|
| 603 |
// case 1: if (verbose_) cout << "bad resy" << endl; break;
|
| 604 |
// case 2: if (verbose_) cout << "bad resz" << endl; break;
|
| 605 |
// case 3: if (verbose_) cout << "bad pullx" << endl; break;
|
| 606 |
// case 4: if (verbose_) cout << "bad pully" << endl; break;
|
| 607 |
// case 5: if (verbose_) cout << "bad pullz" << endl; break;
|
| 608 |
// }
|
| 609 |
}
|
| 610 |
}
|
| 611 |
TAxis *axis0 = hdist->GetXaxis();
|
| 612 |
int nbin = axis0->GetLast();
|
| 613 |
double nOF = hdist->GetBinContent(nbin+1);
|
| 614 |
double nUF = hdist->GetBinContent(0);
|
| 615 |
if ( verbose_ && i==0 ){
|
| 616 |
cout << "Last bin = " << nbin;
|
| 617 |
cout << "; hdist: Overflow Entries = " << nOF;
|
| 618 |
cout << "; Underflow Entries = " << nUF;
|
| 619 |
cout << "; hdist->GetEntries() = " << hdist->GetEntries() << endl;
|
| 620 |
}
|
| 621 |
if (hdist->GetEntries() - nOF - nUF >= 5) { // FIXME: for reasonable Gaussian fit
|
| 622 |
if(i<3){
|
| 623 |
hdist->Fit("gaus","Q0");
|
| 624 |
TF1 *fgaus = hdist->GetFunction("gaus");
|
| 625 |
if (verbose_) cout << "got function" << endl;
|
| 626 |
fpar[0] = ((fgaus->GetParameter(1))?fgaus->GetParameter(1):-999);
|
| 627 |
fpar[1] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999);
|
| 628 |
delete fgaus;
|
| 629 |
}
|
| 630 |
else {
|
| 631 |
hdist->Fit("gausn","Q0");
|
| 632 |
TF1 *fgaus = hdist->GetFunction("gausn");
|
| 633 |
if (verbose_) cout << "got function" << endl;
|
| 634 |
fpar[0] = ((fgaus->GetParameter(1))?fgaus->GetParameter(1):-999);
|
| 635 |
fpar[1] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999);
|
| 636 |
delete fgaus;
|
| 637 |
}
|
| 638 |
if (verbose_) cout << "fpar[2] = (" << fpar[0] << "," << fpar[1] << ")" << endl;
|
| 639 |
switch (i) {
|
| 640 |
case 0:
|
| 641 |
h["resx_"+suffix]->Fit("gaus","Q");
|
| 642 |
data[i] = fpar[1]*cm2um;
|
| 643 |
break;
|
| 644 |
case 1:
|
| 645 |
h["resy_"+suffix]->Fit("gaus","Q");
|
| 646 |
data[i] = fpar[1]*cm2um;
|
| 647 |
break;
|
| 648 |
case 2:
|
| 649 |
h["resz_"+suffix]->Fit("gaus","Q");
|
| 650 |
data[i] = fpar[1]*cm2um;
|
| 651 |
break;
|
| 652 |
case 3:
|
| 653 |
h["pullx_"+suffix]->Fit("gausn","Q");
|
| 654 |
data[i] = fpar[1];
|
| 655 |
break;
|
| 656 |
case 4:
|
| 657 |
h["pully_"+suffix]->Fit("gausn","Q");
|
| 658 |
data[i] = fpar[1];
|
| 659 |
break;
|
| 660 |
case 5:
|
| 661 |
h["pullz_"+suffix]->Fit("gausn","Q");
|
| 662 |
data[i] = fpar[1];
|
| 663 |
break;
|
| 664 |
}
|
| 665 |
if (verbose_ && data[i] > 0) cout << "data[" << i << "] = " << data[i] << (i<3?" micro m":" cm") << endl;
|
| 666 |
}
|
| 667 |
else {
|
| 668 |
if ( verbose_) cout << "Empty Histo! No Fit!!!" << endl;
|
| 669 |
data[i] = -999;
|
| 670 |
}
|
| 671 |
delete hdist;
|
| 672 |
}
|
| 673 |
|
| 674 |
return PVStudy::PVAnaInfo(data[0], data[1], data[2],
|
| 675 |
data[3], data[4], data[5],
|
| 676 |
ntk);
|
| 677 |
}
|
| 678 |
|
| 679 |
//define this as a plug-in
|
| 680 |
DEFINE_FWK_MODULE(PVStudy);
|
