| 1 |
// $Id: VertexTools.cc,v 1.9 2012/03/22 15:54:07 bendavid Exp $
|
| 2 |
|
| 3 |
#include "MitPhysics/Utils/interface/VertexTools.h"
|
| 4 |
#include "MitPhysics/Utils/interface/ElectronTools.h"
|
| 5 |
#include "MitPhysics/Utils/interface/PhotonTools.h"
|
| 6 |
#include "MitAna/DataTree/interface/StableData.h"
|
| 7 |
#include <TFile.h>
|
| 8 |
#include <TVector3.h>
|
| 9 |
#include <TSystem.h>
|
| 10 |
|
| 11 |
ClassImp(mithep::VertexTools)
|
| 12 |
|
| 13 |
using namespace mithep;
|
| 14 |
|
| 15 |
VertexTools* VertexTools::meobject = NULL;
|
| 16 |
|
| 17 |
//--------------------------------------------------------------------------------------------------
|
| 18 |
VertexTools::VertexTools() :
|
| 19 |
fIsInitMvaM(kFALSE),
|
| 20 |
fIsInitMvaP(kFALSE),
|
| 21 |
reader(0),
|
| 22 |
readervtx(0),
|
| 23 |
readerevt(0)
|
| 24 |
{
|
| 25 |
|
| 26 |
}
|
| 27 |
|
| 28 |
//--------------------------------------------------------------------------------------------------
|
| 29 |
void VertexTools::InitM(const char* str)
|
| 30 |
{
|
| 31 |
|
| 32 |
relname = str;
|
| 33 |
reader = new TMVA::Reader( "!Color:!Silent" );
|
| 34 |
reader->AddVariable( "var1", &tmvar1 );
|
| 35 |
reader->AddVariable( "var2", &tmvar2 );
|
| 36 |
reader->AddVariable( "var3", &tmvar3 );
|
| 37 |
reader->AddVariable( "var4", &tmvar4 );
|
| 38 |
reader->AddVariable( "var5", &tmvar5 );
|
| 39 |
reader->AddVariable( "var6", &tmvar6 );
|
| 40 |
reader->BookMVA( "BDTG method",relname + TString("/src/MitPhysics/data/TMVAClassification_BDTG.weights.xml").Data());
|
| 41 |
reader->BookMVA( "BDTD method",relname + TString("/src/MitPhysics/data/TMVAClassification_BDTD.weights.xml" ).Data());
|
| 42 |
//reader->BookMVA( "CFMlpANN method", "/home/maxi/cms/root/TMVAClassification_CFMlpANN.weights.xml" );
|
| 43 |
reader->BookMVA( "MLP method", relname + TString("/src/MitPhysics/data/TMVAClassification_MLP.weights.xml").Data());
|
| 44 |
reader->BookMVA( "MLPBFGS method",relname + TString("/src/MitPhysics/data/TMVAClassification_MLPBFGS.weights.xml" ).Data());
|
| 45 |
|
| 46 |
fIsInitMvaM = kTRUE;
|
| 47 |
|
| 48 |
}
|
| 49 |
|
| 50 |
//--------------------------------------------------------------------------------------------------
|
| 51 |
void VertexTools::InitP()
|
| 52 |
{
|
| 53 |
|
| 54 |
readervtx = new TMVA::Reader( "!Color:!Silent" );
|
| 55 |
readerevt = new TMVA::Reader( "!Color:!Silent" );
|
| 56 |
|
| 57 |
readervtx->AddVariable( "ptbal", &fMvaPVars[0] );
|
| 58 |
readervtx->AddVariable( "ptasym", &fMvaPVars[1] );
|
| 59 |
readervtx->AddVariable( "logsumpt2", &fMvaPVars[2] );
|
| 60 |
readervtx->AddVariable( "limPullToConv", &fMvaPVars[3] );
|
| 61 |
readervtx->AddVariable( "nConv", &fMvaPVars[4] );
|
| 62 |
readervtx->BookMVA( "BDTCat", gSystem->Getenv("CMSSW_BASE") + TString("/src/MitPhysics/data/TMVAClassification_BDTCat_conversions_tmva_407.weights.xml") );
|
| 63 |
|
| 64 |
readerevt->AddVariable( "diphoPt0", &fMvaPEvtVars[0] );
|
| 65 |
readerevt->AddVariable( "nVert", &fMvaPEvtVars[1] );
|
| 66 |
readerevt->AddVariable( "MVA0", &fMvaPEvtVars[2] );
|
| 67 |
readerevt->AddVariable( "MVA1", &fMvaPEvtVars[3] );
|
| 68 |
readerevt->AddVariable( "dZ1", &fMvaPEvtVars[4] );
|
| 69 |
readerevt->AddVariable( "MVA2", &fMvaPEvtVars[5] );
|
| 70 |
readerevt->AddVariable( "dZ2", &fMvaPEvtVars[6] );
|
| 71 |
readerevt->AddVariable( "nConv", &fMvaPEvtVars[7] );
|
| 72 |
readerevt->BookMVA( "BDTEvt", gSystem->Getenv("CMSSW_BASE") + TString("/src/MitPhysics/data/TMVAClassification_evtBDTG_conversions_tmva_407.weights.xml") );
|
| 73 |
|
| 74 |
fIsInitMvaP = kTRUE;
|
| 75 |
}
|
| 76 |
|
| 77 |
|
| 78 |
double VertexTools::NewMass(const Photon* ph1, const Photon* ph2, const BaseVertex* vert)
|
| 79 |
{
|
| 80 |
ThreeVector drv1 = (ThreeVector(ph1->CaloPos()) - vert->Position()).Unit();
|
| 81 |
//ThreeVector drv1 = (ThreeVector(ph1->SCluster()->Point()) - vert->Position()).Unit();
|
| 82 |
FourVector pho1c(drv1.X()*ph1->E(),drv1.Y()*ph1->E(),drv1.Z()*ph1->E(),ph1->E());
|
| 83 |
ThreeVector drv2 = (ThreeVector(ph2->CaloPos()) - vert->Position()).Unit();
|
| 84 |
FourVector pho2c(drv2.X()*ph2->E(),drv2.Y()*ph2->E(),drv2.Z()*ph2->E(),ph2->E());
|
| 85 |
|
| 86 |
FourVector diboso = pho1c+pho2c;
|
| 87 |
return diboso.M();
|
| 88 |
}
|
| 89 |
|
| 90 |
//--------------------------------------------------------------------------------------------------
|
| 91 |
VertexZarray VertexTools::ExtractZarray(const VertexCol* vcol, float zmin, float zmax, const BaseVertex *fBeamSpot)
|
| 92 |
{
|
| 93 |
VertexZarray zs;
|
| 94 |
if(vcol == NULL) return zs;
|
| 95 |
|
| 96 |
for(unsigned vv = 0; vv < vcol->GetEntries(); vv++){
|
| 97 |
const Vertex* vert = vcol->At(vv);
|
| 98 |
double zpos = vert->Z();
|
| 99 |
if(fBeamSpot != NULL)
|
| 100 |
zpos = zpos - fBeamSpot->Z();
|
| 101 |
if(zpos > zmin && zpos > zmin)
|
| 102 |
zs.push_back(zpos);
|
| 103 |
}
|
| 104 |
return zs;
|
| 105 |
}
|
| 106 |
|
| 107 |
VertexZarray VertexTools::ExtractZarray(float zmin, float zmax, float step)
|
| 108 |
{
|
| 109 |
VertexZarray zs;
|
| 110 |
for(float zpos = zmin; zpos < zmax+step; zpos = zpos+step){
|
| 111 |
zs.push_back(zpos);
|
| 112 |
}
|
| 113 |
return zs;
|
| 114 |
}
|
| 115 |
|
| 116 |
const Vertex* VertexTools::BestVtx( const PFCandidateCol *fPFJets, const VertexCol *c,
|
| 117 |
const BaseVertex *fBeamSpot, FourVector diboso) {
|
| 118 |
|
| 119 |
if (!c || !c->GetEntries()) return NULL;
|
| 120 |
|
| 121 |
double bestprob = -100.;
|
| 122 |
const Vertex* bestvert = NULL;
|
| 123 |
for(unsigned vv = 0; vv < c->GetEntries(); vv++){
|
| 124 |
|
| 125 |
const Vertex* vert = c->At(vv);
|
| 126 |
double zpos = vert->Z();
|
| 127 |
double prob = Prob(fPFJets, zpos, fBeamSpot, diboso);
|
| 128 |
if(prob > bestprob){
|
| 129 |
bestprob = prob;
|
| 130 |
bestvert = vert;
|
| 131 |
}
|
| 132 |
}
|
| 133 |
return bestvert;
|
| 134 |
}
|
| 135 |
|
| 136 |
double VertexTools::BestVtx( const PFCandidateCol *fPFJets, VertexZarray zcol,
|
| 137 |
const BaseVertex *fBeamSpot, FourVector diboso){
|
| 138 |
|
| 139 |
double bestprob = -100.;
|
| 140 |
double bestz = -100.;
|
| 141 |
for(unsigned vv = 0; vv < zcol.size(); vv++){
|
| 142 |
double zpos = zcol[vv];
|
| 143 |
double prob = Prob(fPFJets, zpos, fBeamSpot, diboso);
|
| 144 |
if(prob > bestprob){
|
| 145 |
bestprob = prob;
|
| 146 |
bestz = zpos;
|
| 147 |
}
|
| 148 |
}
|
| 149 |
return bestz;
|
| 150 |
|
| 151 |
}
|
| 152 |
|
| 153 |
double VertexTools::Prob(const PFCandidateCol *fPFJets, double zpos,
|
| 154 |
const BaseVertex *fBeamSpot, FourVector diboso){
|
| 155 |
|
| 156 |
double bosophi = diboso.Phi();
|
| 157 |
double bosopt = diboso.Pt();
|
| 158 |
|
| 159 |
Vertex* vert = new Vertex(0,0,zpos);
|
| 160 |
double ZVC = vert->Z()-fBeamSpot->Z();
|
| 161 |
VertexTools* vtool = VertexTools::instance("");
|
| 162 |
|
| 163 |
Double_t sinsum = 0.;
|
| 164 |
Double_t cossum = 0.;
|
| 165 |
Double_t sumpt = 0.;
|
| 166 |
Double_t ntrks = 0.;
|
| 167 |
Double_t ntplus = 0.;
|
| 168 |
Double_t ortminus = 0.;
|
| 169 |
Double_t ortplus = 0.;
|
| 170 |
Double_t bdplus = 0.;
|
| 171 |
Double_t bdminus = 0.;
|
| 172 |
Double_t zmean = 0.;
|
| 173 |
Double_t zmeansq = 0.;
|
| 174 |
Double_t ww = 0.;
|
| 175 |
|
| 176 |
for(unsigned pfj = 0; pfj < fPFJets->GetEntries(); pfj++){
|
| 177 |
const PFCandidate* pfca = fPFJets->At(pfj);
|
| 178 |
if(! pfca->HasTrk()) continue;
|
| 179 |
if(!(pfca->PFType() == PFCandidate::eX ||
|
| 180 |
pfca->PFType() == PFCandidate::eHadron ||
|
| 181 |
pfca->PFType() == PFCandidate::eElectron ||
|
| 182 |
pfca->PFType() == PFCandidate::eMuon) ) continue;
|
| 183 |
const Track *t = pfca->Trk();
|
| 184 |
if(fabs(t->DzCorrected(*fBeamSpot) - ZVC ) > 0.2) continue;
|
| 185 |
|
| 186 |
if(pfca->Pt()<0.3 || pfca->Pt()>200) continue;
|
| 187 |
|
| 188 |
std::vector<const Track*>::iterator itt;
|
| 189 |
itt = find ((vtool->excluded).begin(), (vtool->excluded).end(), t);
|
| 190 |
if(itt != (vtool->excluded).end()) continue;
|
| 191 |
|
| 192 |
zmean = zmean + pfca->Pt()* t->DzCorrected(*fBeamSpot);
|
| 193 |
zmeansq = zmeansq + pfca->Pt()*t->DzCorrected(*fBeamSpot)*t->DzCorrected(*fBeamSpot);
|
| 194 |
ntrks++;
|
| 195 |
sumpt = sumpt + pfca->Pt()*pfca->Pt();
|
| 196 |
ww = ww + pfca->Pt();
|
| 197 |
|
| 198 |
sinsum = sinsum + pfca->Pt()*TMath::Sin(t->Phi());
|
| 199 |
cossum = cossum + pfca->Pt()*TMath::Cos(t->Phi());
|
| 200 |
}
|
| 201 |
if(ntrks < 2 || !(sumpt > 0.) ) return 0;
|
| 202 |
|
| 203 |
Double_t phim = TMath::ATan2(sinsum,cossum);
|
| 204 |
zmean = zmean/ww;
|
| 205 |
zmeansq = sqrt(zmeansq/ww - zmean*zmean);
|
| 206 |
|
| 207 |
//--------------------
|
| 208 |
Double_t bosoproj = 0.;
|
| 209 |
Float_t ymean = 0.;
|
| 210 |
Float_t ymsq = 0.;
|
| 211 |
for(unsigned pfj = 0; pfj < fPFJets->GetEntries(); pfj++){
|
| 212 |
const PFCandidate* pfca = fPFJets->At(pfj);
|
| 213 |
if(! pfca->HasTrk()) continue;
|
| 214 |
if(!(pfca->PFType() == PFCandidate::eX ||
|
| 215 |
pfca->PFType() == PFCandidate::eHadron ||
|
| 216 |
pfca->PFType() == PFCandidate::eElectron ||
|
| 217 |
pfca->PFType() == PFCandidate::eMuon) ) continue;
|
| 218 |
|
| 219 |
const Track *t = pfca->Trk();
|
| 220 |
if(fabs(t->DzCorrected(*fBeamSpot) - ZVC ) > 0.2) continue;
|
| 221 |
//if(fabs(t->DzCorrected(*fBeamSpot) - ZVC ) > 3*zwidth) continue;
|
| 222 |
if(pfca->Pt()<0.3 || pfca->Pt()>200) continue;
|
| 223 |
|
| 224 |
std::vector<const Track*>::iterator itt;
|
| 225 |
itt = find ((vtool->excluded).begin(), (vtool->excluded).end(), t);
|
| 226 |
if(itt != (vtool->excluded).end()) continue;
|
| 227 |
|
| 228 |
//Float_t phid = phim - t->Phi();
|
| 229 |
Float_t phid = bosophi+3.14 - t->Phi();
|
| 230 |
ymean = ymean + pfca->Pt()*TMath::Sin(phid);
|
| 231 |
ymsq = ymsq + pow(pfca->Pt()*TMath::Sin(phid),2);
|
| 232 |
bosoproj = bosoproj + pfca->Pt()*TMath::Cos(bosophi-t->Phi());
|
| 233 |
if(TMath::Sin(phid) > 0.) ortplus = ortplus + pfca->Pt()*TMath::Sin(phid);
|
| 234 |
else ortminus = ortminus + pfca->Pt()*fabs(TMath::Sin(phid));
|
| 235 |
if(TMath::Cos(phid) > 0.) bdplus = bdplus + pfca->Pt()*TMath::Cos(phid);
|
| 236 |
else bdminus = bdminus + pfca->Pt()*fabs(TMath::Cos(phid));
|
| 237 |
if(TMath::Cos(phid) > 0.) ntplus = ntplus + TMath::Cos(phid);
|
| 238 |
}
|
| 239 |
|
| 240 |
Float_t phimsq = sqrt(ymsq/ntrks - pow(ymean/ntrks,2))*180/3.14;
|
| 241 |
// not used... commented out by Fabian
|
| 242 |
//Double_t A2n = (ortplus+ortminus) > 0.01 ? (bdplus+bdminus)/(ortplus+ortminus) : (bdplus+bdminus)/0.01;
|
| 243 |
Double_t A1n = bdplus/bosopt;
|
| 244 |
Double_t A3n = ntplus > 0 ? bdplus/ntplus : 0.;
|
| 245 |
//-------------------
|
| 246 |
Double_t angle = 180/3.14*TVector2::Phi_0_2pi(bosophi-phim);
|
| 247 |
|
| 248 |
vtool->tmvar1 = ntrks;
|
| 249 |
vtool->tmvar2 = sumpt;
|
| 250 |
vtool->tmvar3 = A1n;
|
| 251 |
vtool->tmvar4 = angle;
|
| 252 |
vtool->tmvar5 = phimsq;
|
| 253 |
vtool->tmvar6 = A3n;
|
| 254 |
|
| 255 |
double p1 = vtool->reader->GetProba ( "BDTD method" );
|
| 256 |
double p2 = vtool->reader->GetProba ( "BDTG method" );
|
| 257 |
double p3 = vtool->reader->GetProba ( "MLP method");
|
| 258 |
double p4 = vtool->reader->GetProba ( "MLPBFGS method");
|
| 259 |
|
| 260 |
double retval = p1*p2*p3*p4;
|
| 261 |
|
| 262 |
return retval;
|
| 263 |
}
|
| 264 |
|
| 265 |
double VertexTools::VertexWidth(const Vertex* vert, const BaseVertex *fBeamSpot){
|
| 266 |
double width = 0.;
|
| 267 |
double zmean = 0.;
|
| 268 |
double zmeansq = 0.;
|
| 269 |
double ww = 0.;
|
| 270 |
if(vert == NULL) return width;
|
| 271 |
for(unsigned i = 0; i < vert->NTracks(); i++){
|
| 272 |
const Track *t = vert->Trk(i);
|
| 273 |
if(t->Pt() < 0.3 ) continue;
|
| 274 |
if(t->Pt() > 500.) continue;
|
| 275 |
ww = ww + t->Pt();
|
| 276 |
zmean = zmean + t->Pt()*t->DzCorrected(*fBeamSpot);
|
| 277 |
zmeansq = zmeansq + t->Pt()*t->DzCorrected(*fBeamSpot)*t->DzCorrected(*fBeamSpot);
|
| 278 |
}
|
| 279 |
if( !(ww > 0.) ) return width;
|
| 280 |
zmean = zmean/ww;
|
| 281 |
|
| 282 |
width = sqrt(zmeansq/ww - zmean*zmean);
|
| 283 |
return width;
|
| 284 |
|
| 285 |
}
|
| 286 |
|
| 287 |
void VertexTools::BanThisTrack(const Track* track){
|
| 288 |
VertexTools* vtool = VertexTools::instance("");
|
| 289 |
(vtool->excluded).push_back(track);
|
| 290 |
}
|
| 291 |
void VertexTools::Reset(){
|
| 292 |
VertexTools* vtool = VertexTools::instance("");
|
| 293 |
(vtool->excluded).clear();
|
| 294 |
}
|
| 295 |
|
| 296 |
|
| 297 |
//------------------------------------------------------------------------------------
|
| 298 |
// The below tools are from the H->2photons EPS2011 BaseLine (common) Selection
|
| 299 |
const Vertex* VertexTools::findVtxBasicRanking(const Photon* ph1,
|
| 300 |
const Photon* ph2,
|
| 301 |
const BaseVertex* bsp,
|
| 302 |
const VertexCol* vtcs,
|
| 303 |
const DecayParticleCol* conv, Bool_t useMva, Double_t &vtxProb) {
|
| 304 |
|
| 305 |
//if (useMva) printf("using mva vertex selection\n");
|
| 306 |
|
| 307 |
// check if all input is valid
|
| 308 |
if( !ph1 || !ph2 || !bsp || !vtcs ) return NULL;
|
| 309 |
// CAUTION: We allow for passing NULL for the Conversions, in that case only the simple Ranking is used.
|
| 310 |
|
| 311 |
// here we will store the idx of the best Vtx
|
| 312 |
unsigned int bestIdx = 0;
|
| 313 |
UInt_t bestidxmva = 0;
|
| 314 |
|
| 315 |
// using asd much as possible 'Globe' naming schemes...
|
| 316 |
int* ptbal_rank = new int [vtcs->GetEntries()];
|
| 317 |
int* ptasym_rank = new int [vtcs->GetEntries()];
|
| 318 |
int* total_rank = new int [vtcs->GetEntries()];
|
| 319 |
double* ptbal = new double[vtcs->GetEntries()];
|
| 320 |
double* ptasym = new double[vtcs->GetEntries()];
|
| 321 |
double* sumpt2 = new double[vtcs->GetEntries()];
|
| 322 |
double* limPullToConv = new double[vtcs->GetEntries()];
|
| 323 |
double* mvaval = new double[vtcs->GetEntries()];
|
| 324 |
|
| 325 |
|
| 326 |
unsigned int numVertices = vtcs->GetEntries();
|
| 327 |
double ptgg = 0.; // stored for later in the conversion
|
| 328 |
|
| 329 |
// loop over all the vertices...
|
| 330 |
for(unsigned int iVtx = 0; iVtx < numVertices; ++iVtx) {
|
| 331 |
|
| 332 |
const Vertex* tVtx = vtcs->At(iVtx);
|
| 333 |
ptbal [iVtx] = 0.0;
|
| 334 |
ptasym [iVtx] = 0.0;
|
| 335 |
sumpt2 [iVtx] = 0.0;
|
| 336 |
limPullToConv [iVtx] = -1.0;
|
| 337 |
ptbal_rank [iVtx] = 1;
|
| 338 |
ptasym_rank[iVtx] = 1;
|
| 339 |
|
| 340 |
// compute the photon momenta with respect to this Vtx
|
| 341 |
FourVectorM newMomFst = ph1->MomVtx(tVtx->Position());
|
| 342 |
FourVectorM newMomSec = ph2->MomVtx(tVtx->Position());
|
| 343 |
|
| 344 |
FourVectorM higgsMom = newMomFst+newMomSec;
|
| 345 |
|
| 346 |
double ph1Eta = newMomFst.Eta();
|
| 347 |
double ph2Eta = newMomSec.Eta();
|
| 348 |
|
| 349 |
double ph1Phi = newMomFst.Phi();
|
| 350 |
double ph2Phi = newMomSec.Phi();
|
| 351 |
|
| 352 |
// loop over all tracks and computew variables for ranking...
|
| 353 |
FourVectorM totTrkMom(0,0,0,0);
|
| 354 |
for(unsigned int iTrk = 0; iTrk < tVtx->NTracks(); ++iTrk) {
|
| 355 |
const Track* tTrk = tVtx->Trk(iTrk);
|
| 356 |
|
| 357 |
sumpt2[iVtx] += tTrk->Pt()*tTrk->Pt();
|
| 358 |
|
| 359 |
// compute distance between Trk and the Photons
|
| 360 |
double tEta = tTrk->Eta();
|
| 361 |
double tPhi = tTrk->Phi();
|
| 362 |
double dEta1 = TMath::Abs(tEta-ph1Eta);
|
| 363 |
double dEta2 = TMath::Abs(tEta-ph2Eta);
|
| 364 |
double dPhi1 = TMath::Abs(tPhi-ph1Phi);
|
| 365 |
double dPhi2 = TMath::Abs(tPhi-ph2Phi);
|
| 366 |
if(dPhi1 > M_PI) dPhi1 = 2*M_PI - dPhi1;
|
| 367 |
if(dPhi2 > M_PI) dPhi2 = 2*M_PI - dPhi2;
|
| 368 |
|
| 369 |
double dR1 = TMath::Sqrt(dEta1*dEta1+dPhi1*dPhi1);
|
| 370 |
double dR2 = TMath::Sqrt(dEta2*dEta2+dPhi2*dPhi2);
|
| 371 |
|
| 372 |
if(dR1 < 0.05 || dR2 < 0.05) continue;
|
| 373 |
totTrkMom = totTrkMom + tTrk->Mom4(0);
|
| 374 |
}
|
| 375 |
|
| 376 |
// compute the ranking variables for this Vtx
|
| 377 |
double ptvtx = totTrkMom.Pt();
|
| 378 |
double pthiggs = higgsMom.Pt();
|
| 379 |
if(iVtx == 0) ptgg = pthiggs;
|
| 380 |
ptbal [iVtx] = (totTrkMom.Px()*(newMomFst.Px()+newMomSec.Px()));
|
| 381 |
ptbal [iVtx] += (totTrkMom.Py()*(newMomFst.Py()+newMomSec.Py()));
|
| 382 |
ptbal [iVtx] = -ptbal[iVtx]/pthiggs;
|
| 383 |
ptasym[iVtx] = (ptvtx - pthiggs)/(ptvtx + pthiggs);
|
| 384 |
|
| 385 |
// do the little ranking acrobatics...
|
| 386 |
for(unsigned int cVtx =0; cVtx < iVtx; ++cVtx) {
|
| 387 |
if(ptbal [iVtx] > ptbal [cVtx])
|
| 388 |
ptbal_rank[cVtx]++;
|
| 389 |
else
|
| 390 |
ptbal_rank[iVtx]++;
|
| 391 |
if(ptasym [iVtx] > ptasym [cVtx])
|
| 392 |
ptasym_rank[cVtx]++;
|
| 393 |
else
|
| 394 |
ptasym_rank[iVtx]++;
|
| 395 |
}
|
| 396 |
}
|
| 397 |
|
| 398 |
// loop again over all Vertcices (*sigh*), compute total score and final rank
|
| 399 |
// CAUTION: Total rank starts at '0', so the best ranked Vtx has RANK 0
|
| 400 |
for(unsigned int iVtx = 0; iVtx < numVertices; ++iVtx) {
|
| 401 |
ptasym_rank [iVtx] = ptbal_rank [iVtx]*ptasym_rank [iVtx]*(iVtx+1);
|
| 402 |
total_rank [iVtx] = 0;
|
| 403 |
for(unsigned int cVtx =0; cVtx < iVtx; ++cVtx) {
|
| 404 |
if(ptasym_rank [iVtx] > ptasym_rank [cVtx]) total_rank[iVtx]++;
|
| 405 |
else if(ptasym_rank [iVtx] == ptasym_rank [cVtx]) { // use 'ptbal' as the tie-breaker
|
| 406 |
if(ptbal_rank [iVtx] > ptbal_rank [cVtx]) total_rank[iVtx]++;
|
| 407 |
else total_rank[cVtx]++;
|
| 408 |
}
|
| 409 |
else total_rank[cVtx]++;
|
| 410 |
}
|
| 411 |
}
|
| 412 |
|
| 413 |
|
| 414 |
// find the best ranked Vertex so far....
|
| 415 |
for(unsigned int iVtx = 0; iVtx < numVertices; ++iVtx) {
|
| 416 |
if(total_rank[iVtx] == 0) bestIdx = iVtx;
|
| 417 |
}
|
| 418 |
|
| 419 |
// check if there's a conversion among the pre-selected photons
|
| 420 |
// ...this will return NULL in case conv==NULL
|
| 421 |
const DecayParticle* conv1 = PhotonTools::MatchedCiCConversion(ph1, conv, 0.1, 0.1, 0.1);
|
| 422 |
const DecayParticle* conv2 = PhotonTools::MatchedCiCConversion(ph2, conv, 0.1, 0.1, 0.1);
|
| 423 |
|
| 424 |
double zconv = 0.;
|
| 425 |
double dzconv = 0.;
|
| 426 |
int nConv = 0;
|
| 427 |
|
| 428 |
if (conv1) nConv += 1;
|
| 429 |
if (conv2) nConv += 1;
|
| 430 |
|
| 431 |
|
| 432 |
const double dzpxb = 0.016;
|
| 433 |
const double dztib = 0.331;
|
| 434 |
const double dztob = 1.564;
|
| 435 |
const double dzpxf = 0.082;
|
| 436 |
const double dztid = 0.321;
|
| 437 |
const double dztec = 0.815;
|
| 438 |
|
| 439 |
//--------------------------------------------------------------------
|
| 440 |
// start doing the Conversion acrobatics... 'copied' from the Globe...
|
| 441 |
if(conv1 || conv2) {
|
| 442 |
if( !conv2 ){
|
| 443 |
const mithep::ThreeVector caloPos1(ph1->CaloPos());
|
| 444 |
double zconvsc = conv1->Z0EcalVtxCiC(bsp->Position(), caloPos1);
|
| 445 |
double zconvtrk = conv1->DzCorrected(bsp->Position()) + bsp->Z();
|
| 446 |
if( ph1->IsEB() ) {
|
| 447 |
double rho = conv1->Position().Rho();
|
| 448 |
if ( rho < 15. ) { dzconv = dzpxb; zconv = zconvtrk; }
|
| 449 |
else if( rho < 60. ) { dzconv = dztib; zconv = zconvsc; }
|
| 450 |
else { dzconv = dztob; zconv = zconvsc; }
|
| 451 |
} else {
|
| 452 |
double z = conv1->Position().Z();
|
| 453 |
if ( TMath::Abs(z) < 50. ) { dzconv = dzpxf; zconv = zconvtrk; }
|
| 454 |
else if( TMath::Abs(z) < 100.) { dzconv = dztid; zconv = zconvtrk; }
|
| 455 |
else { dzconv = dztec; zconv = zconvsc; }
|
| 456 |
}
|
| 457 |
} else if( !conv1 ) {
|
| 458 |
const mithep::ThreeVector caloPos2(ph2->CaloPos());
|
| 459 |
double zconvsc = conv2->Z0EcalVtxCiC(bsp->Position(), caloPos2);
|
| 460 |
double zconvtrk = conv2->DzCorrected(bsp->Position()) + bsp->Z();
|
| 461 |
if( ph2->IsEB() ) {
|
| 462 |
double rho = conv2->Position().Rho();
|
| 463 |
if ( rho < 15. ) { dzconv = dzpxb; zconv = zconvtrk; }
|
| 464 |
else if( rho < 60. ) { dzconv = dztib; zconv = zconvsc; }
|
| 465 |
else { dzconv = dztob; zconv = zconvsc; }
|
| 466 |
} else {
|
| 467 |
double z = conv2->Position().Z();
|
| 468 |
if ( TMath::Abs(z) < 50. ) { dzconv = dzpxf; zconv = zconvtrk; }
|
| 469 |
else if( TMath::Abs(z) < 100.) { dzconv = dztid; zconv = zconvtrk; }
|
| 470 |
else { dzconv = dztec; zconv = zconvsc; }
|
| 471 |
}
|
| 472 |
} else {
|
| 473 |
const mithep::ThreeVector caloPos1(ph1->CaloPos());
|
| 474 |
double z1=0.;
|
| 475 |
double z1sc = conv1->Z0EcalVtxCiC(bsp->Position(), caloPos1);
|
| 476 |
double z1trk = conv1->DzCorrected(bsp->Position()) + bsp->Z();
|
| 477 |
double dz1 = 0.;
|
| 478 |
if( ph1->IsEB() ) {
|
| 479 |
double rho = conv1->Position().Rho();
|
| 480 |
if ( rho < 15. ) { dz1 = dzpxb; z1 = z1trk; }
|
| 481 |
else if( rho < 60. ) { dz1 = dztib; z1 = z1sc; }
|
| 482 |
else { dz1 = dztob; z1 = z1sc; }
|
| 483 |
} else {
|
| 484 |
double z = conv1->Position().Z();
|
| 485 |
if ( TMath::Abs(z) < 50. ) { dz1 = dzpxf; z1 = z1trk; }
|
| 486 |
else if( TMath::Abs(z) < 100.) { dz1 = dztid; z1 = z1trk; }
|
| 487 |
else { dz1 = dztec; z1 = z1sc; }
|
| 488 |
}
|
| 489 |
const mithep::ThreeVector caloPos2(ph2->CaloPos());
|
| 490 |
double z2 = 0.;
|
| 491 |
double z2sc = conv2->Z0EcalVtxCiC(bsp->Position(), caloPos2);
|
| 492 |
double z2trk = conv2->DzCorrected(bsp->Position()) + bsp->Z();
|
| 493 |
double dz2 = 0.;
|
| 494 |
if( ph2->IsEB() ) {
|
| 495 |
double rho = conv2->Position().Rho();
|
| 496 |
if ( rho < 15. ) { dz2 = dzpxb; z2 = z2trk; }
|
| 497 |
else if( rho < 60. ) { dz2 = dztib; z2 = z2sc; }
|
| 498 |
else { dz2 = dztob; z2 = z2sc; }
|
| 499 |
} else {
|
| 500 |
double z = conv2->Position().Z();
|
| 501 |
if ( TMath::Abs(z) < 50. ) { dz2 = dzpxf; z2 = z1trk; }
|
| 502 |
else if( TMath::Abs(z) < 100.) { dz2 = dztid; z2 = z1trk; }
|
| 503 |
else { dz2 = dztec; z2 = z1sc; }
|
| 504 |
}
|
| 505 |
|
| 506 |
zconv = ( 1./(1./dz1/dz1 + 1./dz2/dz2 )*(z1/dz1/dz1 + z2/dz2/dz2) ) ; // weighted average
|
| 507 |
dzconv = TMath::Sqrt( 1./(1./dz1/dz1 + 1./dz2/dz2)) ;
|
| 508 |
}
|
| 509 |
|
| 510 |
|
| 511 |
// loop over all ranked Vertices and choose the closest to the Conversion one
|
| 512 |
int maxVertices = ( ptgg > 30 ? 3 : 5);
|
| 513 |
double minDz = -1.;
|
| 514 |
|
| 515 |
|
| 516 |
|
| 517 |
for(unsigned int iVtx =0; iVtx < numVertices; ++iVtx) {
|
| 518 |
|
| 519 |
limPullToConv[iVtx] = TMath::Abs(vtcs->At(iVtx)->Z()-zconv)/dzconv;
|
| 520 |
|
| 521 |
if(total_rank[iVtx] < maxVertices) {
|
| 522 |
const Vertex* tVtx = vtcs->At(iVtx);
|
| 523 |
double tDz = TMath::Abs(zconv - tVtx->Z());
|
| 524 |
|
| 525 |
if( (minDz < 0. || tDz < minDz) && ( tDz < dzconv ) ) {
|
| 526 |
minDz = tDz;
|
| 527 |
bestIdx = iVtx;
|
| 528 |
}
|
| 529 |
}
|
| 530 |
}
|
| 531 |
}
|
| 532 |
// END of Conversion Acrobatics
|
| 533 |
//--------------------------------------------------------------------
|
| 534 |
|
| 535 |
//final loop to compute mva values
|
| 536 |
double mvamax = -1e6;
|
| 537 |
for(unsigned int iVtx =0; iVtx < numVertices; ++iVtx) {
|
| 538 |
double mva = VtxMvaP(ptbal[iVtx],ptasym[iVtx],log(sumpt2[iVtx]),limPullToConv[iVtx],nConv);
|
| 539 |
mvaval[iVtx] = mva;
|
| 540 |
if (mva>mvamax) {
|
| 541 |
mvamax = mva;
|
| 542 |
bestidxmva = iVtx;
|
| 543 |
}
|
| 544 |
//printf("vtx %i: ptbal = %5f, ptasym = %5f, logsumpt2 = %5f, limpulltoconv = %5f, nconv = %i, mva = %5f\n",iVtx,ptbal[iVtx],ptasym[iVtx],log(sumpt2[iVtx]),limPullToConv[iVtx],nConv,mva);
|
| 545 |
}
|
| 546 |
|
| 547 |
// double mvahack = VtxMvaP(4.13519,-0.156296,3.17947,-1.0,0);
|
| 548 |
// printf("mvahack = %5f\n",mvahack);
|
| 549 |
|
| 550 |
//find second and third ranked vertices for event mva;
|
| 551 |
UInt_t mvaidx1 = 0;
|
| 552 |
mvamax = -1e6;
|
| 553 |
for(unsigned int iVtx =0; iVtx < numVertices; ++iVtx) {
|
| 554 |
if (iVtx!=bestidxmva && mvaval[iVtx]>mvamax) {
|
| 555 |
mvamax = mvaval[iVtx];
|
| 556 |
mvaidx1 = iVtx;
|
| 557 |
}
|
| 558 |
}
|
| 559 |
|
| 560 |
UInt_t mvaidx2 = 0;
|
| 561 |
mvamax = -1e6;
|
| 562 |
for(unsigned int iVtx =0; iVtx < numVertices; ++iVtx) {
|
| 563 |
if (iVtx!=bestidxmva && iVtx!=mvaidx1 && mvaval[iVtx]>mvamax) {
|
| 564 |
mvamax = mvaval[iVtx];
|
| 565 |
mvaidx2 = iVtx;
|
| 566 |
}
|
| 567 |
}
|
| 568 |
|
| 569 |
//compute per event mva output
|
| 570 |
FourVectorM newMomFst = ph1->MomVtx(vtcs->At(bestidxmva)->Position());
|
| 571 |
FourVectorM newMomSec = ph2->MomVtx(vtcs->At(bestidxmva)->Position());
|
| 572 |
FourVectorM higgsMom = newMomFst+newMomSec;
|
| 573 |
|
| 574 |
fMvaPEvtVars[0] = higgsMom.Pt();
|
| 575 |
fMvaPEvtVars[1] = numVertices;
|
| 576 |
fMvaPEvtVars[2] = mvaval[bestidxmva];
|
| 577 |
fMvaPEvtVars[3] = mvaval[mvaidx1];
|
| 578 |
fMvaPEvtVars[4] = vtcs->At(mvaidx1)->Z() - vtcs->At(bestidxmva)->Z();
|
| 579 |
fMvaPEvtVars[5] = mvaval[mvaidx2];
|
| 580 |
fMvaPEvtVars[6] = vtcs->At(mvaidx2)->Z() - vtcs->At(bestidxmva)->Z();
|
| 581 |
fMvaPEvtVars[7] = nConv;
|
| 582 |
|
| 583 |
Double_t evtmva = readerevt->EvaluateMVA("BDTEvt");
|
| 584 |
vtxProb = 1.-0.49*(evtmva+1.0);
|
| 585 |
|
| 586 |
// printf("higgspt = %5f, numvert = %5f, mvabest = %5f, mva1 = %5f, dz1 = %5f, mva2 = %5f, dz2 = %5f, nconv = %5f",fMvaPEvtVars[0],fMvaPEvtVars[1],fMvaPEvtVars[2],fMvaPEvtVars[3],fMvaPEvtVars[4],fMvaPEvtVars[5],fMvaPEvtVars[6],fMvaPEvtVars[7]);
|
| 587 |
// printf("vtxmva = %5f, vtxprob = %5f\n",evtmva,vtxProb);
|
| 588 |
//
|
| 589 |
// printf("e1 = %5f, sige1 = %5f\n",ph1->E(),ph1->EnergyErr());
|
| 590 |
// printf("e2 = %5f, sige2 = %5f\n",ph2->E(),ph2->EnergyErr());
|
| 591 |
|
| 592 |
// delete the auxiliary dynamic arrays
|
| 593 |
delete[] ptbal_rank ;
|
| 594 |
delete[] ptasym_rank ;
|
| 595 |
delete[] ptbal ;
|
| 596 |
delete[] ptasym ;
|
| 597 |
delete[] sumpt2 ;
|
| 598 |
delete[] limPullToConv;
|
| 599 |
delete[] mvaval;
|
| 600 |
|
| 601 |
delete[] total_rank ;
|
| 602 |
|
| 603 |
|
| 604 |
if (useMva) return vtcs->At(bestidxmva);
|
| 605 |
else return vtcs->At(bestIdx);
|
| 606 |
}
|
| 607 |
|
| 608 |
//------------------------------------------------------------------------------------
|
| 609 |
double VertexTools::VtxMvaP(float ptbal, float ptasym, float logsumpt2, float limPullToConv, float nConv) const
|
| 610 |
{
|
| 611 |
fMvaPVars[0] = ptbal;
|
| 612 |
fMvaPVars[1] = ptasym;
|
| 613 |
fMvaPVars[2] = logsumpt2;
|
| 614 |
fMvaPVars[3] = limPullToConv;
|
| 615 |
fMvaPVars[4] = nConv;
|
| 616 |
|
| 617 |
return readervtx->EvaluateMVA("BDTCat");
|
| 618 |
|
| 619 |
}
|
| 620 |
|
| 621 |
//------------------------------------------------------------------------------------
|
| 622 |
//Compute contribution to relative uncertainty sigma_m/m from primary vertex location
|
| 623 |
//given ecal shower positions of two photons, plus the vtx z uncertainty (typically sqrt(2)*beamspot width)
|
| 624 |
//code originally from Y. Gershtein
|
| 625 |
Double_t VertexTools::DeltaMassVtx(Double_t xp1, Double_t yp1, Double_t zp1,
|
| 626 |
Double_t xp2, Double_t yp2, Double_t zp2,
|
| 627 |
Double_t xv, Double_t yv, Double_t zv,
|
| 628 |
Double_t dz)
|
| 629 |
{
|
| 630 |
|
| 631 |
Double_t x1 = xp1 - xv;
|
| 632 |
Double_t y1 = yp1 - yv;
|
| 633 |
Double_t z1 = zp1 - zv;
|
| 634 |
|
| 635 |
Double_t x2 = xp2 - xv;
|
| 636 |
Double_t y2 = yp2 - yv;
|
| 637 |
Double_t z2 = zp2 - zv;
|
| 638 |
|
| 639 |
Double_t r1 = sqrt(x1*x1+y1*y1+z1*z1);
|
| 640 |
Double_t r2 = sqrt(x2*x2+y2*y2+z2*z2);
|
| 641 |
Double_t phi1 = atan2(y1,x1);
|
| 642 |
Double_t theta1 = atan2(sqrt(x1*x1+y1*y1),z1);
|
| 643 |
Double_t phi2 = atan2(y2,x2);
|
| 644 |
Double_t theta2 = atan2(sqrt(x2*x2+y2*y2),z2);
|
| 645 |
|
| 646 |
Double_t sech1 = sin(theta1);
|
| 647 |
Double_t tanh1 = cos(theta1);
|
| 648 |
Double_t sech2 = sin(theta2);
|
| 649 |
Double_t tanh2 = cos(theta2);
|
| 650 |
Double_t cos12 = cos(phi1-phi2);
|
| 651 |
|
| 652 |
Double_t rad1 = sech1*(sech1*tanh2-tanh1*sech2*cos12)/(1-tanh1*tanh2-sech1*sech2*cos12);
|
| 653 |
Double_t rad2 = sech2*(sech2*tanh1-tanh2*sech1*cos12)/(1-tanh2*tanh1-sech2*sech1*cos12);
|
| 654 |
|
| 655 |
return dz * 0.5*fabs(rad1/r1 + rad2/r2);
|
| 656 |
|
| 657 |
} |
