| 1 |
#include "MitPhysics/Mods/interface/PhotonCiCMod.h"
|
| 2 |
#include "MitAna/DataTree/interface/PhotonCol.h"
|
| 3 |
#include "MitPhysics/Init/interface/ModNames.h"
|
| 4 |
#include "MitPhysics/Utils/interface/IsolationTools.h"
|
| 5 |
#include "MitPhysics/Utils/interface/PhotonTools.h"
|
| 6 |
#include <TNtuple.h>
|
| 7 |
#include <TRandom3.h>
|
| 8 |
|
| 9 |
using namespace mithep;
|
| 10 |
|
| 11 |
ClassImp(mithep::PhotonCiCMod)
|
| 12 |
|
| 13 |
//--------------------------------------------------------------------------------------------------
|
| 14 |
PhotonCiCMod::PhotonCiCMod(const char *name, const char *title) :
|
| 15 |
BaseMod(name,title),
|
| 16 |
fPhotonBranchName (Names::gkPhotonBrn),
|
| 17 |
fGoodPhotonsName (ModNames::gkGoodPhotonsName),
|
| 18 |
fTrackBranchName (Names::gkTrackBrn),
|
| 19 |
fPileUpDenName (Names::gkPileupEnergyDensityBrn),
|
| 20 |
|
| 21 |
fElectronName ("Electrons"),
|
| 22 |
fPhotonPtMin (20.0),
|
| 23 |
fApplySpikeRemoval (kFALSE),
|
| 24 |
|
| 25 |
fAbsEtaMax (999.99),
|
| 26 |
|
| 27 |
fPhotons(0),
|
| 28 |
fTracks(0),
|
| 29 |
fPileUpDen(0),
|
| 30 |
fElectrons(0),
|
| 31 |
|
| 32 |
fPVName(Names::gkPVBeamSpotBrn),
|
| 33 |
fPV(0),
|
| 34 |
fPVFromBranch(kTRUE),
|
| 35 |
|
| 36 |
fConversions(0),
|
| 37 |
fConversionName(Names::gkMvfConversionBrn),
|
| 38 |
|
| 39 |
fBeamspot(0),
|
| 40 |
|
| 41 |
// May10 ReReco
|
| 42 |
fDataEnCorr_EB_hR9(0),
|
| 43 |
fDataEnCorr_EB_lR9(0),
|
| 44 |
fDataEnCorr_EE_hR9(0),
|
| 45 |
fDataEnCorr_EE_lR9(0),
|
| 46 |
|
| 47 |
fRunStart(0),
|
| 48 |
fRunEnd(0),
|
| 49 |
|
| 50 |
fMCSmear_EB_hR9(0.0089),
|
| 51 |
fMCSmear_EB_lR9(0.0199),
|
| 52 |
fMCSmear_EE_hR9(0.0409),
|
| 53 |
fMCSmear_EE_lR9(0.0246),
|
| 54 |
|
| 55 |
fIsData(false),
|
| 56 |
|
| 57 |
rng(new TRandom3()),
|
| 58 |
|
| 59 |
fMCParticleName(Names::gkMCPartBrn),
|
| 60 |
fMCParticles(0),
|
| 61 |
fPileUpName ("PileupInfo"),
|
| 62 |
fPileUp (0)
|
| 63 |
|
| 64 |
{
|
| 65 |
// Constructor.
|
| 66 |
}
|
| 67 |
|
| 68 |
PhotonCiCMod::~PhotonCiCMod(){
|
| 69 |
if(rng) delete rng;
|
| 70 |
}
|
| 71 |
|
| 72 |
//--------------------------------------------------------------------------------------------------
|
| 73 |
void PhotonCiCMod::Process()
|
| 74 |
{
|
| 75 |
// Process entries of the tree.
|
| 76 |
|
| 77 |
LoadEventObject(fPhotonBranchName, fPhotons);
|
| 78 |
|
| 79 |
PhotonOArr *GoodPhotons = new PhotonOArr;
|
| 80 |
GoodPhotons->SetName(fGoodPhotonsName);
|
| 81 |
GoodPhotons->SetOwner(kTRUE);
|
| 82 |
|
| 83 |
Double_t _tRho = 0.;
|
| 84 |
LoadEventObject(fTrackBranchName, fTracks);
|
| 85 |
LoadEventObject(fPileUpDenName, fPileUpDen);
|
| 86 |
LoadEventObject(fElectronName, fElectrons);
|
| 87 |
LoadEventObject(fPVName, fPV);
|
| 88 |
LoadEventObject(fConversionName, fConversions);
|
| 89 |
LoadBranch("BeamSpot");
|
| 90 |
|
| 91 |
if( !fIsData ) {
|
| 92 |
LoadBranch(fMCParticleName);
|
| 93 |
LoadBranch(fPileUpName);
|
| 94 |
}
|
| 95 |
|
| 96 |
Float_t numPU = -1.;
|
| 97 |
if( !fIsData )
|
| 98 |
numPU = (Float_t) fPileUp->At(0)->GetPU_NumInteractions();
|
| 99 |
|
| 100 |
if(fPileUpDen->GetEntries() > 0)
|
| 101 |
_tRho = (Double_t) fPileUpDen->At(0)->RhoRandomLowEta();
|
| 102 |
|
| 103 |
bool doVtxSelection = true;
|
| 104 |
bool doMCSmear = true;
|
| 105 |
|
| 106 |
const EventHeader* evtHead = this->GetEventHeader();
|
| 107 |
|
| 108 |
unsigned int evtNum = evtHead->EvtNum();
|
| 109 |
Float_t _evtNum1 = (Float_t) ( (int) (evtNum/10000.) );
|
| 110 |
Float_t _evtNum2 = (Float_t) ( (int) (evtNum % 10000) );
|
| 111 |
|
| 112 |
//double evtNumTest = (int) ( ( (double) _evtNum1 )*10000. + (double) _evtNum2 );
|
| 113 |
|
| 114 |
UInt_t runNumber = evtHead->RunNum();
|
| 115 |
Float_t _runNum = (Float_t) runNumber;
|
| 116 |
Float_t _lumiSec = (Float_t) evtHead->LumiSec();
|
| 117 |
|
| 118 |
|
| 119 |
//unsigned int numVertices = fPV->GetEntries();
|
| 120 |
|
| 121 |
const BaseVertex *bsp = dynamic_cast<const BaseVertex*>(fBeamspot->At(0));
|
| 122 |
|
| 123 |
PhotonOArr* preselPh = new PhotonOArr;
|
| 124 |
|
| 125 |
// 1. we do the pre-selection; but keep the non-passing photons in a secont container...
|
| 126 |
for (UInt_t i=0; i<fPhotons->GetEntries(); ++i) {
|
| 127 |
const Photon *ph = fPhotons->At(i);
|
| 128 |
|
| 129 |
if(ph->SCluster()->AbsEta()>=2.5 || (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566)) continue;
|
| 130 |
if(ph->Et() < 20.) continue;
|
| 131 |
if(ph->HadOverEm() > 0.15) continue;
|
| 132 |
if(ph->AbsEta() < 1.5) {
|
| 133 |
if(ph->CoviEtaiEta() > 0.013) continue;
|
| 134 |
} else {
|
| 135 |
if(ph->CoviEtaiEta() > 0.03) continue;
|
| 136 |
}
|
| 137 |
|
| 138 |
Bool_t passSpikeRemovalFilter = kTRUE;
|
| 139 |
|
| 140 |
if (ph->SCluster() && ph->SCluster()->Seed()) {
|
| 141 |
if(ph->SCluster()->Seed()->Energy() > 5.0 &&
|
| 142 |
ph->SCluster()->Seed()->EMax() / ph->SCluster()->Seed()->E3x3() > 0.95 )
|
| 143 |
passSpikeRemovalFilter = kFALSE;
|
| 144 |
}
|
| 145 |
if (fApplySpikeRemoval && !passSpikeRemovalFilter) continue;
|
| 146 |
|
| 147 |
preselPh->Add(ph);
|
| 148 |
}
|
| 149 |
|
| 150 |
// sort both by pt... again ;)
|
| 151 |
preselPh->Sort();
|
| 152 |
|
| 153 |
unsigned int numPairs = 0;
|
| 154 |
if( preselPh->GetEntries() > 0)
|
| 155 |
numPairs = (preselPh->GetEntries()-1)*preselPh->GetEntries()/2;
|
| 156 |
|
| 157 |
// create all possible pairs of pre-selected photons
|
| 158 |
|
| 159 |
std::vector<unsigned int> idxFst;
|
| 160 |
std::vector<unsigned int> idxSec;
|
| 161 |
std::vector<bool> pairPasses;
|
| 162 |
|
| 163 |
if(numPairs > 0) {
|
| 164 |
for(unsigned int iFst = 0; iFst <preselPh->GetEntries() - 1; ++iFst) {
|
| 165 |
for(unsigned int iSec = iFst + 1; iSec <preselPh->GetEntries(); ++iSec) {
|
| 166 |
idxFst.push_back(iFst);
|
| 167 |
idxSec.push_back(iSec);
|
| 168 |
pairPasses.push_back(true);
|
| 169 |
}
|
| 170 |
}
|
| 171 |
}
|
| 172 |
|
| 173 |
// array to store the index of 'chosen Vtx' for each pair
|
| 174 |
const Vertex** theVtx = new const Vertex*[numPairs];
|
| 175 |
UInt_t* theVtxIdx = new UInt_t[numPairs];
|
| 176 |
|
| 177 |
// arays to store the Vtx 'fixed' photons
|
| 178 |
Photon** fixPhFst = new Photon*[numPairs];
|
| 179 |
Photon** fixPhSec = new Photon*[numPairs];
|
| 180 |
|
| 181 |
// sotr pair-indices for pairs passing the selection
|
| 182 |
std::vector<unsigned int> passPairs;
|
| 183 |
passPairs.resize(0);
|
| 184 |
|
| 185 |
unsigned int theChosenVtx = 0;
|
| 186 |
|
| 187 |
float kinPh1[20];
|
| 188 |
float kinPh2[20];
|
| 189 |
|
| 190 |
for(int i=0; i<10; ++i) {
|
| 191 |
kinPh1[i] =-99.;
|
| 192 |
kinPh2[i] =-99.;
|
| 193 |
}
|
| 194 |
|
| 195 |
float ptBefore1 = -99.;
|
| 196 |
float ptBefore2 = -99.;
|
| 197 |
|
| 198 |
bool print = false;
|
| 199 |
if(evtNum == 17031 && false) {
|
| 200 |
std::cout<<" ------------------------------------------- "<<std::endl;
|
| 201 |
std::cout<<" printing info for event #"<<evtNum<<std::endl;
|
| 202 |
print = true;
|
| 203 |
}
|
| 204 |
|
| 205 |
for(unsigned int iPair = 0; iPair < numPairs; ++iPair) {
|
| 206 |
|
| 207 |
// copy the photons for manipulation
|
| 208 |
fixPhFst[iPair] = new Photon(*preselPh->At(idxFst[iPair]));
|
| 209 |
fixPhSec[iPair] = new Photon(*preselPh->At(idxSec[iPair]));
|
| 210 |
|
| 211 |
// if this is Data, scale the energy, if MC smear...
|
| 212 |
FourVectorM scMomFst = fixPhFst[iPair]->Mom();
|
| 213 |
FourVectorM scMomSec = fixPhSec[iPair]->Mom();
|
| 214 |
double scaleFac1 = 1.;
|
| 215 |
double scaleFac2 = 1.;
|
| 216 |
if (fIsData) {
|
| 217 |
if( fRunStart.size() > 0) {
|
| 218 |
// find run in rnage
|
| 219 |
Int_t runRange=-1;
|
| 220 |
for(UInt_t iRun = 0; iRun<fRunStart.size(); ++iRun) {
|
| 221 |
if( runNumber >= fRunStart[iRun] && runNumber <= fRunEnd[iRun]) {
|
| 222 |
runRange = (Int_t) iRun;
|
| 223 |
break;
|
| 224 |
}
|
| 225 |
}
|
| 226 |
if(runRange > -1) {
|
| 227 |
if(fixPhFst[iPair]->IsEB())
|
| 228 |
if(fixPhFst[iPair]->R9() > 0.94)
|
| 229 |
scaleFac1 += fDataEnCorr_EB_hR9[runRange];
|
| 230 |
else
|
| 231 |
scaleFac1 += fDataEnCorr_EB_lR9[runRange];
|
| 232 |
else
|
| 233 |
if(fixPhFst[iPair]->R9() > 0.94)
|
| 234 |
scaleFac1 += fDataEnCorr_EE_hR9[runRange];
|
| 235 |
else
|
| 236 |
scaleFac1 += fDataEnCorr_EE_lR9[runRange];
|
| 237 |
if(fixPhSec[iPair]->IsEB())
|
| 238 |
if(fixPhSec[iPair]->R9() > 0.94)
|
| 239 |
scaleFac2 += fDataEnCorr_EB_hR9[runRange];
|
| 240 |
else
|
| 241 |
scaleFac2 += fDataEnCorr_EB_lR9[runRange];
|
| 242 |
else
|
| 243 |
if(fixPhSec[iPair]->R9() > 0.94)
|
| 244 |
scaleFac2 += fDataEnCorr_EE_hR9[runRange];
|
| 245 |
else
|
| 246 |
scaleFac2 += fDataEnCorr_EE_lR9[runRange];
|
| 247 |
}
|
| 248 |
}
|
| 249 |
} else {
|
| 250 |
// get the smearing for MC photons..
|
| 251 |
UInt_t seedBase = (UInt_t) evtNum + (UInt_t) _runNum + (UInt_t) _lumiSec;
|
| 252 |
UInt_t seed1 = seedBase + (UInt_t) fixPhFst[iPair]->E() + (UInt_t) (TMath::Abs(10.*fixPhFst[iPair]->SCluster()->Eta()));
|
| 253 |
UInt_t seed2 = seedBase + (UInt_t) fixPhSec[iPair]->E() + (UInt_t) (TMath::Abs(10.*fixPhSec[iPair]->SCluster()->Eta()));
|
| 254 |
|
| 255 |
double width1 = 0.;
|
| 256 |
double width2 = 0.;
|
| 257 |
if(fixPhFst[iPair]->IsEB())
|
| 258 |
if(fixPhFst[iPair]->R9() > 0.94)
|
| 259 |
width1 = fMCSmear_EB_hR9;
|
| 260 |
else
|
| 261 |
width1 = fMCSmear_EB_lR9;
|
| 262 |
else
|
| 263 |
if(fixPhFst[iPair]->R9() > 0.94)
|
| 264 |
width1 = fMCSmear_EE_hR9;
|
| 265 |
else
|
| 266 |
width1 = fMCSmear_EE_lR9;
|
| 267 |
|
| 268 |
|
| 269 |
if(fixPhSec[iPair]->IsEB())
|
| 270 |
if(fixPhSec[iPair]->R9() > 0.94)
|
| 271 |
width2 = fMCSmear_EB_hR9;
|
| 272 |
else
|
| 273 |
width2 = fMCSmear_EB_lR9;
|
| 274 |
else
|
| 275 |
if(fixPhSec[iPair]->R9() > 0.94)
|
| 276 |
width2 = fMCSmear_EE_hR9;
|
| 277 |
else
|
| 278 |
width2 = fMCSmear_EE_lR9;
|
| 279 |
|
| 280 |
if(doMCSmear) {
|
| 281 |
rng->SetSeed(seed1);
|
| 282 |
scaleFac1 = rng->Gaus(1.,width1);
|
| 283 |
rng->SetSeed(seed2);
|
| 284 |
scaleFac2 = rng->Gaus(1.,width2);
|
| 285 |
}
|
| 286 |
}
|
| 287 |
|
| 288 |
if(iPair==0) {
|
| 289 |
ptBefore1 = fixPhFst[iPair]->Pt();
|
| 290 |
ptBefore2 = fixPhSec[iPair]->Pt();
|
| 291 |
}
|
| 292 |
|
| 293 |
if(print && false) {
|
| 294 |
std::cout<<" Photon Pair #"<<iPair+1<<std::endl;
|
| 295 |
std::cout<<" Ph1 px = "<<fixPhFst[iPair]->Mom().X()<<std::endl;
|
| 296 |
std::cout<<" py = "<<fixPhFst[iPair]->Mom().Y()<<std::endl;
|
| 297 |
std::cout<<" pz = "<<fixPhFst[iPair]->Mom().Z()<<std::endl;
|
| 298 |
std::cout<<" E = "<<fixPhFst[iPair]->Mom().E()<<std::endl;
|
| 299 |
std::cout<<" M = "<<fixPhFst[iPair]->Mom().M()<<std::endl;
|
| 300 |
std::cout<<" Ph2 px = "<<fixPhSec[iPair]->Mom().X()<<std::endl;
|
| 301 |
std::cout<<" py = "<<fixPhSec[iPair]->Mom().Y()<<std::endl;
|
| 302 |
std::cout<<" pz = "<<fixPhSec[iPair]->Mom().Z()<<std::endl;
|
| 303 |
std::cout<<" E = "<<fixPhSec[iPair]->Mom().E()<<std::endl;
|
| 304 |
std::cout<<" M = "<<fixPhSec[iPair]->Mom().M()<<std::endl;
|
| 305 |
}
|
| 306 |
|
| 307 |
fixPhFst[iPair]->SetMom(scaleFac1*scMomFst.X(), scaleFac1*scMomFst.Y(), scaleFac1*scMomFst.Z(), scaleFac1*scMomFst.E());
|
| 308 |
fixPhSec[iPair]->SetMom(scaleFac2*scMomSec.X(), scaleFac2*scMomSec.Y(), scaleFac2*scMomSec.Z(), scaleFac2*scMomSec.E());
|
| 309 |
|
| 310 |
if(print && false) {
|
| 311 |
std::cout<<" SF 1 = "<<scaleFac1<<std::endl;
|
| 312 |
std::cout<<" SF 2 = "<<scaleFac2<<std::endl;
|
| 313 |
|
| 314 |
std::cout<<" Photon Pair #"<<iPair+1<<std::endl;
|
| 315 |
std::cout<<" Ph1 px = "<<fixPhFst[iPair]->Mom().X()<<std::endl;
|
| 316 |
std::cout<<" py = "<<fixPhFst[iPair]->Mom().Y()<<std::endl;
|
| 317 |
std::cout<<" pz = "<<fixPhFst[iPair]->Mom().Z()<<std::endl;
|
| 318 |
std::cout<<" E = "<<fixPhFst[iPair]->Mom().E()<<std::endl;
|
| 319 |
std::cout<<" M = "<<fixPhFst[iPair]->Mom().M()<<std::endl;
|
| 320 |
std::cout<<" Ph2 px = "<<fixPhSec[iPair]->Mom().X()<<std::endl;
|
| 321 |
std::cout<<" py = "<<fixPhSec[iPair]->Mom().Y()<<std::endl;
|
| 322 |
std::cout<<" pz = "<<fixPhSec[iPair]->Mom().Z()<<std::endl;
|
| 323 |
std::cout<<" E = "<<fixPhSec[iPair]->Mom().E()<<std::endl;
|
| 324 |
std::cout<<" M = "<<fixPhSec[iPair]->Mom().M()<<std::endl;
|
| 325 |
}
|
| 326 |
|
| 327 |
// store the vertex for this pair
|
| 328 |
if(doVtxSelection) {
|
| 329 |
unsigned int iVtx = findBestVertex(fixPhFst[iPair],fixPhSec[iPair],bsp, print);
|
| 330 |
theVtx[iPair] = fPV->At(iVtx);
|
| 331 |
theVtxIdx[iPair] = iVtx;
|
| 332 |
if(iPair == 0) theChosenVtx = iVtx;
|
| 333 |
} else
|
| 334 |
theVtx[iPair] = fPV->At(0);
|
| 335 |
|
| 336 |
|
| 337 |
// fix the kinematics for both events
|
| 338 |
FourVectorM newMomFst = fixPhFst[iPair]->MomVtx(theVtx[iPair]->Position());
|
| 339 |
FourVectorM newMomSec = fixPhSec[iPair]->MomVtx(theVtx[iPair]->Position());
|
| 340 |
fixPhFst[iPair]->SetMom(newMomFst.X(), newMomFst.Y(), newMomFst.Z(), newMomFst.E());
|
| 341 |
fixPhSec[iPair]->SetMom(newMomSec.X(), newMomSec.Y(), newMomSec.Z(), newMomSec.E());
|
| 342 |
|
| 343 |
|
| 344 |
if(print && false) {
|
| 345 |
std::cout<<" Vtx = "<<theChosenVtx<<std::endl;
|
| 346 |
std::cout<<" Photon Pair #"<<iPair+1<<std::endl;
|
| 347 |
std::cout<<" Ph1 px = "<<fixPhFst[iPair]->Mom().X()<<std::endl;
|
| 348 |
std::cout<<" py = "<<fixPhFst[iPair]->Mom().Y()<<std::endl;
|
| 349 |
std::cout<<" pz = "<<fixPhFst[iPair]->Mom().Z()<<std::endl;
|
| 350 |
std::cout<<" E = "<<fixPhFst[iPair]->Mom().E()<<std::endl;
|
| 351 |
std::cout<<" M = "<<fixPhFst[iPair]->Mom().M()<<std::endl;
|
| 352 |
std::cout<<" Ph2 px = "<<fixPhSec[iPair]->Mom().X()<<std::endl;
|
| 353 |
std::cout<<" py = "<<fixPhSec[iPair]->Mom().Y()<<std::endl;
|
| 354 |
std::cout<<" pz = "<<fixPhSec[iPair]->Mom().Z()<<std::endl;
|
| 355 |
std::cout<<" E = "<<fixPhSec[iPair]->Mom().E()<<std::endl;
|
| 356 |
std::cout<<" M = "<<fixPhSec[iPair]->Mom().M()<<std::endl;
|
| 357 |
}
|
| 358 |
|
| 359 |
|
| 360 |
|
| 361 |
if(iPair != 0) {
|
| 362 |
// check if both photons pass the CiC selection
|
| 363 |
bool pass1 = PhotonTools::PassCiCSelection(fixPhFst[iPair], theVtx[iPair], fTracks, fElectrons, fPV, _tRho, 40., true, false);
|
| 364 |
bool pass2 = PhotonTools::PassCiCSelection(fixPhSec[iPair], theVtx[iPair], fTracks, fElectrons, fPV, _tRho, 30., true, false);
|
| 365 |
if( pass1 && pass2 )
|
| 366 |
passPairs.push_back(iPair);
|
| 367 |
} else {
|
| 368 |
bool pass1 = PhotonTools::PassCiCSelection(fixPhFst[iPair], theVtx[iPair], fTracks, fElectrons, fPV, _tRho, 40., true, false, kinPh1);
|
| 369 |
bool pass2 = PhotonTools::PassCiCSelection(fixPhSec[iPair], theVtx[iPair], fTracks, fElectrons, fPV, _tRho, 30., true, false, kinPh2);
|
| 370 |
|
| 371 |
if( pass1 && pass2 )
|
| 372 |
passPairs.push_back(iPair);
|
| 373 |
}
|
| 374 |
}
|
| 375 |
|
| 376 |
// loop over all passing pairs and find the one with the highest sum Et
|
| 377 |
Photon* phHard = NULL;
|
| 378 |
Photon* phSoft = NULL;
|
| 379 |
|
| 380 |
const Vertex* _theVtx = NULL;
|
| 381 |
|
| 382 |
double maxSumEt = 0.;
|
| 383 |
for(unsigned int iPair=0; iPair<passPairs.size(); ++iPair){
|
| 384 |
double sumEt = fixPhFst[passPairs[iPair]]->Et();
|
| 385 |
sumEt += fixPhSec[passPairs[iPair]]->Et();
|
| 386 |
if( sumEt > maxSumEt ) {
|
| 387 |
maxSumEt = sumEt;
|
| 388 |
phHard = fixPhFst[passPairs[iPair]];
|
| 389 |
phSoft = fixPhSec[passPairs[iPair]];
|
| 390 |
_theVtx = theVtx[iPair];
|
| 391 |
theChosenVtx = theVtxIdx[iPair];
|
| 392 |
}
|
| 393 |
}
|
| 394 |
|
| 395 |
Float_t _theVtxZ = -999.;
|
| 396 |
Float_t catPh1 = 0.;
|
| 397 |
Float_t catPh2 = 0.;
|
| 398 |
Float_t catEvt = -1.;
|
| 399 |
|
| 400 |
bool ph1IsLowR9 = false;
|
| 401 |
bool ph2IsLowR9 = false;
|
| 402 |
|
| 403 |
if(phHard && phSoft) {
|
| 404 |
GoodPhotons->AddOwned(phHard);
|
| 405 |
GoodPhotons->AddOwned(phSoft);
|
| 406 |
if(_theVtx)
|
| 407 |
_theVtxZ=_theVtx->Position().Z();
|
| 408 |
|
| 409 |
catPh1=1.;
|
| 410 |
catPh2=1.;
|
| 411 |
catEvt=0.;
|
| 412 |
|
| 413 |
if(phHard->SCluster()->AbsEta()>1.5)
|
| 414 |
catPh1=3;
|
| 415 |
if(phHard->R9() < 0.94) {
|
| 416 |
catPh1=catPh1+1;
|
| 417 |
ph1IsLowR9 = true;
|
| 418 |
}
|
| 419 |
|
| 420 |
if(phSoft->SCluster()->AbsEta()>1.5)
|
| 421 |
catPh2=3;
|
| 422 |
|
| 423 |
if(phSoft->R9() < 0.94) {
|
| 424 |
catPh2=catPh2+1;
|
| 425 |
ph2IsLowR9 = true;
|
| 426 |
}
|
| 427 |
|
| 428 |
if(catPh1 > 2.5 || catPh2 > 2.5)
|
| 429 |
catEvt=2.;
|
| 430 |
if(ph1IsLowR9 || ph2IsLowR9)
|
| 431 |
catEvt=catEvt+1.;
|
| 432 |
}
|
| 433 |
|
| 434 |
// sort according to pt
|
| 435 |
GoodPhotons->Sort();
|
| 436 |
|
| 437 |
// add to event for other modules to use
|
| 438 |
AddObjThisEvt(GoodPhotons);
|
| 439 |
|
| 440 |
delete preselPh;
|
| 441 |
|
| 442 |
bool doFill = (phHard && phSoft);
|
| 443 |
Float_t _mass = ( doFill ? (phHard->Mom()+phSoft->Mom()).M() : -100.);
|
| 444 |
Float_t _ptgg = ( doFill ? (phHard->Mom()+phSoft->Mom()).Pt() : -100.);
|
| 445 |
if(_ptgg < 40. && doFill) catEvt = catEvt+4.;
|
| 446 |
|
| 447 |
Float_t _pth = -100.;
|
| 448 |
Float_t _decayZ = -100.;
|
| 449 |
if( !fIsData ) findHiggsPtAndZ(_pth, _decayZ);
|
| 450 |
|
| 451 |
Float_t fillEvent[] = { (float)_tRho,
|
| 452 |
_pth,
|
| 453 |
_decayZ,
|
| 454 |
_theVtxZ,
|
| 455 |
numPU,
|
| 456 |
_mass,
|
| 457 |
_ptgg,
|
| 458 |
_evtNum1,
|
| 459 |
_evtNum2,
|
| 460 |
_runNum,
|
| 461 |
_lumiSec,
|
| 462 |
(float) theChosenVtx,
|
| 463 |
(float) numPairs,
|
| 464 |
catPh1,
|
| 465 |
catPh2,
|
| 466 |
catEvt,
|
| 467 |
kinPh1[0],
|
| 468 |
kinPh1[1],
|
| 469 |
kinPh1[2],
|
| 470 |
kinPh1[3],
|
| 471 |
kinPh1[4],
|
| 472 |
kinPh1[5],
|
| 473 |
kinPh1[6],
|
| 474 |
kinPh1[7],
|
| 475 |
kinPh1[8],
|
| 476 |
kinPh1[9],
|
| 477 |
kinPh1[10],
|
| 478 |
kinPh1[11],
|
| 479 |
kinPh1[12],
|
| 480 |
kinPh1[13],
|
| 481 |
kinPh1[14],
|
| 482 |
kinPh1[15],
|
| 483 |
kinPh1[16],
|
| 484 |
kinPh1[17],
|
| 485 |
kinPh1[18],
|
| 486 |
kinPh1[19],
|
| 487 |
kinPh2[0],
|
| 488 |
kinPh2[1],
|
| 489 |
kinPh2[2],
|
| 490 |
kinPh2[3],
|
| 491 |
kinPh2[4],
|
| 492 |
kinPh2[5],
|
| 493 |
kinPh2[6],
|
| 494 |
kinPh2[7],
|
| 495 |
kinPh2[8],
|
| 496 |
kinPh2[9],
|
| 497 |
kinPh2[10],
|
| 498 |
kinPh2[11],
|
| 499 |
kinPh2[12],
|
| 500 |
kinPh2[13],
|
| 501 |
kinPh2[14],
|
| 502 |
kinPh2[15],
|
| 503 |
kinPh2[16],
|
| 504 |
kinPh2[17],
|
| 505 |
kinPh2[18],
|
| 506 |
kinPh2[19],
|
| 507 |
ptBefore1,
|
| 508 |
ptBefore2
|
| 509 |
};
|
| 510 |
|
| 511 |
|
| 512 |
|
| 513 |
if(_mass > 0.) {
|
| 514 |
//std::cout<<catPh1<<" "<<catPh2<<" "<<catEvt<<" "<<_mass<<std::endl;
|
| 515 |
hCiCTuple->Fill(fillEvent);
|
| 516 |
}
|
| 517 |
|
| 518 |
delete[] theVtx;
|
| 519 |
delete[] theVtxIdx;
|
| 520 |
|
| 521 |
return;
|
| 522 |
|
| 523 |
}
|
| 524 |
|
| 525 |
//--------------------------------------------------------------------------------------------------
|
| 526 |
void PhotonCiCMod::SlaveBegin()
|
| 527 |
{
|
| 528 |
// Run startup code on the computer (slave) doing the actual analysis. Here,
|
| 529 |
// we just request the photon collection branch.
|
| 530 |
|
| 531 |
ReqEventObject(fPhotonBranchName, fPhotons, kTRUE);
|
| 532 |
ReqEventObject(fTrackBranchName, fTracks, kTRUE);
|
| 533 |
ReqEventObject(fElectronName, fElectrons, kTRUE);
|
| 534 |
ReqEventObject(fPileUpDenName, fPileUpDen, kTRUE);
|
| 535 |
ReqEventObject(fPVName, fPV, fPVFromBranch);
|
| 536 |
ReqEventObject(fConversionName, fConversions,kTRUE);
|
| 537 |
ReqBranch("BeamSpot",fBeamspot);
|
| 538 |
|
| 539 |
if (!fIsData) {
|
| 540 |
ReqBranch(fPileUpName, fPileUp);
|
| 541 |
ReqBranch(Names::gkMCPartBrn,fMCParticles);
|
| 542 |
}
|
| 543 |
|
| 544 |
hCiCTuple = new TNtuple("hCiCTuple","hCiCTuple","rho:higgspt:higgsZ:vtxZ:numPU:mass:ptgg:evtnum1:evtnum2:runnum:lumisec:ivtx:npairs:ph1Cat:ph2Cat:evtCat:ph1Iso1:ph1Iso2:ph1Iso3:ph1Cov:ph1HoE:ph1R9:ph1DR:ph1Pt:ph1Eta:ph1Phi:ph1Eiso3:ph1Eiso4:ph1Hiso4:ph1TisoA:ph1TisoW:ph1Tiso:ph1Et:ph1E:ph1Pass:ph1CatDebug:ph2Iso1:ph2Iso2:ph2Iso3:ph2Cov:ph2HoE:ph2R9:ph2DR:ph2Pt:ph2Eta:ph2Phi:ph2Eiso3:ph2Eiso4:ph2Hiso4:ph2TisoA:ph2TisoW:ph2Tiso:ph2Et:ph2E:ph2Pass:ph2CatDebug:ph1UPt:ph2UPt");
|
| 545 |
|
| 546 |
AddOutput(hCiCTuple);
|
| 547 |
|
| 548 |
}
|
| 549 |
|
| 550 |
// return the index of the bext vertex
|
| 551 |
unsigned int PhotonCiCMod::findBestVertex(Photon* ph1, Photon* ph2, const BaseVertex *bsp, bool print) {
|
| 552 |
|
| 553 |
// loop over all vertices and assigne the ranks
|
| 554 |
int* ptbal_rank = new int[fPV->GetEntries()];
|
| 555 |
int* ptasym_rank = new int[fPV->GetEntries()];
|
| 556 |
int* total_rank = new int[fPV->GetEntries()];
|
| 557 |
double* ptbal = new double[fPV->GetEntries()];
|
| 558 |
double* ptasym = new double[fPV->GetEntries()];
|
| 559 |
|
| 560 |
unsigned int numVertices = fPV->GetEntries();
|
| 561 |
|
| 562 |
double ptgg = 0.; // stored for later in the conversion
|
| 563 |
|
| 564 |
if(print && false) {
|
| 565 |
std::cout<<" --------------------------------- "<<std::endl;
|
| 566 |
std::cout<<" looking for Vtx for photon pair "<<std::endl;
|
| 567 |
std::cout<<" pt 1 = "<<ph1->Et()<<std::endl;
|
| 568 |
std::cout<<" pt 2 = "<<ph2->Et()<<std::endl;
|
| 569 |
std::cout<<" among #"<<numVertices<<" Vertices."<<std::endl;
|
| 570 |
}
|
| 571 |
|
| 572 |
|
| 573 |
for(unsigned int iVtx = 0; iVtx < numVertices; ++iVtx) {
|
| 574 |
if(print)
|
| 575 |
std::cout<<std::endl<<" Vertex #"<<iVtx<<std::endl;
|
| 576 |
|
| 577 |
const Vertex* tVtx = fPV->At(iVtx);
|
| 578 |
ptbal [iVtx] = 0.0;
|
| 579 |
ptasym[iVtx] = 0.0;
|
| 580 |
ptbal_rank [iVtx] = 1;
|
| 581 |
ptasym_rank[iVtx] = 1;
|
| 582 |
|
| 583 |
// compute the photon momenta with respect to this Vtx
|
| 584 |
FourVectorM newMomFst = ph1->MomVtx(tVtx->Position());
|
| 585 |
FourVectorM newMomSec = ph2->MomVtx(tVtx->Position());
|
| 586 |
|
| 587 |
FourVectorM higgsMom = newMomFst+newMomSec;
|
| 588 |
|
| 589 |
double ph1Eta = newMomFst.Eta();
|
| 590 |
double ph2Eta = newMomSec.Eta();
|
| 591 |
|
| 592 |
double ph1Phi = newMomFst.Phi();
|
| 593 |
double ph2Phi = newMomSec.Phi();
|
| 594 |
|
| 595 |
if(print && iVtx == 0 && false ) {
|
| 596 |
std::cout<<" new momenta Et1 = "<<newMomFst.Et()<<std::endl;
|
| 597 |
std::cout<<" Eta = "<<newMomFst.Eta()<<std::endl;
|
| 598 |
std::cout<<" Phi = "<<newMomFst.Phi()<<std::endl;
|
| 599 |
std::cout<<" Px = "<<newMomFst.Px()<<std::endl;
|
| 600 |
std::cout<<" Py = "<<newMomFst.Py()<<std::endl;
|
| 601 |
std::cout<<" new momenta Et2 = "<<newMomSec.Et()<<std::endl;
|
| 602 |
std::cout<<" Eta = "<<newMomSec.Eta()<<std::endl;
|
| 603 |
std::cout<<" Phi = "<<newMomSec.Phi()<<std::endl;
|
| 604 |
std::cout<<" Px = "<<newMomSec.Px()<<std::endl;
|
| 605 |
std::cout<<" Py = "<<newMomSec.Py()<<std::endl;
|
| 606 |
}
|
| 607 |
|
| 608 |
FourVectorM totTrkMom;
|
| 609 |
for(unsigned int iTrk = 0; iTrk < tVtx->NTracks(); ++iTrk) {
|
| 610 |
const Track* tTrk = tVtx->Trk(iTrk);
|
| 611 |
//if(tTrk->Pt()<1.) continue;
|
| 612 |
double tEta = tTrk->Eta();
|
| 613 |
double tPhi = tTrk->Phi();
|
| 614 |
double dEta1 = TMath::Abs(tEta-ph1Eta);
|
| 615 |
double dEta2 = TMath::Abs(tEta-ph2Eta);
|
| 616 |
double dPhi1 = TMath::Abs(tPhi-ph1Phi);
|
| 617 |
double dPhi2 = TMath::Abs(tPhi-ph2Phi);
|
| 618 |
if(dPhi1 > M_PI) dPhi1 = 2*M_PI - dPhi1;
|
| 619 |
if(dPhi2 > M_PI) dPhi2 = 2*M_PI - dPhi2;
|
| 620 |
|
| 621 |
double dR1 = TMath::Sqrt(dEta1*dEta1+dPhi1*dPhi1);
|
| 622 |
double dR2 = TMath::Sqrt(dEta2*dEta2+dPhi2*dPhi2);
|
| 623 |
|
| 624 |
if( ( iVtx == 0 || iVtx == 1 ) && print && false) {
|
| 625 |
std::cout<<" Track #"<<iTrk<<std::endl;
|
| 626 |
std::cout<<" pt = "<<tTrk->Pt()<<std::endl;
|
| 627 |
std::cout<<" eta = "<<tTrk->Eta()<<std::endl;
|
| 628 |
std::cout<<" phi = "<<tTrk->Phi()<<std::endl;
|
| 629 |
std::cout<<" px = "<<tTrk->Px()<<std::endl;
|
| 630 |
std::cout<<" py = "<<tTrk->Py()<<std::endl;
|
| 631 |
std::cout<<" dR1 = "<<dR1<<std::endl;
|
| 632 |
std::cout<<" dR2 = "<<dR2<<std::endl;
|
| 633 |
}
|
| 634 |
|
| 635 |
if(dR1 < 0.05 || dR2 < 0.05) continue;
|
| 636 |
|
| 637 |
if(iTrk == 0) totTrkMom = tTrk->Mom4(0);
|
| 638 |
else totTrkMom = totTrkMom + tTrk->Mom4(0);
|
| 639 |
|
| 640 |
}
|
| 641 |
|
| 642 |
if(iVtx ==0 && print && false) {
|
| 643 |
std::cout<<" Trk passes cuts: "<<std::endl;
|
| 644 |
std::cout<<" px tot = "<<totTrkMom.Px()<<std::endl;
|
| 645 |
std::cout<<" py tot = "<<totTrkMom.Py()<<std::endl;
|
| 646 |
}
|
| 647 |
|
| 648 |
double ptvtx = totTrkMom.Pt();
|
| 649 |
if(iVtx ==0 && print && false)
|
| 650 |
std::cout<<" Total TkPt = "<<ptvtx<<std::endl;
|
| 651 |
double pthiggs = higgsMom.Pt();
|
| 652 |
if(iVtx ==0 && print && false)
|
| 653 |
std::cout<<" Total H Pt = "<<pthiggs<<std::endl;
|
| 654 |
if(iVtx == 0) ptgg = pthiggs;
|
| 655 |
ptbal [iVtx] = (totTrkMom.Px()*(newMomFst.Px()+newMomSec.Px()));
|
| 656 |
ptbal [iVtx] += (totTrkMom.Py()*(newMomFst.Py()+newMomSec.Py()));
|
| 657 |
ptbal [iVtx] = -ptbal[iVtx]/pthiggs;
|
| 658 |
ptasym[iVtx] = (ptvtx - pthiggs)/(ptvtx + pthiggs);
|
| 659 |
|
| 660 |
if(iVtx ==0 && print && false) {
|
| 661 |
std::cout<<" Results: ptbal = "<<ptbal [iVtx]<<std::endl;
|
| 662 |
std::cout<<" ptasym = "<<ptasym[iVtx]<<std::endl;
|
| 663 |
}
|
| 664 |
|
| 665 |
for(unsigned int cVtx =0; cVtx < iVtx; ++cVtx) {
|
| 666 |
if(ptbal [iVtx] > ptbal [cVtx])
|
| 667 |
ptbal_rank[cVtx]++;
|
| 668 |
else
|
| 669 |
ptbal_rank[iVtx]++;
|
| 670 |
if(ptasym [iVtx] > ptasym [cVtx])
|
| 671 |
ptasym_rank[cVtx]++;
|
| 672 |
else
|
| 673 |
ptasym_rank[iVtx]++;
|
| 674 |
}
|
| 675 |
}
|
| 676 |
|
| 677 |
|
| 678 |
// compute the total rank
|
| 679 |
for(unsigned int iVtx = 0; iVtx < numVertices; ++iVtx) {
|
| 680 |
if(print && false) {
|
| 681 |
std::cout<<" Vertex #"<<iVtx<<" has rank PTB "<<ptbal_rank[iVtx]<<" ("<<ptbal[iVtx]<<")"<<std::endl;
|
| 682 |
std::cout<<" Vertex #"<<iVtx<<" has rank PTSYM "<<ptasym_rank[iVtx]<<" ("<<ptasym[iVtx]<<")"<<std::endl;
|
| 683 |
}
|
| 684 |
ptasym_rank [iVtx] = ptbal_rank [iVtx]*ptasym_rank [iVtx]*(iVtx+1);
|
| 685 |
total_rank [iVtx] = 0;
|
| 686 |
for(unsigned int cVtx =0; cVtx < iVtx; ++cVtx) {
|
| 687 |
if(ptasym_rank [iVtx] > ptasym_rank [cVtx])
|
| 688 |
total_rank[iVtx]++;
|
| 689 |
else if(ptasym_rank [iVtx] == ptasym_rank [cVtx]) {
|
| 690 |
if(ptbal_rank [iVtx] > ptbal_rank [cVtx])
|
| 691 |
total_rank[iVtx]++;
|
| 692 |
else
|
| 693 |
total_rank[cVtx]++;
|
| 694 |
}
|
| 695 |
else
|
| 696 |
total_rank[cVtx]++;
|
| 697 |
}
|
| 698 |
}
|
| 699 |
|
| 700 |
if(print && false) std::cout<<std::endl;
|
| 701 |
|
| 702 |
unsigned int bestIdx = 0;
|
| 703 |
for(unsigned int iVtx = 0; iVtx < numVertices; ++iVtx) {
|
| 704 |
if(total_rank[iVtx] == 0) bestIdx = iVtx;
|
| 705 |
if(print && false) {
|
| 706 |
std::cout<<" Vertex #"<<iVtx<<" has rank "<<total_rank[iVtx]<<std::endl;
|
| 707 |
}
|
| 708 |
}
|
| 709 |
|
| 710 |
delete[] ptbal_rank ;
|
| 711 |
delete[] ptasym_rank ;
|
| 712 |
delete[] ptbal ;
|
| 713 |
delete[] ptasym ;
|
| 714 |
|
| 715 |
//return bestIdx;
|
| 716 |
|
| 717 |
// check if there's a conversion among the pre-selected photons
|
| 718 |
const DecayParticle* conv1 = PhotonTools::MatchedCiCConversion(ph1, fConversions, 0.1, 0.1, 0.1, print);
|
| 719 |
const DecayParticle* conv2 = PhotonTools::MatchedCiCConversion(ph2, fConversions, 0.1, 0.1, 0.1, print);
|
| 720 |
|
| 721 |
if(print && false) {
|
| 722 |
if (conv1) {
|
| 723 |
std::cout<<" Photon 1 has has conversion with P = "<<conv1->Prob()<<std::endl;
|
| 724 |
std::cout<<" Rho = "<<conv1->Position().Rho()<<std::endl;
|
| 725 |
std::cout<<" Z = "<<conv1->Position().Z()<<std::endl;
|
| 726 |
}
|
| 727 |
if (conv2) {
|
| 728 |
std::cout<<" Photon 2 has has conversion with P = "<<conv2->Prob()<<std::endl;
|
| 729 |
std::cout<<" Rho = "<<conv2->Position().Rho()<<std::endl;
|
| 730 |
std::cout<<" Z = "<<conv2->Position().Z()<<std::endl;
|
| 731 |
}
|
| 732 |
}
|
| 733 |
|
| 734 |
if( conv1 && ( conv1->Prob() < 0.0005) ) conv1 = NULL;
|
| 735 |
if( conv2 && ( conv2->Prob() < 0.0005) ) conv2 = NULL;
|
| 736 |
|
| 737 |
double zconv = 0.;
|
| 738 |
double dzconv = 0.;
|
| 739 |
|
| 740 |
if(conv1 || conv2) {
|
| 741 |
if( !conv2 ){
|
| 742 |
//const mithep::ThreeVector caloPos1(ph1->SCluster()->Point());
|
| 743 |
const mithep::ThreeVector caloPos1(ph1->CaloPos());
|
| 744 |
zconv = conv1->Z0EcalVtx(bsp->Position(), caloPos1);
|
| 745 |
if( ph1->IsEB() ) {
|
| 746 |
double rho = conv1->Position().Rho();
|
| 747 |
if ( rho < 15. ) dzconv = 0.06;
|
| 748 |
else if( rho < 60. ) dzconv = 0.67;
|
| 749 |
else dzconv = 2.04;
|
| 750 |
} else {
|
| 751 |
double z = conv1->Position().Z();
|
| 752 |
if ( TMath::Abs(z) < 50. ) dzconv = 0.18;
|
| 753 |
else if( TMath::Abs(z) < 100.) dzconv = 0.61;
|
| 754 |
else dzconv = 0.99;
|
| 755 |
}
|
| 756 |
} else if( !conv1 ) {
|
| 757 |
//const mithep::ThreeVector caloPos2(ph2->SCluster()->Point());
|
| 758 |
const mithep::ThreeVector caloPos2(ph2->CaloPos());
|
| 759 |
zconv = conv2->Z0EcalVtx(bsp->Position(), caloPos2);
|
| 760 |
if( ph2->IsEB() ) {
|
| 761 |
double rho = conv2->Position().Rho();
|
| 762 |
if ( rho < 15. ) dzconv = 0.06;
|
| 763 |
else if( rho < 60. ) dzconv = 0.67;
|
| 764 |
else dzconv = 2.04;
|
| 765 |
} else {
|
| 766 |
double z = conv2->Position().Z();
|
| 767 |
if ( TMath::Abs(z) < 50. ) dzconv = 0.18;
|
| 768 |
else if( TMath::Abs(z) < 100.) dzconv = 0.61;
|
| 769 |
else dzconv = 0.99;
|
| 770 |
}
|
| 771 |
} else {
|
| 772 |
//const mithep::ThreeVector caloPos1(ph1->SCluster()->Point());
|
| 773 |
const mithep::ThreeVector caloPos1(ph1->CaloPos());
|
| 774 |
double z1 = conv1->Z0EcalVtx(bsp->Position(), caloPos1);
|
| 775 |
double dz1 = 0.;
|
| 776 |
if( ph1->IsEB() ) {
|
| 777 |
double rho = conv1->Position().Rho();
|
| 778 |
if ( rho < 15. ) dz1 = 0.06;
|
| 779 |
else if( rho < 60. ) dz1 = 0.67;
|
| 780 |
else dz1 = 2.04;
|
| 781 |
} else {
|
| 782 |
double z = conv1->Position().Z();
|
| 783 |
if ( TMath::Abs(z) < 50. ) dz1 = 0.18;
|
| 784 |
else if( TMath::Abs(z) < 100.) dz1 = 0.61;
|
| 785 |
else dz1 = 0.99;
|
| 786 |
}
|
| 787 |
//const mithep::ThreeVector caloPos2(ph2->SCluster()->Point());
|
| 788 |
const mithep::ThreeVector caloPos2(ph2->CaloPos());
|
| 789 |
double z2 = conv2->Z0EcalVtx(bsp->Position(), caloPos2);
|
| 790 |
double dz2 = 0.;
|
| 791 |
if( ph2->IsEB() ) {
|
| 792 |
double rho = conv2->Position().Rho();
|
| 793 |
if ( rho < 15. ) dz2 = 0.06;
|
| 794 |
else if( rho < 60. ) dz2 = 0.67;
|
| 795 |
else dz2 = 2.04;
|
| 796 |
} else {
|
| 797 |
double z = conv2->Position().Z();
|
| 798 |
if ( TMath::Abs(z) < 50. ) dz2 = 0.18;
|
| 799 |
else if( TMath::Abs(z) < 100.) dz2 = 0.61;
|
| 800 |
else dz2 = 0.99;
|
| 801 |
}
|
| 802 |
|
| 803 |
if(print) {
|
| 804 |
std::cout<<" z1 = "<<z1<<std::endl;
|
| 805 |
std::cout<<" dz1 = "<<dz1<<std::endl;
|
| 806 |
std::cout<<" z2 = "<<z2<<std::endl;
|
| 807 |
std::cout<<" dz2 = "<<dz2<<std::endl;
|
| 808 |
}
|
| 809 |
|
| 810 |
zconv = ( 1./(1./dz1/dz1 + 1./dz2/dz2 )*(z1/dz1/dz1 + z2/dz2/dz2) ) ; // weighted average
|
| 811 |
dzconv = TMath::Sqrt( 1./(1./dz1/dz1 + 1./dz2/dz2)) ;
|
| 812 |
}
|
| 813 |
|
| 814 |
if(print) {
|
| 815 |
std::cout<<" Conversion Z = "<<zconv<<std::endl;
|
| 816 |
std::cout<<" dZ = "<<dzconv<<std::endl;
|
| 817 |
}
|
| 818 |
|
| 819 |
|
| 820 |
if(true) {
|
| 821 |
|
| 822 |
// loop over all ranked Vertices and choose the closest to the Conversion one
|
| 823 |
int maxVertices = ( ptgg > 30 ? 3 : 5);
|
| 824 |
double minDz = -1.;
|
| 825 |
|
| 826 |
if(print) std::cout<<std::endl<<" looping over vertices... "<<ptgg<<" "<<maxVertices<<std::endl;
|
| 827 |
|
| 828 |
for(unsigned int iVtx =0; iVtx < numVertices; ++iVtx) {
|
| 829 |
|
| 830 |
if(print) std::cout<<" "<<iVtx<<" has rank "<<total_rank[iVtx]<<std::endl;
|
| 831 |
|
| 832 |
if(total_rank[iVtx] < maxVertices) {
|
| 833 |
const Vertex* tVtx = fPV->At(iVtx);
|
| 834 |
double tDz = TMath::Abs(zconv - tVtx->Z());
|
| 835 |
if(print) std::cout<<" is considered with tDz = "<<tDz<<std::endl;
|
| 836 |
if( (minDz < 0. || tDz < minDz) && ( tDz < dzconv ) ) {
|
| 837 |
minDz = tDz;
|
| 838 |
bestIdx = iVtx;
|
| 839 |
if(print) std::cout<<" and is the best now."<<std::endl;
|
| 840 |
}
|
| 841 |
}
|
| 842 |
}
|
| 843 |
} else {
|
| 844 |
unsigned int bestIdxTmp = bestIdx;
|
| 845 |
|
| 846 |
// loop over all ranked Vertices and choose the closest to the Conversion one
|
| 847 |
double minDz = -1.;
|
| 848 |
int maxVertices = ( ptgg > 30 ? 3 : 5);
|
| 849 |
|
| 850 |
if(print) std::cout<<std::endl<<" looping over vertices... "<<ptgg<<" "<<maxVertices<<std::endl;
|
| 851 |
|
| 852 |
for(unsigned int iVtx =0; iVtx < numVertices; ++iVtx) {
|
| 853 |
|
| 854 |
if(print) std::cout<<" "<<iVtx<<" has rank "<<total_rank[iVtx]<<std::endl;
|
| 855 |
|
| 856 |
const Vertex* tVtx = fPV->At(iVtx);
|
| 857 |
double tDz = TMath::Abs(zconv - tVtx->Z());
|
| 858 |
if(print) std::cout<<" is considered with tDz = "<<tDz<<std::endl;
|
| 859 |
if( (minDz < 0. || tDz < minDz) && ( tDz < dzconv ) ) {
|
| 860 |
minDz = tDz;
|
| 861 |
bestIdxTmp = iVtx;
|
| 862 |
if(print) std::cout<<" and is the best now."<<std::endl;
|
| 863 |
}
|
| 864 |
}
|
| 865 |
|
| 866 |
// check if best Vtx is among higest ranked ones
|
| 867 |
if(total_rank[bestIdxTmp] < maxVertices)
|
| 868 |
bestIdx = bestIdxTmp;
|
| 869 |
}
|
| 870 |
}
|
| 871 |
|
| 872 |
delete[] total_rank ;
|
| 873 |
return bestIdx;
|
| 874 |
}
|
| 875 |
|
| 876 |
void PhotonCiCMod::findHiggsPtAndZ(Float_t& pt, Float_t& decayZ) {
|
| 877 |
|
| 878 |
pt = -999.;
|
| 879 |
decayZ = -999.;
|
| 880 |
|
| 881 |
// loop over all GEN particles and look for status 1 photons
|
| 882 |
for(UInt_t i=0; i<fMCParticles->GetEntries(); ++i) {
|
| 883 |
const MCParticle* p = fMCParticles->At(i);
|
| 884 |
if( !(p->Is(MCParticle::kH)) ) continue;
|
| 885 |
pt=p->Pt();
|
| 886 |
decayZ = p->DecayVertex().Z();
|
| 887 |
break;
|
| 888 |
}
|
| 889 |
|
| 890 |
return;
|
| 891 |
}
|
