| 1 |
fabstoec |
1.5 |
// $Id: VertexTools.cc,v 1.4 2011/07/04 13:48:50 fabstoec Exp $
|
| 2 |
bendavid |
1.1 |
|
| 3 |
|
|
#include "MitPhysics/Utils/interface/VertexTools.h"
|
| 4 |
|
|
#include "MitPhysics/Utils/interface/ElectronTools.h"
|
| 5 |
fabstoec |
1.5 |
#include "MitPhysics/Utils/interface/PhotonTools.h"
|
| 6 |
bendavid |
1.1 |
#include "MitAna/DataTree/interface/StableData.h"
|
| 7 |
|
|
#include <TFile.h>
|
| 8 |
maxi |
1.2 |
#include <TVector3.h>
|
| 9 |
bendavid |
1.1 |
|
| 10 |
|
|
ClassImp(mithep::VertexTools)
|
| 11 |
|
|
|
| 12 |
|
|
using namespace mithep;
|
| 13 |
|
|
|
| 14 |
maxi |
1.2 |
VertexTools* VertexTools::meobject = NULL;
|
| 15 |
|
|
|
| 16 |
bendavid |
1.1 |
//--------------------------------------------------------------------------------------------------
|
| 17 |
maxi |
1.2 |
VertexTools::VertexTools(const char* str)
|
| 18 |
bendavid |
1.1 |
{
|
| 19 |
|
|
// Constructor.
|
| 20 |
maxi |
1.2 |
relname = str;
|
| 21 |
|
|
reader = new TMVA::Reader( "!Color:!Silent" );
|
| 22 |
|
|
reader->AddVariable( "var1", &tmvar1 );
|
| 23 |
|
|
reader->AddVariable( "var2", &tmvar2 );
|
| 24 |
|
|
reader->AddVariable( "var3", &tmvar3 );
|
| 25 |
|
|
reader->AddVariable( "var4", &tmvar4 );
|
| 26 |
|
|
reader->AddVariable( "var5", &tmvar5 );
|
| 27 |
|
|
reader->AddVariable( "var6", &tmvar6 );
|
| 28 |
|
|
reader->BookMVA( "BDTG method",relname + TString("/src/MitPhysics/data/TMVAClassification_BDTG.weights.xml").Data());
|
| 29 |
|
|
reader->BookMVA( "BDTD method",relname + TString("/src/MitPhysics/data/TMVAClassification_BDTD.weights.xml" ).Data());
|
| 30 |
|
|
//reader->BookMVA( "CFMlpANN method", "/home/maxi/cms/root/TMVAClassification_CFMlpANN.weights.xml" );
|
| 31 |
|
|
reader->BookMVA( "MLP method", relname + TString("/src/MitPhysics/data/TMVAClassification_MLP.weights.xml").Data());
|
| 32 |
|
|
reader->BookMVA( "MLPBFGS method",relname + TString("/src/MitPhysics/data/TMVAClassification_MLPBFGS.weights.xml" ).Data());
|
| 33 |
|
|
}
|
| 34 |
|
|
|
| 35 |
|
|
double VertexTools::NewMass(const Photon* ph1, const Photon* ph2, const BaseVertex* vert)
|
| 36 |
|
|
{
|
| 37 |
maxi |
1.3 |
ThreeVector drv1 = (ThreeVector(ph1->CaloPos()) - vert->Position()).Unit();
|
| 38 |
|
|
//ThreeVector drv1 = (ThreeVector(ph1->SCluster()->Point()) - vert->Position()).Unit();
|
| 39 |
maxi |
1.2 |
FourVector pho1c(drv1.X()*ph1->E(),drv1.Y()*ph1->E(),drv1.Z()*ph1->E(),ph1->E());
|
| 40 |
maxi |
1.3 |
ThreeVector drv2 = (ThreeVector(ph2->CaloPos()) - vert->Position()).Unit();
|
| 41 |
maxi |
1.2 |
FourVector pho2c(drv2.X()*ph2->E(),drv2.Y()*ph2->E(),drv2.Z()*ph2->E(),ph2->E());
|
| 42 |
|
|
|
| 43 |
|
|
FourVector diboso = pho1c+pho2c;
|
| 44 |
|
|
return diboso.M();
|
| 45 |
bendavid |
1.1 |
}
|
| 46 |
|
|
|
| 47 |
|
|
//--------------------------------------------------------------------------------------------------
|
| 48 |
maxi |
1.2 |
VertexZarray VertexTools::ExtractZarray(const VertexCol* vcol, float zmin, float zmax, const BaseVertex *fBeamSpot)
|
| 49 |
|
|
{
|
| 50 |
|
|
VertexZarray zs;
|
| 51 |
|
|
if(vcol == NULL) return zs;
|
| 52 |
|
|
|
| 53 |
|
|
for(unsigned vv = 0; vv < vcol->GetEntries(); vv++){
|
| 54 |
|
|
const Vertex* vert = vcol->At(vv);
|
| 55 |
|
|
double zpos = vert->Z();
|
| 56 |
|
|
if(fBeamSpot != NULL)
|
| 57 |
|
|
zpos = zpos - fBeamSpot->Z();
|
| 58 |
|
|
if(zpos > zmin && zpos > zmin)
|
| 59 |
|
|
zs.push_back(zpos);
|
| 60 |
|
|
}
|
| 61 |
|
|
return zs;
|
| 62 |
|
|
}
|
| 63 |
|
|
|
| 64 |
|
|
VertexZarray VertexTools::ExtractZarray(float zmin, float zmax, float step)
|
| 65 |
|
|
{
|
| 66 |
|
|
VertexZarray zs;
|
| 67 |
|
|
for(float zpos = zmin; zpos < zmax+step; zpos = zpos+step){
|
| 68 |
|
|
zs.push_back(zpos);
|
| 69 |
|
|
}
|
| 70 |
|
|
return zs;
|
| 71 |
|
|
}
|
| 72 |
|
|
|
| 73 |
|
|
const Vertex* VertexTools::BestVtx( const PFCandidateCol *fPFJets, const VertexCol *c,
|
| 74 |
|
|
const BaseVertex *fBeamSpot, FourVector diboso) {
|
| 75 |
|
|
|
| 76 |
|
|
if (!c || !c->GetEntries()) return NULL;
|
| 77 |
|
|
|
| 78 |
|
|
double bestprob = -100.;
|
| 79 |
|
|
const Vertex* bestvert = NULL;
|
| 80 |
|
|
for(unsigned vv = 0; vv < c->GetEntries(); vv++){
|
| 81 |
|
|
|
| 82 |
|
|
const Vertex* vert = c->At(vv);
|
| 83 |
|
|
double zpos = vert->Z();
|
| 84 |
|
|
double prob = Prob(fPFJets, zpos, fBeamSpot, diboso);
|
| 85 |
|
|
if(prob > bestprob){
|
| 86 |
|
|
bestprob = prob;
|
| 87 |
|
|
bestvert = vert;
|
| 88 |
|
|
}
|
| 89 |
|
|
}
|
| 90 |
|
|
return bestvert;
|
| 91 |
|
|
}
|
| 92 |
|
|
|
| 93 |
|
|
double VertexTools::BestVtx( const PFCandidateCol *fPFJets, VertexZarray zcol,
|
| 94 |
|
|
const BaseVertex *fBeamSpot, FourVector diboso){
|
| 95 |
|
|
|
| 96 |
|
|
double bestprob = -100.;
|
| 97 |
|
|
double bestz = -100.;
|
| 98 |
|
|
for(unsigned vv = 0; vv < zcol.size(); vv++){
|
| 99 |
|
|
double zpos = zcol[vv];
|
| 100 |
|
|
double prob = Prob(fPFJets, zpos, fBeamSpot, diboso);
|
| 101 |
|
|
if(prob > bestprob){
|
| 102 |
|
|
bestprob = prob;
|
| 103 |
|
|
bestz = zpos;
|
| 104 |
|
|
}
|
| 105 |
|
|
}
|
| 106 |
|
|
return bestz;
|
| 107 |
|
|
|
| 108 |
|
|
}
|
| 109 |
bendavid |
1.1 |
|
| 110 |
maxi |
1.2 |
double VertexTools::Prob(const PFCandidateCol *fPFJets, double zpos,
|
| 111 |
|
|
const BaseVertex *fBeamSpot, FourVector diboso){
|
| 112 |
|
|
|
| 113 |
|
|
double bosophi = diboso.Phi();
|
| 114 |
|
|
double bosopt = diboso.Pt();
|
| 115 |
|
|
|
| 116 |
|
|
Vertex* vert = new Vertex(0,0,zpos);
|
| 117 |
|
|
double ZVC = vert->Z()-fBeamSpot->Z();
|
| 118 |
maxi |
1.3 |
VertexTools* vtool = VertexTools::instance("");
|
| 119 |
maxi |
1.2 |
|
| 120 |
|
|
Double_t sinsum = 0.;
|
| 121 |
|
|
Double_t cossum = 0.;
|
| 122 |
|
|
Double_t sumpt = 0.;
|
| 123 |
|
|
Double_t ntrks = 0.;
|
| 124 |
|
|
Double_t ntplus = 0.;
|
| 125 |
|
|
Double_t ortminus = 0.;
|
| 126 |
|
|
Double_t ortplus = 0.;
|
| 127 |
|
|
Double_t bdplus = 0.;
|
| 128 |
|
|
Double_t bdminus = 0.;
|
| 129 |
|
|
Double_t zmean = 0.;
|
| 130 |
|
|
Double_t zmeansq = 0.;
|
| 131 |
|
|
Double_t ww = 0.;
|
| 132 |
|
|
|
| 133 |
|
|
for(unsigned pfj = 0; pfj < fPFJets->GetEntries(); pfj++){
|
| 134 |
|
|
const PFCandidate* pfca = fPFJets->At(pfj);
|
| 135 |
|
|
if(! pfca->HasTrk()) continue;
|
| 136 |
|
|
if(!(pfca->PFType() == PFCandidate::eX ||
|
| 137 |
|
|
pfca->PFType() == PFCandidate::eHadron ||
|
| 138 |
|
|
pfca->PFType() == PFCandidate::eElectron ||
|
| 139 |
|
|
pfca->PFType() == PFCandidate::eMuon) ) continue;
|
| 140 |
|
|
const Track *t = pfca->Trk();
|
| 141 |
|
|
if(fabs(t->DzCorrected(*fBeamSpot) - ZVC ) > 0.2) continue;
|
| 142 |
|
|
|
| 143 |
|
|
if(pfca->Pt()<0.3 || pfca->Pt()>200) continue;
|
| 144 |
maxi |
1.3 |
|
| 145 |
|
|
std::vector<const Track*>::iterator itt;
|
| 146 |
|
|
itt = find ((vtool->excluded).begin(), (vtool->excluded).end(), t);
|
| 147 |
|
|
if(itt != (vtool->excluded).end()) continue;
|
| 148 |
maxi |
1.2 |
|
| 149 |
|
|
zmean = zmean + pfca->Pt()* t->DzCorrected(*fBeamSpot);
|
| 150 |
|
|
zmeansq = zmeansq + pfca->Pt()*t->DzCorrected(*fBeamSpot)*t->DzCorrected(*fBeamSpot);
|
| 151 |
|
|
ntrks++;
|
| 152 |
|
|
sumpt = sumpt + pfca->Pt()*pfca->Pt();
|
| 153 |
|
|
ww = ww + pfca->Pt();
|
| 154 |
|
|
|
| 155 |
|
|
sinsum = sinsum + pfca->Pt()*TMath::Sin(t->Phi());
|
| 156 |
|
|
cossum = cossum + pfca->Pt()*TMath::Cos(t->Phi());
|
| 157 |
|
|
}
|
| 158 |
|
|
if(ntrks < 2 || !(sumpt > 0.) ) return 0;
|
| 159 |
|
|
|
| 160 |
|
|
Double_t phim = TMath::ATan2(sinsum,cossum);
|
| 161 |
|
|
zmean = zmean/ww;
|
| 162 |
|
|
zmeansq = sqrt(zmeansq/ww - zmean*zmean);
|
| 163 |
|
|
|
| 164 |
|
|
//--------------------
|
| 165 |
|
|
Double_t bosoproj = 0.;
|
| 166 |
|
|
Float_t ymean = 0.;
|
| 167 |
|
|
Float_t ymsq = 0.;
|
| 168 |
|
|
for(unsigned pfj = 0; pfj < fPFJets->GetEntries(); pfj++){
|
| 169 |
|
|
const PFCandidate* pfca = fPFJets->At(pfj);
|
| 170 |
|
|
if(! pfca->HasTrk()) continue;
|
| 171 |
|
|
if(!(pfca->PFType() == PFCandidate::eX ||
|
| 172 |
|
|
pfca->PFType() == PFCandidate::eHadron ||
|
| 173 |
|
|
pfca->PFType() == PFCandidate::eElectron ||
|
| 174 |
|
|
pfca->PFType() == PFCandidate::eMuon) ) continue;
|
| 175 |
|
|
|
| 176 |
|
|
const Track *t = pfca->Trk();
|
| 177 |
|
|
if(fabs(t->DzCorrected(*fBeamSpot) - ZVC ) > 0.2) continue;
|
| 178 |
|
|
//if(fabs(t->DzCorrected(*fBeamSpot) - ZVC ) > 3*zwidth) continue;
|
| 179 |
|
|
if(pfca->Pt()<0.3 || pfca->Pt()>200) continue;
|
| 180 |
maxi |
1.3 |
|
| 181 |
|
|
std::vector<const Track*>::iterator itt;
|
| 182 |
|
|
itt = find ((vtool->excluded).begin(), (vtool->excluded).end(), t);
|
| 183 |
|
|
if(itt != (vtool->excluded).end()) continue;
|
| 184 |
maxi |
1.2 |
|
| 185 |
|
|
//Float_t phid = phim - t->Phi();
|
| 186 |
|
|
Float_t phid = bosophi+3.14 - t->Phi();
|
| 187 |
|
|
ymean = ymean + pfca->Pt()*TMath::Sin(phid);
|
| 188 |
|
|
ymsq = ymsq + pow(pfca->Pt()*TMath::Sin(phid),2);
|
| 189 |
|
|
bosoproj = bosoproj + pfca->Pt()*TMath::Cos(bosophi-t->Phi());
|
| 190 |
|
|
if(TMath::Sin(phid) > 0.) ortplus = ortplus + pfca->Pt()*TMath::Sin(phid);
|
| 191 |
|
|
else ortminus = ortminus + pfca->Pt()*fabs(TMath::Sin(phid));
|
| 192 |
|
|
if(TMath::Cos(phid) > 0.) bdplus = bdplus + pfca->Pt()*TMath::Cos(phid);
|
| 193 |
|
|
else bdminus = bdminus + pfca->Pt()*fabs(TMath::Cos(phid));
|
| 194 |
|
|
if(TMath::Cos(phid) > 0.) ntplus = ntplus + TMath::Cos(phid);
|
| 195 |
|
|
}
|
| 196 |
bendavid |
1.1 |
|
| 197 |
maxi |
1.2 |
Float_t phimsq = sqrt(ymsq/ntrks - pow(ymean/ntrks,2))*180/3.14;
|
| 198 |
fabstoec |
1.4 |
// not used... commented out by Fabian
|
| 199 |
|
|
//Double_t A2n = (ortplus+ortminus) > 0.01 ? (bdplus+bdminus)/(ortplus+ortminus) : (bdplus+bdminus)/0.01;
|
| 200 |
maxi |
1.2 |
Double_t A1n = bdplus/bosopt;
|
| 201 |
|
|
Double_t A3n = ntplus > 0 ? bdplus/ntplus : 0.;
|
| 202 |
|
|
//-------------------
|
| 203 |
|
|
Double_t angle = 180/3.14*TVector2::Phi_0_2pi(bosophi-phim);
|
| 204 |
bendavid |
1.1 |
|
| 205 |
maxi |
1.2 |
vtool->tmvar1 = ntrks;
|
| 206 |
|
|
vtool->tmvar2 = sumpt;
|
| 207 |
|
|
vtool->tmvar3 = A1n;
|
| 208 |
|
|
vtool->tmvar4 = angle;
|
| 209 |
|
|
vtool->tmvar5 = phimsq;
|
| 210 |
|
|
vtool->tmvar6 = A3n;
|
| 211 |
|
|
|
| 212 |
|
|
double p1 = vtool->reader->GetProba ( "BDTD method" );
|
| 213 |
|
|
double p2 = vtool->reader->GetProba ( "BDTG method" );
|
| 214 |
|
|
double p3 = vtool->reader->GetProba ( "MLP method");
|
| 215 |
|
|
double p4 = vtool->reader->GetProba ( "MLPBFGS method");
|
| 216 |
|
|
|
| 217 |
|
|
double retval = p1*p2*p3*p4;
|
| 218 |
|
|
|
| 219 |
|
|
return retval;
|
| 220 |
|
|
}
|
| 221 |
|
|
|
| 222 |
|
|
double VertexTools::VertexWidth(const Vertex* vert, const BaseVertex *fBeamSpot){
|
| 223 |
|
|
double width = 0.;
|
| 224 |
|
|
double zmean = 0.;
|
| 225 |
|
|
double zmeansq = 0.;
|
| 226 |
|
|
double ww = 0.;
|
| 227 |
|
|
if(vert == NULL) return width;
|
| 228 |
|
|
for(unsigned i = 0; i < vert->NTracks(); i++){
|
| 229 |
|
|
const Track *t = vert->Trk(i);
|
| 230 |
|
|
if(t->Pt() < 0.3 ) continue;
|
| 231 |
|
|
if(t->Pt() > 500.) continue;
|
| 232 |
|
|
ww = ww + t->Pt();
|
| 233 |
|
|
zmean = zmean + t->Pt()*t->DzCorrected(*fBeamSpot);
|
| 234 |
|
|
zmeansq = zmeansq + t->Pt()*t->DzCorrected(*fBeamSpot)*t->DzCorrected(*fBeamSpot);
|
| 235 |
|
|
}
|
| 236 |
|
|
if( !(ww > 0.) ) return width;
|
| 237 |
|
|
zmean = zmean/ww;
|
| 238 |
|
|
|
| 239 |
|
|
width = sqrt(zmeansq/ww - zmean*zmean);
|
| 240 |
|
|
return width;
|
| 241 |
|
|
|
| 242 |
|
|
}
|
| 243 |
maxi |
1.3 |
|
| 244 |
|
|
void VertexTools::BanThisTrack(const Track* track){
|
| 245 |
|
|
VertexTools* vtool = VertexTools::instance("");
|
| 246 |
|
|
(vtool->excluded).push_back(track);
|
| 247 |
|
|
}
|
| 248 |
|
|
void VertexTools::Reset(){
|
| 249 |
|
|
VertexTools* vtool = VertexTools::instance("");
|
| 250 |
|
|
(vtool->excluded).clear();
|
| 251 |
|
|
}
|
| 252 |
fabstoec |
1.5 |
|
| 253 |
|
|
|
| 254 |
|
|
//------------------------------------------------------------------------------------
|
| 255 |
|
|
// The below tools are from the H->2photons EPS2011 BaseLine (common) Selection
|
| 256 |
|
|
const Vertex* VertexTools::findVtxBasicRanking(const Photon* ph1,
|
| 257 |
|
|
const Photon* ph2,
|
| 258 |
|
|
const BaseVertex* bsp,
|
| 259 |
|
|
const VertexCol* vtcs,
|
| 260 |
|
|
const DecayParticleCol* conv) {
|
| 261 |
|
|
|
| 262 |
|
|
// check if all input is valid
|
| 263 |
|
|
if( !ph1 || !ph2 || !bsp || !vtcs ) return NULL;
|
| 264 |
|
|
// CAUTION: We allow for passing NULL for the Conversions, in that case only the simple Ranking is used.
|
| 265 |
|
|
|
| 266 |
|
|
// here we will store the idx of the best Vtx
|
| 267 |
|
|
unsigned int bestIdx = 0;
|
| 268 |
|
|
|
| 269 |
|
|
// using asd much as possible 'Globe' naming schemes...
|
| 270 |
|
|
int* ptbal_rank = new int [vtcs->GetEntries()];
|
| 271 |
|
|
int* ptasym_rank = new int [vtcs->GetEntries()];
|
| 272 |
|
|
int* total_rank = new int [vtcs->GetEntries()];
|
| 273 |
|
|
double* ptbal = new double[vtcs->GetEntries()];
|
| 274 |
|
|
double* ptasym = new double[vtcs->GetEntries()];
|
| 275 |
|
|
|
| 276 |
|
|
unsigned int numVertices = vtcs->GetEntries();
|
| 277 |
|
|
double ptgg = 0.; // stored for later in the conversion
|
| 278 |
|
|
|
| 279 |
|
|
// loop over all the vertices...
|
| 280 |
|
|
for(unsigned int iVtx = 0; iVtx < numVertices; ++iVtx) {
|
| 281 |
|
|
|
| 282 |
|
|
const Vertex* tVtx = vtcs->At(iVtx);
|
| 283 |
|
|
ptbal [iVtx] = 0.0;
|
| 284 |
|
|
ptasym [iVtx] = 0.0;
|
| 285 |
|
|
ptbal_rank [iVtx] = 1;
|
| 286 |
|
|
ptasym_rank[iVtx] = 1;
|
| 287 |
|
|
|
| 288 |
|
|
// compute the photon momenta with respect to this Vtx
|
| 289 |
|
|
FourVectorM newMomFst = ph1->MomVtx(tVtx->Position());
|
| 290 |
|
|
FourVectorM newMomSec = ph2->MomVtx(tVtx->Position());
|
| 291 |
|
|
|
| 292 |
|
|
FourVectorM higgsMom = newMomFst+newMomSec;
|
| 293 |
|
|
|
| 294 |
|
|
double ph1Eta = newMomFst.Eta();
|
| 295 |
|
|
double ph2Eta = newMomSec.Eta();
|
| 296 |
|
|
|
| 297 |
|
|
double ph1Phi = newMomFst.Phi();
|
| 298 |
|
|
double ph2Phi = newMomSec.Phi();
|
| 299 |
|
|
|
| 300 |
|
|
// loop over all tracks and computew variables for ranking...
|
| 301 |
|
|
FourVectorM totTrkMom(0,0,0,0);
|
| 302 |
|
|
for(unsigned int iTrk = 0; iTrk < tVtx->NTracks(); ++iTrk) {
|
| 303 |
|
|
const Track* tTrk = tVtx->Trk(iTrk);
|
| 304 |
|
|
|
| 305 |
|
|
// compute distance between Trk and the Photons
|
| 306 |
|
|
double tEta = tTrk->Eta();
|
| 307 |
|
|
double tPhi = tTrk->Phi();
|
| 308 |
|
|
double dEta1 = TMath::Abs(tEta-ph1Eta);
|
| 309 |
|
|
double dEta2 = TMath::Abs(tEta-ph2Eta);
|
| 310 |
|
|
double dPhi1 = TMath::Abs(tPhi-ph1Phi);
|
| 311 |
|
|
double dPhi2 = TMath::Abs(tPhi-ph2Phi);
|
| 312 |
|
|
if(dPhi1 > M_PI) dPhi1 = 2*M_PI - dPhi1;
|
| 313 |
|
|
if(dPhi2 > M_PI) dPhi2 = 2*M_PI - dPhi2;
|
| 314 |
|
|
|
| 315 |
|
|
double dR1 = TMath::Sqrt(dEta1*dEta1+dPhi1*dPhi1);
|
| 316 |
|
|
double dR2 = TMath::Sqrt(dEta2*dEta2+dPhi2*dPhi2);
|
| 317 |
|
|
|
| 318 |
|
|
if(dR1 < 0.05 || dR2 < 0.05) continue;
|
| 319 |
|
|
totTrkMom = totTrkMom + tTrk->Mom4(0);
|
| 320 |
|
|
}
|
| 321 |
|
|
|
| 322 |
|
|
// compute the ranking variables for this Vtx
|
| 323 |
|
|
double ptvtx = totTrkMom.Pt();
|
| 324 |
|
|
double pthiggs = higgsMom.Pt();
|
| 325 |
|
|
ptbal [iVtx] = (totTrkMom.Px()*(newMomFst.Px()+newMomSec.Px()));
|
| 326 |
|
|
ptbal [iVtx] += (totTrkMom.Py()*(newMomFst.Py()+newMomSec.Py()));
|
| 327 |
|
|
ptbal [iVtx] = -ptbal[iVtx]/pthiggs;
|
| 328 |
|
|
ptasym[iVtx] = (ptvtx - pthiggs)/(ptvtx + pthiggs);
|
| 329 |
|
|
|
| 330 |
|
|
// do the little ranking acrobatics...
|
| 331 |
|
|
for(unsigned int cVtx =0; cVtx < iVtx; ++cVtx) {
|
| 332 |
|
|
if(ptbal [iVtx] > ptbal [cVtx])
|
| 333 |
|
|
ptbal_rank[cVtx]++;
|
| 334 |
|
|
else
|
| 335 |
|
|
ptbal_rank[iVtx]++;
|
| 336 |
|
|
if(ptasym [iVtx] > ptasym [cVtx])
|
| 337 |
|
|
ptasym_rank[cVtx]++;
|
| 338 |
|
|
else
|
| 339 |
|
|
ptasym_rank[iVtx]++;
|
| 340 |
|
|
}
|
| 341 |
|
|
}
|
| 342 |
|
|
|
| 343 |
|
|
// loop again over all Vertcices (*sigh*), compute total score and final rank
|
| 344 |
|
|
// CAUTION: Total rank starts at '0', so the best ranked Vtx has RANK 0
|
| 345 |
|
|
for(unsigned int iVtx = 0; iVtx < numVertices; ++iVtx) {
|
| 346 |
|
|
ptasym_rank [iVtx] = ptbal_rank [iVtx]*ptasym_rank [iVtx]*(iVtx+1);
|
| 347 |
|
|
total_rank [iVtx] = 0;
|
| 348 |
|
|
for(unsigned int cVtx =0; cVtx < iVtx; ++cVtx) {
|
| 349 |
|
|
if(ptasym_rank [iVtx] > ptasym_rank [cVtx]) total_rank[iVtx]++;
|
| 350 |
|
|
else if(ptasym_rank [iVtx] == ptasym_rank [cVtx]) { // use 'ptbal' as the tie-breaker
|
| 351 |
|
|
if(ptbal_rank [iVtx] > ptbal_rank [cVtx]) total_rank[iVtx]++;
|
| 352 |
|
|
else total_rank[cVtx]++;
|
| 353 |
|
|
}
|
| 354 |
|
|
else total_rank[cVtx]++;
|
| 355 |
|
|
}
|
| 356 |
|
|
}
|
| 357 |
|
|
|
| 358 |
|
|
// delete the auxiliary dynamic arrays
|
| 359 |
|
|
delete[] ptbal_rank ;
|
| 360 |
|
|
delete[] ptasym_rank ;
|
| 361 |
|
|
delete[] ptbal ;
|
| 362 |
|
|
delete[] ptasym ;
|
| 363 |
|
|
|
| 364 |
|
|
// find the best ranked Vertex so far....
|
| 365 |
|
|
for(unsigned int iVtx = 0; iVtx < numVertices; ++iVtx) {
|
| 366 |
|
|
if(total_rank[iVtx] == 0) bestIdx = iVtx;
|
| 367 |
|
|
}
|
| 368 |
|
|
|
| 369 |
|
|
// check if there's a conversion among the pre-selected photons
|
| 370 |
|
|
// ...this will return NULL in case conv==NULL
|
| 371 |
|
|
const DecayParticle* conv1 = PhotonTools::MatchedCiCConversion(ph1, conv, 0.1, 0.1, 0.1);
|
| 372 |
|
|
const DecayParticle* conv2 = PhotonTools::MatchedCiCConversion(ph2, conv, 0.1, 0.1, 0.1);
|
| 373 |
|
|
|
| 374 |
|
|
if( conv1 && ( conv1->Prob() < 0.0005) ) conv1 = NULL;
|
| 375 |
|
|
if( conv2 && ( conv2->Prob() < 0.0005) ) conv2 = NULL;
|
| 376 |
|
|
|
| 377 |
|
|
double zconv = 0.;
|
| 378 |
|
|
double dzconv = 0.;
|
| 379 |
|
|
|
| 380 |
|
|
//--------------------------------------------------------------------
|
| 381 |
|
|
// start doing the Conversion acrobatics... 'copied' from the Globe...
|
| 382 |
|
|
if(conv1 || conv2) {
|
| 383 |
|
|
if( !conv2 ){
|
| 384 |
|
|
const mithep::ThreeVector caloPos1(ph1->CaloPos());
|
| 385 |
|
|
zconv = conv1->Z0EcalVtx(bsp->Position(), caloPos1);
|
| 386 |
|
|
if( ph1->IsEB() ) {
|
| 387 |
|
|
double rho = conv1->Position().Rho();
|
| 388 |
|
|
if ( rho < 15. ) dzconv = 0.06;
|
| 389 |
|
|
else if( rho < 60. ) dzconv = 0.67;
|
| 390 |
|
|
else dzconv = 2.04;
|
| 391 |
|
|
} else {
|
| 392 |
|
|
double z = conv1->Position().Z();
|
| 393 |
|
|
if ( TMath::Abs(z) < 50. ) dzconv = 0.18;
|
| 394 |
|
|
else if( TMath::Abs(z) < 100.) dzconv = 0.61;
|
| 395 |
|
|
else dzconv = 0.99;
|
| 396 |
|
|
}
|
| 397 |
|
|
} else if( !conv1 ) {
|
| 398 |
|
|
const mithep::ThreeVector caloPos2(ph2->CaloPos());
|
| 399 |
|
|
zconv = conv2->Z0EcalVtx(bsp->Position(), caloPos2);
|
| 400 |
|
|
if( ph2->IsEB() ) {
|
| 401 |
|
|
double rho = conv2->Position().Rho();
|
| 402 |
|
|
if ( rho < 15. ) dzconv = 0.06;
|
| 403 |
|
|
else if( rho < 60. ) dzconv = 0.67;
|
| 404 |
|
|
else dzconv = 2.04;
|
| 405 |
|
|
} else {
|
| 406 |
|
|
double z = conv2->Position().Z();
|
| 407 |
|
|
if ( TMath::Abs(z) < 50. ) dzconv = 0.18;
|
| 408 |
|
|
else if( TMath::Abs(z) < 100.) dzconv = 0.61;
|
| 409 |
|
|
else dzconv = 0.99;
|
| 410 |
|
|
}
|
| 411 |
|
|
} else {
|
| 412 |
|
|
const mithep::ThreeVector caloPos1(ph1->CaloPos());
|
| 413 |
|
|
double z1 = conv1->Z0EcalVtx(bsp->Position(), caloPos1);
|
| 414 |
|
|
double dz1 = 0.;
|
| 415 |
|
|
if( ph1->IsEB() ) {
|
| 416 |
|
|
double rho = conv1->Position().Rho();
|
| 417 |
|
|
if ( rho < 15. ) dz1 = 0.06;
|
| 418 |
|
|
else if( rho < 60. ) dz1 = 0.67;
|
| 419 |
|
|
else dz1 = 2.04;
|
| 420 |
|
|
} else {
|
| 421 |
|
|
double z = conv1->Position().Z();
|
| 422 |
|
|
if ( TMath::Abs(z) < 50. ) dz1 = 0.18;
|
| 423 |
|
|
else if( TMath::Abs(z) < 100.) dz1 = 0.61;
|
| 424 |
|
|
else dz1 = 0.99;
|
| 425 |
|
|
}
|
| 426 |
|
|
const mithep::ThreeVector caloPos2(ph2->CaloPos());
|
| 427 |
|
|
double z2 = conv2->Z0EcalVtx(bsp->Position(), caloPos2);
|
| 428 |
|
|
double dz2 = 0.;
|
| 429 |
|
|
if( ph2->IsEB() ) {
|
| 430 |
|
|
double rho = conv2->Position().Rho();
|
| 431 |
|
|
if ( rho < 15. ) dz2 = 0.06;
|
| 432 |
|
|
else if( rho < 60. ) dz2 = 0.67;
|
| 433 |
|
|
else dz2 = 2.04;
|
| 434 |
|
|
} else {
|
| 435 |
|
|
double z = conv2->Position().Z();
|
| 436 |
|
|
if ( TMath::Abs(z) < 50. ) dz2 = 0.18;
|
| 437 |
|
|
else if( TMath::Abs(z) < 100.) dz2 = 0.61;
|
| 438 |
|
|
else dz2 = 0.99;
|
| 439 |
|
|
}
|
| 440 |
|
|
zconv = ( 1./(1./dz1/dz1 + 1./dz2/dz2 )*(z1/dz1/dz1 + z2/dz2/dz2) ) ; // weighted average
|
| 441 |
|
|
dzconv = TMath::Sqrt( 1./(1./dz1/dz1 + 1./dz2/dz2)) ;
|
| 442 |
|
|
}
|
| 443 |
|
|
|
| 444 |
|
|
// loop over all ranked Vertices and choose the closest to the Conversion one
|
| 445 |
|
|
int maxVertices = ( ptgg > 30 ? 3 : 5);
|
| 446 |
|
|
double minDz = -1.;
|
| 447 |
|
|
for(unsigned int iVtx =0; iVtx < numVertices; ++iVtx) {
|
| 448 |
|
|
if(total_rank[iVtx] < maxVertices) {
|
| 449 |
|
|
const Vertex* tVtx = vtcs->At(iVtx);
|
| 450 |
|
|
double tDz = TMath::Abs(zconv - tVtx->Z());
|
| 451 |
|
|
if( (minDz < 0. || tDz < minDz) && ( tDz < dzconv ) ) {
|
| 452 |
|
|
minDz = tDz;
|
| 453 |
|
|
bestIdx = iVtx;
|
| 454 |
|
|
}
|
| 455 |
|
|
}
|
| 456 |
|
|
}
|
| 457 |
|
|
}
|
| 458 |
|
|
// END of Conversion Acrobatics
|
| 459 |
|
|
//--------------------------------------------------------------------
|
| 460 |
|
|
|
| 461 |
|
|
delete[] total_rank ;
|
| 462 |
|
|
return vtcs->At(bestIdx);
|
| 463 |
|
|
}
|
