| 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 |
< |
|
| 53 |
> |
#include <TROOT.h> |
| 54 |
> |
#include <TF1.h> |
| 55 |
> |
#include <TString.h> |
| 56 |
> |
#include <TStyle.h> |
| 57 |
> |
#include <TPaveStats.h> |
| 58 |
> |
#include <TPad.h> |
| 59 |
|
|
| 60 |
|
using namespace std; |
| 61 |
|
|
| 81 |
|
ftree_->Branch("nTrk",&fntrk,"fntrk/I"); |
| 82 |
|
|
| 83 |
|
edm::Service<TFileService> fs; |
| 84 |
< |
// TFileDirectory subDir = fs->mkdir( "Summary" ); |
| 85 |
< |
// TFileDirectory subSubDir = subDir.mkdir( "mySubSubDirectory" ); |
| 86 |
< |
|
| 87 |
< |
h_nTrk = fs->make<TH1F>("nTrk","Num of total tracks",300,0,300); |
| 88 |
< |
h_nTrkPV = fs->make<TH1F>("nTrkPV","Num of Tracks from PV",300,0,300); |
| 89 |
< |
h_trkPt = fs->make<TH1D>("trkPt","Pt of all tracks",100,0,100); |
| 90 |
< |
h_trkPtPV = fs->make<TH1D>("trkPtPV","Pt dist of tracks from PV",100,0,100); |
| 84 |
> |
TFileDirectory subDir = fs->mkdir( "Summary" ); |
| 85 |
> |
TFileDirectory subDir1 = fs->mkdir( "Others" ); |
| 86 |
> |
// TFileDirectory subSubDir = subDir.mkdir( "mySubSubDirectory" ); |
| 87 |
> |
|
| 88 |
> |
setRootStyle(); |
| 89 |
> |
|
| 90 |
> |
h["nTrk"] = subDir.make<TH1F>("nTrk","Num of total rec tracks",300,0,300); |
| 91 |
> |
h["nTrkPV"] = subDir.make<TH1F>("nTrkPV","Num of rec tracks of PVs",300,0,300); |
| 92 |
> |
h["trkPt"] = subDir.make<TH1D>("trkPt","Pt of all rec tracks",100,0,100); |
| 93 |
> |
h["trkPtPV"] = subDir.make<TH1D>("trkPtPV","Pt dist of rec tracks from PV",100,0,100); |
| 94 |
> |
h["genPart_T"] = subDir.make<TH1D>("genPart_T","t component of gen particles",300,-0.5,0.5); |
| 95 |
> |
h["genPart_T"]->GetXaxis()->SetTitle("t (nanosecond)"); |
| 96 |
> |
h["genPart_cT"] = subDir.make<TH1D>("genPart_cT","ct component of gen particles",300,-150.,150.); |
| 97 |
> |
h["genPart_cT"]->GetXaxis()->SetTitle("ct (mm)"); |
| 98 |
|
|
| 99 |
|
if (!realData_) { |
| 100 |
< |
h_deltax = fs->make<TH1D>("deltax","x-distance (Rec - Sim)",800,-0.04,0.04); |
| 101 |
< |
h_deltay = fs->make<TH1D>("deltay","y-distance (Rec - Sim)",800,-0.04,0.04); |
| 102 |
< |
h_deltaz = fs->make<TH1D>("deltaz","z-distance (Rec - Sim)",800,-0.04,0.04); |
| 103 |
< |
h_pullx = fs->make<TH1D>("pullx","x-pull (Rec - Sim)",1000,-25.,25.); |
| 104 |
< |
h_pully = fs->make<TH1D>("pully","y-pull (Rec - Sim)",1000,-25.,25.); |
| 105 |
< |
h_pullz = fs->make<TH1D>("pullz","z-pull (Rec - Sim)",1000,-25.,25.); |
| 100 |
> |
h["deltax"] = subDir.make<TH1D>("deltax","x-distance (Rec - Sim)",800,-0.04,0.04); |
| 101 |
> |
h["deltay"] = subDir.make<TH1D>("deltay","y-distance (Rec - Sim)",800,-0.04,0.04); |
| 102 |
> |
h["deltaz"] = subDir.make<TH1D>("deltaz","z-distance (Rec - Sim)",800,-0.04,0.04); |
| 103 |
> |
h["pullx"] = subDir.make<TH1D>("pullx","x-pull (Rec - Sim)",1000,-25.,25.); |
| 104 |
> |
h["pully"] = subDir.make<TH1D>("pully","y-pull (Rec - Sim)",1000,-25.,25.); |
| 105 |
> |
h["pullz"] = subDir.make<TH1D>("pullz","z-pull (Rec - Sim)",1000,-25.,25.); |
| 106 |
> |
h["errPVx"] = subDir.make<TH1D>("errPVx","X vertex error",100,0.,0.02); |
| 107 |
> |
h["errPVy"] = subDir.make<TH1D>("errPVy","Y vertex error",100,0.,0.02); |
| 108 |
> |
h["errPVz"] = subDir.make<TH1D>("errPVz","Z vertex error",100,0.,0.02); |
| 109 |
|
|
| 110 |
|
// Summary of analysis: |
| 111 |
|
if (analyze_) { |
| 112 |
< |
h_resx_nTrk = fs->make<TH2D>("resx_ntrk","x-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200); |
| 113 |
< |
h_resx_nTrk->SetMarkerStyle(21); |
| 114 |
< |
h_resx_nTrk->SetMarkerColor(4); |
| 115 |
< |
h_resx_nTrk->GetXaxis()->SetTitle("Num of tracks"); |
| 116 |
< |
h_resx_nTrk->GetYaxis()->SetTitle("#mum"); |
| 117 |
< |
h_resy_nTrk = fs->make<TH2D>("resy_ntrk","y-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200); |
| 118 |
< |
h_resy_nTrk->SetMarkerStyle(21); |
| 119 |
< |
h_resy_nTrk->SetMarkerColor(4); |
| 120 |
< |
h_resy_nTrk->GetXaxis()->SetTitle("Num of tracks"); |
| 121 |
< |
h_resy_nTrk->GetYaxis()->SetTitle("#mum"); |
| 122 |
< |
h_resz_nTrk = fs->make<TH2D>("resz_ntrk","z-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200); |
| 123 |
< |
h_resz_nTrk->SetMarkerStyle(21); |
| 124 |
< |
h_resz_nTrk->SetMarkerColor(4); |
| 125 |
< |
h_resz_nTrk->GetXaxis()->SetTitle("Num of tracks"); |
| 126 |
< |
h_resz_nTrk->GetYaxis()->SetTitle("#mum"); |
| 127 |
< |
h_pullx_nTrk = fs->make<TH2D>("pullx_ntrk","x-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.); |
| 128 |
< |
h_pullx_nTrk->SetMarkerStyle(21); |
| 129 |
< |
h_pullx_nTrk->SetMarkerColor(4); |
| 130 |
< |
h_pullx_nTrk->SetBit(TH1::kCanRebin); |
| 131 |
< |
h_pullx_nTrk->GetXaxis()->SetTitle("Num of tracks"); |
| 132 |
< |
h_pully_nTrk = fs->make<TH2D>("pully_ntrk","y-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.); |
| 133 |
< |
h_pully_nTrk->SetMarkerStyle(21); |
| 134 |
< |
h_pully_nTrk->SetMarkerColor(4); |
| 135 |
< |
h_pully_nTrk->SetBit(TH1::kCanRebin); |
| 136 |
< |
h_pully_nTrk->GetXaxis()->SetTitle("Num of tracks"); |
| 137 |
< |
h_pullz_nTrk = fs->make<TH2D>("pullz_ntrk","z-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.); |
| 138 |
< |
h_pullz_nTrk->SetMarkerStyle(21); |
| 139 |
< |
h_pullz_nTrk->SetMarkerColor(4); |
| 140 |
< |
h_pullz_nTrk->SetBit(TH1::kCanRebin); |
| 141 |
< |
h_pullz_nTrk->GetXaxis()->SetTitle("Num of tracks"); |
| 112 |
> |
h2["resx_nTrk"]= subDir.make<TH2D>("resx_ntrk","x-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200); |
| 113 |
> |
h2["resx_nTrk"]->SetMarkerStyle(21); |
| 114 |
> |
h2["resx_nTrk"]->SetMarkerColor(4); |
| 115 |
> |
h2["resx_nTrk"]->GetXaxis()->SetTitle("Num of tracks"); |
| 116 |
> |
h2["resx_nTrk"]->GetYaxis()->SetTitle("#mum"); |
| 117 |
> |
h2["resy_nTrk"]= subDir.make<TH2D>("resy_ntrk","y-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200); |
| 118 |
> |
h2["resy_nTrk"]->SetMarkerStyle(21); |
| 119 |
> |
h2["resy_nTrk"]->SetMarkerColor(4); |
| 120 |
> |
h2["resy_nTrk"]->GetXaxis()->SetTitle("Num of tracks"); |
| 121 |
> |
h2["resy_nTrk"]->GetYaxis()->SetTitle("#mum"); |
| 122 |
> |
h2["resz_nTrk"]= subDir.make<TH2D>("resz_ntrk","z-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200); |
| 123 |
> |
h2["resz_nTrk"]->SetMarkerStyle(21); |
| 124 |
> |
h2["resz_nTrk"]->SetMarkerColor(4); |
| 125 |
> |
h2["resz_nTrk"]->GetXaxis()->SetTitle("Num of tracks"); |
| 126 |
> |
h2["resz_nTrk"]->GetYaxis()->SetTitle("#mum"); |
| 127 |
> |
h2["pullx_nTrk"]= subDir.make<TH2D>("pullx_ntrk","x-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.); |
| 128 |
> |
h2["pullx_nTrk"]->SetMarkerStyle(21); |
| 129 |
> |
h2["pullx_nTrk"]->SetMarkerColor(4); |
| 130 |
> |
h2["pullx_nTrk"]->SetBit(TH1::kCanRebin); |
| 131 |
> |
h2["pullx_nTrk"]->GetXaxis()->SetTitle("Num of tracks"); |
| 132 |
> |
h2["pully_nTrk"]= subDir.make<TH2D>("pully_ntrk","y-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.); |
| 133 |
> |
h2["pully_nTrk"]->SetMarkerStyle(21); |
| 134 |
> |
h2["pully_nTrk"]->SetMarkerColor(4); |
| 135 |
> |
h2["pully_nTrk"]->SetBit(TH1::kCanRebin); |
| 136 |
> |
h2["pully_nTrk"]->GetXaxis()->SetTitle("Num of tracks"); |
| 137 |
> |
h2["pullz_nTrk"]= subDir.make<TH2D>("pullz_ntrk","z-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.); |
| 138 |
> |
h2["pullz_nTrk"]->SetMarkerStyle(21); |
| 139 |
> |
h2["pullz_nTrk"]->SetMarkerColor(4); |
| 140 |
> |
h2["pullz_nTrk"]->SetBit(TH1::kCanRebin); |
| 141 |
> |
h2["pullz_nTrk"]->GetXaxis()->SetTitle("Num of tracks"); |
| 142 |
|
} |
| 143 |
|
for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) { |
| 144 |
|
stringstream ss; |
| 145 |
|
ss << ntrk; |
| 146 |
|
string suffix = ss.str(); |
| 147 |
< |
h["resx_"+suffix] = fs->make<TH1D> (TString("deltax_"+suffix),"x-distance (Rec - Sim)",100,-0.02,0.02); |
| 147 |
> |
h["resx_"+suffix] = fs->make<TH1D>(TString("deltax_"+suffix),"x-distance (Rec - Sim)",100,-0.02,0.02); |
| 148 |
|
h["resx_"+suffix]->GetXaxis()->SetTitle("cm"); |
| 149 |
< |
h["resy_"+suffix] = fs->make<TH1D> (TString("deltay_"+suffix),"y-distance (Rec - Sim)",100,-0.02,0.02); |
| 149 |
> |
h["resy_"+suffix] = fs->make<TH1D>(TString("deltay_"+suffix),"y-distance (Rec - Sim)",100,-0.02,0.02); |
| 150 |
|
h["resy_"+suffix]->GetXaxis()->SetTitle("cm"); |
| 151 |
< |
h["resz_"+suffix] = fs->make<TH1D> (TString("deltaz_"+suffix),"z-distance (Rec - Sim)",100,-0.02,0.02); |
| 151 |
> |
h["resz_"+suffix] = fs->make<TH1D>(TString("deltaz_"+suffix),"z-distance (Rec - Sim)",100,-0.02,0.02); |
| 152 |
|
h["resz_"+suffix]->GetXaxis()->SetTitle("cm"); |
| 153 |
< |
h["pullx_"+suffix] = fs->make<TH1D> (TString("pullx_"+suffix),"x-pull (Rec - Sim)",100,-10.,10.); |
| 154 |
< |
h["pully_"+suffix] = fs->make<TH1D> (TString("pully_"+suffix),"y-pull (Rec - Sim)",100,-10.,10.); |
| 155 |
< |
h["pullz_"+suffix] = fs->make<TH1D> (TString("pullz_"+suffix),"z-pull (Rec - Sim)",100,-10.,10.); |
| 153 |
> |
h["pullx_"+suffix] = fs->make<TH1D>(TString("pullx_"+suffix),"x-pull (Rec - Sim)",100,-10.,10.); |
| 154 |
> |
h["pully_"+suffix] = fs->make<TH1D>(TString("pully_"+suffix),"y-pull (Rec - Sim)",100,-10.,10.); |
| 155 |
> |
h["pullz_"+suffix] = fs->make<TH1D>(TString("pullz_"+suffix),"z-pull (Rec - Sim)",100,-10.,10.); |
| 156 |
> |
h["errPVx_"+suffix] = subDir1.make<TH1D>(TString("errPVx_"+suffix),"X vertex error",100,0.,0.02); |
| 157 |
> |
h["errPVy_"+suffix] = subDir1.make<TH1D>(TString("errPVy_"+suffix),"Y vertex error",100,0.,0.02); |
| 158 |
> |
h["errPVz_"+suffix] = subDir1.make<TH1D>(TString("errPVz_"+suffix),"Z vertex error",100,0.,0.02); |
| 159 |
|
suffix.clear(); |
| 160 |
|
} |
| 161 |
|
} |
| 170 |
|
|
| 171 |
|
} |
| 172 |
|
|
| 173 |
< |
|
| 173 |
> |
void PVStudy::setRootStyle() { |
| 174 |
> |
// |
| 175 |
> |
gROOT->SetStyle("Plain"); |
| 176 |
> |
gStyle->SetFillColor(1); |
| 177 |
> |
gStyle->SetOptDate(); |
| 178 |
> |
// gStyle->SetOptStat(1111110); |
| 179 |
> |
// gStyle->SetOptFit(0111); |
| 180 |
> |
// gStyle->SetPadTickX(1); |
| 181 |
> |
// gStyle->SetPadTickY(1); |
| 182 |
> |
gStyle->SetMarkerSize(0.5); |
| 183 |
> |
gStyle->SetMarkerStyle(8); |
| 184 |
> |
gStyle->SetGridStyle(3); |
| 185 |
> |
//gStyle->SetPadGridX(1); |
| 186 |
> |
//gStyle->SetPadGridY(1); |
| 187 |
> |
gStyle->SetPaperSize(TStyle::kA4); |
| 188 |
> |
gStyle->SetStatW(0.25); // width of statistics box; default is 0.19 |
| 189 |
> |
// gStyle->SetStatH(0.10); // height of statistics box; default is 0.1 |
| 190 |
> |
gStyle->SetStatFormat("6.4g"); // leave default format for now |
| 191 |
> |
gStyle->SetTitleSize(0.055, ""); // size for pad title; default is 0.02 |
| 192 |
> |
gStyle->SetLabelSize(0.03, "XYZ"); // size for axis labels; default is 0.04 |
| 193 |
> |
gStyle->SetStatFontSize(0.08); // size for stat. box |
| 194 |
> |
gStyle->SetTitleFont(32, "XYZ"); // times-bold-italic font (p. 153) for axes |
| 195 |
> |
gStyle->SetTitleFont(32, ""); // same for pad title |
| 196 |
> |
gStyle->SetLabelFont(32, "XYZ"); // same for axis labels |
| 197 |
> |
gStyle->SetStatFont(32); // same for stat. box |
| 198 |
> |
gStyle->SetLabelOffset(0.006, "Y"); // default is 0.005 |
| 199 |
> |
// |
| 200 |
> |
return; |
| 201 |
> |
} |
| 202 |
|
// |
| 203 |
|
// member functions |
| 204 |
|
// |
| 214 |
|
cout << "number of vertices " << evt->vertices_size() << endl; |
| 215 |
|
} |
| 216 |
|
|
| 217 |
< |
int idx=0; |
| 217 |
> |
int idx=0; int npv=0; |
| 218 |
|
for(HepMC::GenEvent::vertex_const_iterator vitr= evt->vertices_begin(); |
| 219 |
< |
vitr != evt->vertices_end(); ++vitr ) |
| 220 |
< |
{ // loop for vertex ... |
| 221 |
< |
HepMC::FourVector pos = (*vitr)->position(); |
| 222 |
< |
//HepLorentzVector pos = (*vitr)->position(); |
| 223 |
< |
if (pos.t()>0) { continue;} // ? |
| 224 |
< |
|
| 225 |
< |
bool hasMotherVertex=false; |
| 226 |
< |
// if(verbose_) cout << "mothers" << endl; |
| 227 |
< |
for ( HepMC::GenVertex::particle_iterator mother = (*vitr)->particles_begin(HepMC::parents); |
| 228 |
< |
mother != (*vitr)->particles_end(HepMC::parents); ++mother ) { |
| 229 |
< |
HepMC::GenVertex * mv=(*mother)->production_vertex(); |
| 230 |
< |
if (mv) {hasMotherVertex=true;} |
| 231 |
< |
// if(verbose_) { |
| 232 |
< |
// cout << "\t"; |
| 191 |
< |
// (*mother)->print(); |
| 192 |
< |
// } |
| 193 |
< |
} |
| 194 |
< |
|
| 195 |
< |
if(hasMotherVertex) {continue;} |
| 196 |
< |
|
| 197 |
< |
// could be a new vertex, check all primaries found so far to avoid multiple entries |
| 198 |
< |
const double mm2cm=0.1; |
| 199 |
< |
simPrimaryVertex sv(pos.x()*mm2cm,pos.y()*mm2cm,pos.z()*mm2cm); // sim unit mm, rec unit cm |
| 200 |
< |
simPrimaryVertex *vp=NULL; // will become non-NULL if a vertex is found and then point to it |
| 201 |
< |
for(vector<simPrimaryVertex>::iterator v0=simpv.begin(); |
| 202 |
< |
v0!=simpv.end(); v0++){ |
| 203 |
< |
if( (fabs(sv.x-v0->x)<1e-5) && (fabs(sv.y-v0->y)<1e-5) && (fabs(sv.z-v0->z)<1e-5)){ |
| 204 |
< |
vp=&(*v0); |
| 205 |
< |
break; |
| 219 |
> |
vitr != evt->vertices_end(); ++vitr ) { // loop for vertex ... |
| 220 |
> |
HepMC::FourVector pos = (*vitr)->position(); |
| 221 |
> |
//HepLorentzVector pos = (*vitr)->position(); |
| 222 |
> |
|
| 223 |
> |
// t component of PV: |
| 224 |
> |
for ( HepMC::GenVertex::particle_iterator mother = (*vitr)->particles_begin(HepMC::parents); |
| 225 |
> |
mother != (*vitr)->particles_end(HepMC::parents); ++mother ) { |
| 226 |
> |
// std::cout << "Status = " << (*mother)->status() << std::endl; |
| 227 |
> |
HepMC::GenVertex * mv=(*mother)->production_vertex(); |
| 228 |
> |
if( ((*mother)->status() == 3) && (!mv)) { |
| 229 |
> |
// std::cout << "npv= " << npv << std::endl; |
| 230 |
> |
if (npv == 0) { |
| 231 |
> |
h["genPart_cT"]->Fill(pos.t()); // mm |
| 232 |
> |
h["genPart_T"]->Fill(pos.t()/299.792458); // ns |
| 233 |
|
} |
| 234 |
+ |
npv++; |
| 235 |
|
} |
| 236 |
+ |
} |
| 237 |
+ |
// if (pos.t()>0) { continue;} // for 22X when t of PV was not smeared |
| 238 |
+ |
|
| 239 |
+ |
bool hasMotherVertex=false; |
| 240 |
+ |
if(verbose_) cout << "mothers of vertex[" << ++idx << "]: " << endl; |
| 241 |
+ |
for ( HepMC::GenVertex::particle_iterator mother = (*vitr)->particles_begin(HepMC::parents); |
| 242 |
+ |
mother != (*vitr)->particles_end(HepMC::parents); ++mother ) { |
| 243 |
+ |
HepMC::GenVertex * mv=(*mother)->production_vertex(); |
| 244 |
+ |
// std::cout << "Status = " << (*mother)->status() << std::endl; |
| 245 |
+ |
if (mv) { |
| 246 |
+ |
hasMotherVertex=true; |
| 247 |
+ |
if(!verbose_) break; //if verbose_, print all particles of gen vertices |
| 248 |
+ |
} |
| 249 |
+ |
if(verbose_) { |
| 250 |
+ |
cout << "\t"; |
| 251 |
+ |
(*mother)->print(); |
| 252 |
+ |
} |
| 253 |
+ |
} |
| 254 |
|
|
| 255 |
< |
if(!vp){ |
| 256 |
< |
// this is a new vertex |
| 257 |
< |
if(verbose_) cout << "this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl; |
| 258 |
< |
simpv.push_back(sv); |
| 259 |
< |
vp=&simpv.back(); |
| 260 |
< |
}else{ |
| 261 |
< |
if(verbose_) cout << "this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl; |
| 262 |
< |
} |
| 263 |
< |
vp->genVertex.push_back((*vitr)->barcode()); |
| 264 |
< |
// collect final state descendants |
| 265 |
< |
for ( HepMC::GenVertex::particle_iterator |
| 266 |
< |
daughter = (*vitr)->particles_begin(HepMC::descendants); |
| 267 |
< |
daughter != (*vitr)->particles_end(HepMC::descendants); |
| 268 |
< |
++daughter ) { |
| 269 |
< |
if (isFinalstateParticle(*daughter)){ |
| 270 |
< |
if ( find(vp->finalstateParticles.begin(), vp->finalstateParticles.end(),(*daughter)->barcode()) |
| 271 |
< |
== vp->finalstateParticles.end()){ |
| 272 |
< |
vp->finalstateParticles.push_back((*daughter)->barcode()); |
| 273 |
< |
HepMC::FourVector m=(*daughter)->momentum(); |
| 274 |
< |
// the next four lines used to be "vp->ptot+=m;" in the days of CLHEP::HepLorentzVector |
| 275 |
< |
// but adding FourVectors seems not to be foreseen |
| 276 |
< |
vp->ptot.setPx(vp->ptot.px()+m.px()); |
| 277 |
< |
vp->ptot.setPy(vp->ptot.py()+m.py()); |
| 278 |
< |
vp->ptot.setPz(vp->ptot.pz()+m.pz()); |
| 279 |
< |
vp->ptot.setE(vp->ptot.e()+m.e()); |
| 280 |
< |
vp->ptsq+=(m.perp())*(m.perp()); |
| 281 |
< |
if ( (m.perp()>0.8) && (fabs(m.pseudoRapidity())<2.5) && isCharged( *daughter ) ){ |
| 282 |
< |
vp->nGenTrk++; |
| 283 |
< |
} |
| 255 |
> |
if(hasMotherVertex) continue; |
| 256 |
> |
|
| 257 |
> |
// could be a new vertex, check all primaries found so far to avoid multiple entries |
| 258 |
> |
const double mm2cm=0.1; |
| 259 |
> |
simPrimaryVertex sv(pos.x()*mm2cm,pos.y()*mm2cm,pos.z()*mm2cm); // sim unit mm, rec unit cm |
| 260 |
> |
simPrimaryVertex *vp=NULL; // will become non-NULL if a vertex is found and then point to it |
| 261 |
> |
for(vector<simPrimaryVertex>::iterator v0=simpv.begin(); |
| 262 |
> |
v0!=simpv.end(); v0++){ |
| 263 |
> |
if( (fabs(sv.x-v0->x)<1e-5) && (fabs(sv.y-v0->y)<1e-5) && (fabs(sv.z-v0->z)<1e-5)){ |
| 264 |
> |
vp=&(*v0); |
| 265 |
> |
break; |
| 266 |
> |
} |
| 267 |
> |
} |
| 268 |
> |
|
| 269 |
> |
if(!vp){ |
| 270 |
> |
// this is a new vertex |
| 271 |
> |
if(verbose_) cout << "this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl; |
| 272 |
> |
simpv.push_back(sv); |
| 273 |
> |
vp=&simpv.back(); |
| 274 |
> |
}else{ |
| 275 |
> |
if(verbose_) cout << "this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl; |
| 276 |
> |
} |
| 277 |
> |
vp->genVertex.push_back((*vitr)->barcode()); |
| 278 |
> |
// collect final state descendants |
| 279 |
> |
for ( HepMC::GenVertex::particle_iterator daughter = (*vitr)->particles_begin(HepMC::descendants); |
| 280 |
> |
daughter != (*vitr)->particles_end(HepMC::descendants); |
| 281 |
> |
++daughter ) { |
| 282 |
> |
if (isFinalstateParticle(*daughter)){ |
| 283 |
> |
if ( find(vp->finalstateParticles.begin(), vp->finalstateParticles.end(),(*daughter)->barcode()) |
| 284 |
> |
== vp->finalstateParticles.end()){ |
| 285 |
> |
vp->finalstateParticles.push_back((*daughter)->barcode()); |
| 286 |
> |
HepMC::FourVector m=(*daughter)->momentum(); |
| 287 |
> |
// the next four lines used to be "vp->ptot+=m;" in the days of CLHEP::HepLorentzVector |
| 288 |
> |
// but adding FourVectors seems not to be foreseen |
| 289 |
> |
vp->ptot.setPx(vp->ptot.px()+m.px()); |
| 290 |
> |
vp->ptot.setPy(vp->ptot.py()+m.py()); |
| 291 |
> |
vp->ptot.setPz(vp->ptot.pz()+m.pz()); |
| 292 |
> |
vp->ptot.setE(vp->ptot.e()+m.e()); |
| 293 |
> |
vp->ptsq+=(m.perp())*(m.perp()); |
| 294 |
> |
if ( (m.perp()>0.8) && (fabs(m.pseudoRapidity())<2.5) && isCharged( *daughter ) ){ |
| 295 |
> |
vp->nGenTrk++; |
| 296 |
|
} |
| 297 |
|
} |
| 298 |
|
} |
| 241 |
– |
idx++; |
| 299 |
|
} |
| 300 |
+ |
} |
| 301 |
|
} |
| 302 |
|
return simpv; |
| 303 |
|
} |
| 461 |
|
ferror[1] = vsim->recVtx->yError(); |
| 462 |
|
ferror[2] = vsim->recVtx->zError(); |
| 463 |
|
fntrk = nrectrk; |
| 406 |
– |
ftree_->Fill(); |
| 464 |
|
|
| 465 |
< |
h_deltax->Fill( vsim->recVtx->x()-vsim->x ); |
| 466 |
< |
h_deltay->Fill( vsim->recVtx->y()-vsim->y ); |
| 467 |
< |
h_deltaz->Fill( vsim->recVtx->z()-vsim->z ); |
| 468 |
< |
h_pullx->Fill( (vsim->recVtx->x()-vsim->x )/vsim->recVtx->xError() ); |
| 469 |
< |
h_pully->Fill( (vsim->recVtx->y()-vsim->y )/vsim->recVtx->yError() ); |
| 470 |
< |
h_pullz->Fill( (vsim->recVtx->z()-vsim->z )/vsim->recVtx->zError() ); |
| 471 |
< |
pvinfo.push_back(PVStudy::PVInfo(res(vsim->recVtx->x()-vsim->x, |
| 472 |
< |
vsim->recVtx->y()-vsim->y, |
| 473 |
< |
vsim->recVtx->z()-vsim->z), |
| 474 |
< |
error(vsim->recVtx->xError(), |
| 475 |
< |
vsim->recVtx->yError(), |
| 419 |
< |
vsim->recVtx->zError()), |
| 465 |
> |
h["deltax"]->Fill( fres[0] ); |
| 466 |
> |
h["deltay"]->Fill( fres[1] ); |
| 467 |
> |
h["deltaz"]->Fill( fres[2] ); |
| 468 |
> |
h["pullx"]->Fill( fres[0]/ferror[0] ); |
| 469 |
> |
h["pully"]->Fill( fres[1]/ferror[1] ); |
| 470 |
> |
h["pullz"]->Fill( fres[2]/ferror[2] ); |
| 471 |
> |
h["errPVx"]->Fill( ferror[0] ); |
| 472 |
> |
h["errPVy"]->Fill( ferror[1] ); |
| 473 |
> |
h["errPVz"]->Fill( ferror[2] ); |
| 474 |
> |
pvinfo.push_back(PVStudy::PVInfo(res(fres[0], fres[1], fres[2]), |
| 475 |
> |
error(ferror[0], ferror[1], ferror[2]), |
| 476 |
|
nrectrk)); |
| 477 |
|
// Fill histo according to its track multiplicity |
| 478 |
< |
fillHisto(res(vsim->recVtx->x()-vsim->x, |
| 479 |
< |
vsim->recVtx->y()-vsim->y, |
| 424 |
< |
vsim->recVtx->z()-vsim->z), |
| 425 |
< |
pull((vsim->recVtx->x()-vsim->x )/vsim->recVtx->xError(), |
| 426 |
< |
(vsim->recVtx->y()-vsim->y )/vsim->recVtx->yError(), |
| 427 |
< |
(vsim->recVtx->z()-vsim->z )/vsim->recVtx->zError()), |
| 478 |
> |
fillHisto(res(fres[0], fres[1], fres[2]), |
| 479 |
> |
error(ferror[0], ferror[1], ferror[2]), |
| 480 |
|
nrectrk); |
| 481 |
+ |
|
| 482 |
+ |
ftree_->Fill(); |
| 483 |
|
} |
| 484 |
|
else{ // no rec vertex found for this simvertex |
| 485 |
|
if(verbose_) { |
| 498 |
|
t!=v->tracks_end(); t++) { |
| 499 |
|
// illegal charge |
| 500 |
|
if ( (**t).charge() < -1 || (**t).charge() > 1 ) { |
| 501 |
< |
// h_trkPtPV->Fill(0.); |
| 501 |
> |
// h["trkPtPV"]->Fill(0.); |
| 502 |
|
} |
| 503 |
|
else { |
| 504 |
< |
h_trkPtPV->Fill((**t).pt()); |
| 504 |
> |
h["trkPtPV"]->Fill((**t).pt()); |
| 505 |
|
} |
| 506 |
|
} |
| 507 |
|
} |
| 508 |
|
catch (...) { |
| 509 |
|
// exception thrown when trying to use linked track |
| 510 |
< |
// h_trkPtPV->Fill(0.); |
| 510 |
> |
// h["trkPtPV"]->Fill(0.); |
| 511 |
|
} |
| 512 |
|
|
| 513 |
< |
h_nTrkPV->Fill(v->tracksSize()); |
| 513 |
> |
h["nTrkPV"]->Fill(v->tracksSize()); |
| 514 |
|
|
| 515 |
|
} |
| 516 |
|
|
| 517 |
< |
h_nTrk->Fill(tracks->size()); |
| 517 |
> |
h["nTrk"]->Fill(tracks->size()); |
| 518 |
|
|
| 519 |
|
for(TrackCollection::const_iterator itTrack = tracks->begin(); |
| 520 |
|
itTrack != tracks->end(); |
| 521 |
|
++itTrack) { |
| 522 |
|
int charge = 0; |
| 523 |
|
charge = itTrack->charge(); |
| 524 |
< |
h_trkPt->Fill(itTrack->pt()); |
| 524 |
> |
h["trkPt"]->Fill(itTrack->pt()); |
| 525 |
|
} |
| 526 |
|
|
| 527 |
|
} |
| 543 |
|
if (verbose_) cout << "Fill point of ntrk = " << ntrk << endl; |
| 544 |
|
PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk); |
| 545 |
|
if (analyze_) { |
| 546 |
< |
// if ( ipv.res_.x() > 0 ) h_resx_nTrk->Fill(ntrk,ipv.res_.x()/sqt2, 1.); |
| 547 |
< |
// if ( ipv.res_.y() > 0 ) h_resy_nTrk->Fill(ntrk,ipv.res_.y()/sqt2, 1.); |
| 548 |
< |
// if ( ipv.res_.z() > 0 ) h_resz_nTrk->Fill(ntrk,ipv.res_.z()/sqt2, 1.); |
| 549 |
< |
if ( ipv.res_.x() > 0 ) h_resx_nTrk->Fill(ntrk,ipv.res_.x(), 1.); |
| 550 |
< |
if ( ipv.res_.y() > 0 ) h_resy_nTrk->Fill(ntrk,ipv.res_.y(), 1.); |
| 551 |
< |
if ( ipv.res_.z() > 0 ) h_resz_nTrk->Fill(ntrk,ipv.res_.z(), 1.); |
| 552 |
< |
if ( ipv.pull_.x() > 0 ) h_pullx_nTrk->Fill(ntrk,ipv.pull_.x(), 1.); |
| 553 |
< |
if ( ipv.pull_.y() > 0 ) h_pully_nTrk->Fill(ntrk,ipv.pull_.y(), 1.); |
| 554 |
< |
if ( ipv.pull_.z() > 0 ) h_pullz_nTrk->Fill(ntrk,ipv.pull_.z(), 1.); |
| 546 |
> |
// if ( ipv.res_.x() > 0 ) h2["resx_nTrk"]->Fill(ntrk,ipv.res_.x()/sqt2, 1.); |
| 547 |
> |
// if ( ipv.res_.y() > 0 ) h2["resy_nTrk"]->Fill(ntrk,ipv.res_.y()/sqt2, 1.); |
| 548 |
> |
// if ( ipv.res_.z() > 0 ) h2["resz_nTrk"]->Fill(ntrk,ipv.res_.z()/sqt2, 1.); |
| 549 |
> |
if ( ipv.res_.x() > 0 ) h2["resx_nTrk"]->Fill(ntrk,ipv.res_.x(), 1.); |
| 550 |
> |
if ( ipv.res_.y() > 0 ) h2["resy_nTrk"]->Fill(ntrk,ipv.res_.y(), 1.); |
| 551 |
> |
if ( ipv.res_.z() > 0 ) h2["resz_nTrk"]->Fill(ntrk,ipv.res_.z(), 1.); |
| 552 |
> |
if ( ipv.pull_.x() > 0 ) h2["pullx_nTrk"]->Fill(ntrk,ipv.pull_.x(), 1.); |
| 553 |
> |
if ( ipv.pull_.y() > 0 ) h2["pully_nTrk"]->Fill(ntrk,ipv.pull_.y(), 1.); |
| 554 |
> |
if ( ipv.pull_.z() > 0 ) h2["pullz_nTrk"]->Fill(ntrk,ipv.pull_.z(), 1.); |
| 555 |
|
} |
| 556 |
|
} |
| 557 |
|
} |
| 637 |
|
} |
| 638 |
|
} |
| 639 |
|
|
| 640 |
< |
void PVStudy::fillHisto(res r, pull p, int ntk){ |
| 640 |
> |
void PVStudy::fillHisto(res r, error er, int ntk){ |
| 641 |
|
stringstream ss; |
| 642 |
|
ss << ntk; |
| 643 |
|
string suffix = ss.str(); |
| 646 |
|
h["resx_"+suffix]->Fill(r.x()); |
| 647 |
|
h["resy_"+suffix]->Fill(r.y()); |
| 648 |
|
h["resz_"+suffix]->Fill(r.z()); |
| 649 |
< |
h["pullx_"+suffix]->Fill(p.x()); |
| 650 |
< |
h["pully_"+suffix]->Fill(p.y()); |
| 651 |
< |
h["pullz_"+suffix]->Fill(p.z()); |
| 649 |
> |
h["pullx_"+suffix]->Fill(r.x()/er.x()); |
| 650 |
> |
h["pully_"+suffix]->Fill(r.y()/er.y()); |
| 651 |
> |
h["pullz_"+suffix]->Fill(r.z()/er.z()); |
| 652 |
> |
h["errPVx_"+suffix]->Fill(er.x()); |
| 653 |
> |
h["errPVy_"+suffix]->Fill(er.y()); |
| 654 |
> |
h["errPVz_"+suffix]->Fill(er.z()); |
| 655 |
|
} |
| 656 |
|
|
| 657 |
|
PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk) { |
| 658 |
|
map<int,double> data; |
| 659 |
|
data.clear(); |
| 660 |
< |
double fpar[2]; |
| 660 |
> |
double fpar[2] = {-999,-999}; |
| 661 |
|
double cm2um = 10000; |
| 662 |
|
stringstream ss; |
| 663 |
|
ss << ntk; |
| 665 |
|
|
| 666 |
|
if(analyze_) { |
| 667 |
|
for ( int i = 0; i < 6; ++i) { |
| 611 |
– |
cout << "Here" << endl; |
| 668 |
|
switch (i) { |
| 669 |
|
case 0: |
| 670 |
|
fit(h["resx_"+suffix], fpar); |
| 706 |
|
ntk); |
| 707 |
|
} |
| 708 |
|
|
| 709 |
< |
void PVStudy::fit(TH1 *&hdist, double data[]){ |
| 709 |
> |
void PVStudy::fit(TH1 *&hdist, double fitPars[]){ |
| 710 |
|
TAxis *axis0 = hdist->GetXaxis(); |
| 711 |
|
int nbin = axis0->GetLast(); |
| 712 |
|
double nOF = hdist->GetBinContent(nbin+1); |
| 713 |
|
double nUF = hdist->GetBinContent(0); |
| 714 |
+ |
// double fitRange = axis0->GetBinUpEdge(nbin); |
| 715 |
+ |
double fitRange = axis0->GetXmax(); |
| 716 |
+ |
double sigMax[2] = {0.,0.}; |
| 717 |
+ |
|
| 718 |
|
if ( verbose_ ){ |
| 719 |
|
cout << "Last bin = " << nbin; |
| 720 |
|
cout << "; hdist: Overflow Entries = " << nOF; |
| 721 |
|
cout << "; Underflow Entries = " << nUF; |
| 722 |
< |
cout << "; hdist->GetEntries() = " << hdist->GetEntries() << endl; |
| 722 |
> |
cout << "; hdist->GetEntries() = " << hdist->GetEntries(); |
| 723 |
> |
cout << "; fitting range = " << -fitRange << " to " << fitRange << endl; |
| 724 |
|
} |
| 725 |
< |
if (hdist->GetEntries() - nOF - nUF >= 5) { // FIXME: for reasonable Gaussian fit |
| 726 |
< |
hdist->Fit("gaus","Q"); |
| 727 |
< |
TF1 *fgaus = hdist->GetFunction("gaus"); |
| 725 |
> |
if (hdist->GetEntries() - nOF - nUF >= 10) { // FIXME: for reasonable Gaussian fit |
| 726 |
> |
for (int bi = 0; bi < nbin; bi++) { |
| 727 |
> |
if ( (axis0->GetBinLowEdge(bi) < 0) && (hdist->GetBinContent(bi) > 0) ) { |
| 728 |
> |
sigMax[0] = axis0->GetBinLowEdge(bi); |
| 729 |
> |
if ( abs(sigMax[0]) > abs(sigMax[1]) ) sigMax[1] = abs(sigMax[0]); |
| 730 |
> |
} |
| 731 |
> |
if ( (axis0->GetBinLowEdge(bi) >= 0) && (hdist->GetBinContent(bi) > 0) ) { |
| 732 |
> |
sigMax[0] = axis0->GetBinUpEdge(bi); |
| 733 |
> |
if ( abs(sigMax[0]) > abs(sigMax[1]) ) sigMax[1] = abs(sigMax[0]); |
| 734 |
> |
} |
| 735 |
> |
} |
| 736 |
> |
if (verbose_) cout << "Fit sigMax = " << sqrt(2.)*sigMax[1] << endl; |
| 737 |
> |
|
| 738 |
> |
TF1 *fgaus = new TF1("fgaus","gaus",-fitRange, fitRange); |
| 739 |
> |
fgaus->SetParameter(1, 0.); |
| 740 |
> |
fgaus->SetParLimits(1, -fitRange/40., fitRange/40.); |
| 741 |
> |
fgaus->SetParLimits(2, 0., sqrt(2.)*sigMax[1]); |
| 742 |
> |
fgaus->SetLineColor(4); |
| 743 |
> |
hdist->Fit(fgaus,"Q"); |
| 744 |
> |
gPad->Update(); |
| 745 |
> |
TPaveStats *s = (TPaveStats*) hdist->GetListOfFunctions()->FindObject("stats"); |
| 746 |
> |
s->SetOptStat(1111111); |
| 747 |
> |
s->SetOptFit(0111); |
| 748 |
> |
gPad->Update(); |
| 749 |
|
if (verbose_) cout << "got function" << endl; |
| 750 |
< |
data[0] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999); |
| 750 |
> |
fitPars[0] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999); |
| 751 |
|
} |
| 752 |
+ |
else |
| 753 |
+ |
fitPars[0] = -999; |
| 754 |
|
} |
| 755 |
|
|
| 756 |
|
//define this as a plug-in |
