| 1 |
+ |
#include <sstream> |
| 2 |
|
#include "MitPhysics/Mods/interface/PhotonTreeWriter.h" |
| 3 |
|
#include "MitAna/DataTree/interface/PhotonCol.h" |
| 4 |
|
#include "MitAna/DataTree/interface/PFCandidateCol.h" |
| 107 |
|
fFillClusterArrays (kFALSE), |
| 108 |
|
fFillVertexTree (kFALSE), |
| 109 |
|
fDo2012LepTag (kFALSE), |
| 110 |
+ |
fVerbosityLevel (0), |
| 111 |
|
fPhFixDataFile (gSystem->Getenv("CMSSW_BASE") + |
| 112 |
|
TString("/src/MitPhysics/data/PhotonFixSTART42V13.dat")), |
| 113 |
|
fBeamspotWidth (5.8), |
| 114 |
+ |
fTmpFile (0), |
| 115 |
+ |
|
| 116 |
+ |
// JV: moved up the initializtion of fTupleName to avoid compilation warning |
| 117 |
+ |
fTupleName ("hPhotonTree"), |
| 118 |
|
|
| 119 |
|
fElectronIDMVA(0), |
| 120 |
|
fElectronMVAWeights_Subdet0Pt10To20(""), |
| 123 |
|
fElectronMVAWeights_Subdet0Pt20ToInf(""), |
| 124 |
|
fElectronMVAWeights_Subdet1Pt20ToInf(""), |
| 125 |
|
fElectronMVAWeights_Subdet2Pt20ToInf(""), |
| 120 |
– |
fTheRhoType(RhoUtilities::DEFAULT), |
| 126 |
|
|
| 127 |
< |
fTupleName ("hPhotonTree") |
| 127 |
> |
fTheRhoType(RhoUtilities::DEFAULT) |
| 128 |
> |
|
| 129 |
|
{ |
| 130 |
|
// Constructor |
| 131 |
|
} |
| 133 |
|
PhotonTreeWriter::~PhotonTreeWriter() |
| 134 |
|
{ |
| 135 |
|
// Destructor |
| 136 |
+ |
// Deal with the temporary file here? |
| 137 |
+ |
// fTmpFile->Write(); |
| 138 |
+ |
fTmpFile->Close(); |
| 139 |
+ |
TString shellcmd = TString("rm ") + TString(fTmpFile->GetName()); |
| 140 |
+ |
delete fTmpFile; |
| 141 |
+ |
cout << shellcmd.Data() << endl; |
| 142 |
+ |
gSystem->Exec(shellcmd.Data()); |
| 143 |
+ |
|
| 144 |
|
} |
| 145 |
|
|
| 146 |
|
//-------------------------------------------------------------------------------------------------- |
| 1075 |
|
fDiphotonEvent->vbfTag = 1; |
| 1076 |
|
} |
| 1077 |
|
} |
| 1078 |
< |
} |
| 1079 |
< |
|
| 1078 |
> |
} // End of VBF tag stuff |
| 1079 |
> |
|
| 1080 |
|
// ttH tag stuff |
| 1081 |
|
fDiphotonEvent->tthTag = -1; |
| 1082 |
< |
if (fApplyTTHTag) { |
| 1083 |
< |
ApplyTTHTag(phHard,phSoft,selvtx); |
| 1082 |
> |
if (fApplyTTHTag && phHard && phSoft && selvtx) { |
| 1083 |
> |
ApplyTTHTag(phHard, phSoft, selvtx); |
| 1084 |
> |
PrintTTHDebugInfo(); |
| 1085 |
|
} |
| 1086 |
|
|
| 1087 |
|
//printf("vbfbdt:%f\n",fDiphotonEvent->vbfbdt); |
| 1257 |
|
fDiphotonEvent = new PhotonTreeWriterDiphotonEvent; |
| 1258 |
|
fSinglePhoton = new PhotonTreeWriterPhoton<16>; |
| 1259 |
|
|
| 1260 |
< |
TFile *ftmp = TFile::Open(TString::Format("%s_tmp.root",GetName()),"RECREATE"); |
| 1260 |
> |
fTmpFile = TFile::Open(TString::Format("%s_tmp.root",GetName()),"RECREATE"); |
| 1261 |
|
|
| 1262 |
|
if (fWriteDiphotonTree) { |
| 1263 |
|
hCiCTuple = new TTree(fTupleName.Data(),fTupleName.Data()); |
| 1510 |
|
// ----------------------------------------------------- |
| 1511 |
|
// PF-CiC4 Debug Stuff |
| 1512 |
|
std::vector<double> debugVals; |
| 1513 |
< |
bool tmpPass = PhotonTools::PassCiCPFIsoSelection(p, vtx, fPFCands, vtxCol, rho, 20., &debugVals); |
| 1513 |
> |
// bool tmpPass = PhotonTools::PassCiCPFIsoSelection(p, vtx, fPFCands, vtxCol, rho, 20., &debugVals); |
| 1514 |
> |
PhotonTools::PassCiCPFIsoSelection(p, vtx, fPFCands, vtxCol, rho, 20., &debugVals); |
| 1515 |
|
if( debugVals.size() == 13 ) { |
| 1516 |
|
pfcic4_tIso1 = debugVals[0]; |
| 1517 |
|
pfcic4_tIso2 = debugVals[1]; |
| 2260 |
|
|
| 2261 |
|
} |
| 2262 |
|
|
| 2263 |
+ |
|
| 2264 |
|
//_____________________________________________________________________________ |
| 2265 |
< |
void PhotonTreeWriter::ApplyTTHTag(const Photon *phHard, const Photon *phSoft,const Vertex *selvtx) |
| 2265 |
> |
// Applies the ttH tag given precelected leading and trailing photons |
| 2266 |
> |
// phHard and phSoft and the corresponding (pre?) selected vertex selvtx. |
| 2267 |
> |
// The result is stored as an integer value of the tthTag variable |
| 2268 |
> |
// entering the diphoton event record. |
| 2269 |
> |
void PhotonTreeWriter::ApplyTTHTag(const Photon *phHard, |
| 2270 |
> |
const Photon *phSoft, |
| 2271 |
> |
const Vertex *selvtx) |
| 2272 |
|
{ |
| 2273 |
|
// ttH tag = -1 .. not applied |
| 2274 |
|
// 0 .. not tagged |
| 2276 |
|
// 2 .. tagged as a hadronic ttH event |
| 2277 |
|
fDiphotonEvent->tthTag = 0; |
| 2278 |
|
|
| 2279 |
< |
// Selection taken from the AN2012_480_V6 of 24 April 2013, further |
| 2280 |
< |
// refferred to as "the AN" |
| 2279 |
> |
// The selection taken from the AN2012_480_V6 of 24 April 2013, further |
| 2280 |
> |
// referred to as "the AN", and the approval slides at |
| 2281 |
> |
// https://twiki.cern.ch/twiki/pub/CMS/JanVeverkaHgg/micheli_ttH_approval_20130508_2_reduced.pdf, |
| 2282 |
> |
// further referred to as "the slides" or "the slide n" (n = 1, 2, ..) |
| 2283 |
|
|
| 2284 |
|
// Check the pt of the photons, see L141 and L142 of the AN |
| 2285 |
< |
if (phHard->Pt() < 33.) return; |
| 2286 |
< |
if (phSoft->Pt() < 25.) return; |
| 2285 |
> |
if (fDiphotonEvent->photons[0].Pt() < 33.) return; |
| 2286 |
> |
if (fDiphotonEvent->photons[1].Pt() < 25.) return; |
| 2287 |
|
|
| 2288 |
|
// Init final-state object counters |
| 2289 |
|
UInt_t nJets = 0; |
| 2291 |
|
UInt_t nElectrons = 0; |
| 2292 |
|
UInt_t nMuons = 0; |
| 2293 |
|
|
| 2294 |
+ |
// Get the selected vertex |
| 2295 |
+ |
|
| 2296 |
|
// Loop over jets, count those passing selection. |
| 2297 |
|
// No Delta R(gamma, j) requirement!? |
| 2298 |
|
for(UInt_t ijet=0; ijet < fPFJets->GetEntries(); ++ijet){ |
| 2307 |
|
// (Aram & Valentina say so) |
| 2308 |
|
if (!JetTools::passPFLooseId(pfjet)) continue; |
| 2309 |
|
// Apply the jet ID as given in Table 4 |
| 2310 |
< |
// TODO: JetTools::betaStar or JetTools::betaStarClassic? |
| 2310 |
> |
// TODO: JetTools::betaStar or JetTools::betaStarClassic? (yes, classic) |
| 2311 |
|
// Go for the latter now as it is used in the JetIDMVA::passCut(..) which in turn |
| 2312 |
|
// is used for the VBF tag. |
| 2313 |
|
Double_t betaStar = JetTools::betaStarClassic(pfjet, selvtx, fPV); |
| 2314 |
|
// TODO: Really 0.67? JetIDMVA::passCut(..) and AN-13-008 use 0.64 insted. |
| 2315 |
|
if (betaStar > 0.2 * log(fPV->GetEntries() - 0.67)) continue; |
| 2316 |
< |
// TODO: is the second argument iPFType = -1 correct? |
| 2316 |
> |
// TODO: is the second argument iPFType = -1 correct? yes |
| 2317 |
|
if (JetTools::dR2Mean(pfjet, -1) > 0.065) continue; |
| 2318 |
|
// this jet passes, count it in |
| 2319 |
|
++nJets; |
| 2351 |
|
for (UInt_t iele=0; iele < fLeptonTagElectrons->GetEntries(); ++iele) { |
| 2352 |
|
const Electron *ele = fLeptonTagElectrons->At(iele); |
| 2353 |
|
// Apply kinematic cuts, see L133 and L134 of the AN |
| 2354 |
< |
if (ele->Pt() < 20. || ele->AbsEta() < 2.5) continue; |
| 2354 |
> |
if (ele->Pt() < 20. || ele->AbsEta() > 2.5) continue; |
| 2355 |
> |
// Require separation between this electron and both photons, |
| 2356 |
> |
// see the slide 7, bullet 5 |
| 2357 |
> |
if (MathUtils::DeltaR(ele, phHard) < 1.0) continue; |
| 2358 |
> |
if (MathUtils::DeltaR(ele, phSoft) < 1.0) continue; |
| 2359 |
|
++nElectrons; |
| 2360 |
|
} |
| 2361 |
|
|
| 2363 |
|
for (UInt_t imu=0; imu < fLeptonTagMuons->GetEntries(); ++imu) { |
| 2364 |
|
const Muon *mu = fLeptonTagMuons->At(imu); |
| 2365 |
|
// Apply kinematic cuts, see L132 and L134 of the AN |
| 2366 |
< |
if (mu->Pt() < 20. || mu->AbsEta() < 2.4) continue; |
| 2366 |
> |
if (mu->Pt() < 20. || mu->AbsEta() > 2.4) continue; |
| 2367 |
> |
// Require separation between this electron and both photons, |
| 2368 |
> |
// see the slide 7, bullet 2 |
| 2369 |
> |
// TODO: Check with authors that this is indeed applied. |
| 2370 |
> |
if (MathUtils::DeltaR(mu, phHard) < 1.0) continue; |
| 2371 |
> |
if (MathUtils::DeltaR(mu, phSoft) < 1.0) continue; |
| 2372 |
|
++nMuons; |
| 2373 |
|
} |
| 2374 |
|
|
| 2390 |
|
|
| 2391 |
|
} // void PhotonTreeWriter::ApplyTTHTag() |
| 2392 |
|
|
| 2393 |
+ |
|
| 2394 |
+ |
//_____________________________________________________________________________ |
| 2395 |
+ |
void PhotonTreeWriter::Terminate() |
| 2396 |
+ |
{ |
| 2397 |
+ |
// Run finishing code on the computer (slave) that did the analysis |
| 2398 |
+ |
} |
| 2399 |
+ |
|
| 2400 |
+ |
|
| 2401 |
+ |
//_____________________________________________________________________________ |
| 2402 |
+ |
void PhotonTreeWriter::PrintTTHDebugInfo() |
| 2403 |
+ |
{ |
| 2404 |
+ |
if (fVerbosityLevel > 0) { |
| 2405 |
+ |
cout << "JV: Processing run " << fDiphotonEvent->run |
| 2406 |
+ |
<< " event " << fDiphotonEvent->evt << " .. " << endl; |
| 2407 |
+ |
PrintTTHDecay(); |
| 2408 |
+ |
} |
| 2409 |
+ |
|
| 2410 |
+ |
if (fVerbosityLevel > 1) { |
| 2411 |
+ |
PrintGenElectrons(); |
| 2412 |
+ |
PrintElectrons("reco electrons", fElectrons); |
| 2413 |
+ |
PrintElectrons("lepton-tag electrons", fLeptonTagElectrons); |
| 2414 |
+ |
} |
| 2415 |
+ |
} // void PhotonTreeWriter::PrintTTHDebugInfo() |
| 2416 |
+ |
|
| 2417 |
+ |
|
| 2418 |
+ |
//_____________________________________________________________________________ |
| 2419 |
+ |
void PhotonTreeWriter::PrintTTHDecay() |
| 2420 |
+ |
{ |
| 2421 |
+ |
// Initialize W decay type counters |
| 2422 |
+ |
UInt_t nLeptonic = 0; |
| 2423 |
+ |
UInt_t nHadronic = 0; |
| 2424 |
+ |
Int_t WplusDau1PID = 0; |
| 2425 |
+ |
Int_t WminusDau1PID = 0; |
| 2426 |
+ |
// loop over all GEN particles and look for W bosons |
| 2427 |
+ |
for(UInt_t i=0; i<fMCParticles->GetEntries(); ++i) { |
| 2428 |
+ |
const MCParticle* p = fMCParticles->At(i); |
| 2429 |
+ |
if (p->Is(MCParticle::kW) && p->Status() == 3 && p->NDaughters() > 0) { |
| 2430 |
+ |
// found a W in the ME |
| 2431 |
+ |
cout << "JV: " << p->PdgEntry()->GetName() |
| 2432 |
+ |
<< "[" << i << "] -> "; |
| 2433 |
+ |
|
| 2434 |
+ |
// Set the PID of the first daughter |
| 2435 |
+ |
if (p->PdgId() > 0) { |
| 2436 |
+ |
WplusDau1PID = p->Daughter(0)->PdgId(); |
| 2437 |
+ |
} else { |
| 2438 |
+ |
WminusDau1PID = p->Daughter(0)->PdgId(); |
| 2439 |
+ |
} |
| 2440 |
+ |
|
| 2441 |
+ |
// loop over W daughters |
| 2442 |
+ |
for (UInt_t j=0; j < p->NDaughters(); ++j) { |
| 2443 |
+ |
const MCParticle *d = p->Daughter(j); |
| 2444 |
+ |
if (d->Status() == 2) continue; |
| 2445 |
+ |
cout << d->PdgEntry()->GetName() << " "; |
| 2446 |
+ |
} // loop over W daughters |
| 2447 |
+ |
cout << endl; |
| 2448 |
+ |
|
| 2449 |
+ |
// Check the decay type |
| 2450 |
+ |
if (p->NDaughters() > 0 && p->Daughter(0)->IsQuark()) { |
| 2451 |
+ |
++nHadronic; |
| 2452 |
+ |
} else { |
| 2453 |
+ |
++nLeptonic; |
| 2454 |
+ |
} |
| 2455 |
+ |
} // found the W |
| 2456 |
+ |
} // loop over all GEN particles |
| 2457 |
+ |
|
| 2458 |
+ |
cout << "JV: ttbar decay is "; |
| 2459 |
+ |
switch (nHadronic) { |
| 2460 |
+ |
case 0: cout << "leptonic"; break; |
| 2461 |
+ |
case 1: cout << "semi-leptonic"; break; |
| 2462 |
+ |
case 2: cout << "hadronic"; break; |
| 2463 |
+ |
default: cout << "ERROR"; |
| 2464 |
+ |
} |
| 2465 |
+ |
cout << endl; |
| 2466 |
+ |
|
| 2467 |
+ |
cout << "JV: tthTag WplusDau1PID WminusDau1PID:\t" |
| 2468 |
+ |
<< fDiphotonEvent->tthTag << "\t" |
| 2469 |
+ |
<< WplusDau1PID << "\t" |
| 2470 |
+ |
<< WminusDau1PID << endl; |
| 2471 |
+ |
return; |
| 2472 |
+ |
} // void PrintTTHDecay() |
| 2473 |
+ |
|
| 2474 |
+ |
|
| 2475 |
+ |
//_____________________________________________________________________________ |
| 2476 |
+ |
void PhotonTreeWriter::PrintGenElectrons() |
| 2477 |
+ |
{ |
| 2478 |
+ |
ostringstream report; |
| 2479 |
+ |
report << "JV: gen electrons" << endl; |
| 2480 |
+ |
UInt_t nEl = 0; |
| 2481 |
+ |
for (UInt_t i=0; i < fMCParticles->GetEntries(); ++i){ |
| 2482 |
+ |
const MCParticle *p = fMCParticles->At(i); |
| 2483 |
+ |
if (p->Is(MCParticle::kEl) && |
| 2484 |
+ |
p->Status() == 1 && |
| 2485 |
+ |
p->Pt() > 15 && |
| 2486 |
+ |
p->AbsEta() < 2.6) { |
| 2487 |
+ |
++nEl; |
| 2488 |
+ |
report << " " << p->PdgEntry()->GetName() << ": " |
| 2489 |
+ |
<< "pt=" << p->Pt() << ", " |
| 2490 |
+ |
<< "eta=" << p->AbsEta() << endl; |
| 2491 |
+ |
} // Electron preselection |
| 2492 |
+ |
} // End of loop over MC particles |
| 2493 |
+ |
if (nEl > 0) cout << report.str(); |
| 2494 |
+ |
} // void PhotonTreeWriter::PrintGenElectrons() |
| 2495 |
+ |
|
| 2496 |
+ |
|
| 2497 |
+ |
//_____________________________________________________________________________ |
| 2498 |
+ |
void PhotonTreeWriter::PrintElectrons(const char *tag, |
| 2499 |
+ |
const ElectronCol *electrons) |
| 2500 |
+ |
{ |
| 2501 |
+ |
ostringstream report; |
| 2502 |
+ |
report << "JV: " << tag << endl; |
| 2503 |
+ |
UInt_t nEl = 0; |
| 2504 |
+ |
for (UInt_t i=0; i < electrons->GetEntries(); ++i){ |
| 2505 |
+ |
const Electron *e = electrons->At(i); |
| 2506 |
+ |
if (e->Pt() > 20 && e->AbsEta() < 2.5) { |
| 2507 |
+ |
++nEl; |
| 2508 |
+ |
report << " q=" << e->Charge() << ", " |
| 2509 |
+ |
<< "pt=" << e->Pt() << ", " |
| 2510 |
+ |
<< "eta=" << e->AbsEta() << endl; |
| 2511 |
+ |
} // Electron preselection |
| 2512 |
+ |
} // End of loop over reco electrons |
| 2513 |
+ |
if (nEl > 0) cout << report.str(); |
| 2514 |
+ |
} // void PhotonTreeWriter::PrintRecoElectrons() |
| 2515 |
+ |
|
| 2516 |
+ |
|
