--- UserCode/kiesel/TreeWriter/treeWriter.cc 2013/04/10 20:46:33 1.8 +++ UserCode/kiesel/TreeWriter/treeWriter.cc 2013/04/17 09:57:52 1.14 @@ -8,33 +8,35 @@ using namespace std; -TreeWriter::TreeWriter(TString inputName, TString outputName, int loggingVerbosity_ ) { +TreeWriter::TreeWriter( TString inputName, TString outputName, int loggingVerbosity_ ) { // read the input file inputTree = new TChain("susyTree"); if (loggingVerbosity_ > 0) std::cout << "Add files to chain" << std::endl; inputTree->Add( inputName ); + Init( outputName, loggingVerbosity_ ); +} + +TreeWriter::TreeWriter( TChain* inputTree_, TString outputName, int loggingVerbosity_ ) { + inputTree = inputTree_; + Init( outputName, loggingVerbosity_ ); +} + + +void TreeWriter::Init( TString outputName, int loggingVerbosity_ ) { if (loggingVerbosity_ > 0) std::cout << "Set Branch Address of susy::Event" << std::endl; event = new susy::Event; inputTree->SetBranchAddress("susyEvent", &event); - if (outputName == "") { - // Here the output filename is generated automaticly - try{ - string fName = (string) inputName; - // Search for 'susyEvents_', and get the unique characters afterwards. - // eg /path/susyEvents_123_Myl.root -> susyTree_123_Mly.root - outputName = "susyTree_" + fName.substr(fName.find("susyEvents_")+11,-1); - } catch (int e ) { - outputName = "susyTree.root"; - } - } + // Here the number of proceeded events will be stored. For plotting, simply use L*sigma/eventNumber + eventNumbers = new TH1F("eventNumbers", "Histogram containing number of generated Events", 1, 0, 1); + eventNumbers->GetXaxis()->SetBinLabel(1,"Number of generated Events"); + // open the output file if (loggingVerbosity_>0) std::cout << "Open file " << outputName << " for writing." << std::endl; - // open the output file outFile = new TFile( outputName, "recreate" ); tree = new TTree("susyTree","Tree for single photon analysis"); @@ -43,12 +45,18 @@ TreeWriter::TreeWriter(TString inputName reportEvery = 1000; loggingVerbosity = loggingVerbosity_; skim = true; + pileupHisto = 0; +} +void TreeWriter::PileUpWeightFile( string pileupFileName ) { + TFile *puFile = new TFile( pileupFileName.c_str() ); + pileupHisto = (TH1F*) puFile->Get("pileup"); } TreeWriter::~TreeWriter() { - if (!inputTree) return; - delete inputTree->GetCurrentFile(); + if (pileupHisto != 0 ) + delete pileupHisto; + inputTree->GetCurrentFile()->Close(); } // useful functions @@ -56,63 +64,92 @@ float TreeWriter::deltaR( TLorentzVector return sqrt(pow(v1.Eta() - v2.Eta(), 2) + pow(v1.Phi() - v2.Phi(), 2) ); } +float effectiveAreaElectron( float eta ) { + // see https://twiki.cern.ch/twiki/bin/view/CMS/EgammaEARhoCorrection + // only for Delta R = 0.3 on 2012 Data + eta = fabs( eta ); + float ea; + if( eta < 1.0 ) ea = 0.13; + else if( eta < 1.479 ) ea = 0.14; + else if( eta < 2.0 ) ea = 0.07; + else if( eta < 2.2 ) ea = 0.09; + else if( eta < 2.3 ) ea = 0.11; + else if( eta < 2.4 ) ea = 0.11; + else ea = 0.14; + return ea; +} + // correct iso, see https://twiki.cern.ch/twiki/bin/view/CMS/CutBasedPhotonID2012 float chargedHadronIso_corrected(susy::Photon gamma, float rho) { - float eta = fabs(gamma.caloPosition.Eta()); - float ea; + float eta = fabs(gamma.caloPosition.Eta()); + float ea; - if(eta < 1.0) ea = 0.012; - else if(eta < 1.479) ea = 0.010; - else if(eta < 2.0) ea = 0.014; - else if(eta < 2.2) ea = 0.012; - else if(eta < 2.3) ea = 0.016; - else if(eta < 2.4) ea = 0.020; - else ea = 0.012; + if(eta < 1.0) ea = 0.012; + else if(eta < 1.479) ea = 0.010; + else if(eta < 2.0) ea = 0.014; + else if(eta < 2.2) ea = 0.012; + else if(eta < 2.3) ea = 0.016; + else if(eta < 2.4) ea = 0.020; + else ea = 0.012; - float iso = gamma.chargedHadronIso; - iso = max(iso - rho*ea, (float)0.); + float iso = gamma.chargedHadronIso; + iso = max(iso - rho*ea, (float)0.); - return iso; + return iso; } float neutralHadronIso_corrected(susy::Photon gamma, float rho) { - float eta = fabs(gamma.caloPosition.Eta()); - float ea; + float eta = fabs(gamma.caloPosition.Eta()); + float ea; - if(eta < 1.0) ea = 0.030; - else if(eta < 1.479) ea = 0.057; - else if(eta < 2.0) ea = 0.039; - else if(eta < 2.2) ea = 0.015; - else if(eta < 2.3) ea = 0.024; - else if(eta < 2.4) ea = 0.039; - else ea = 0.072; + if(eta < 1.0) ea = 0.030; + else if(eta < 1.479) ea = 0.057; + else if(eta < 2.0) ea = 0.039; + else if(eta < 2.2) ea = 0.015; + else if(eta < 2.3) ea = 0.024; + else if(eta < 2.4) ea = 0.039; + else ea = 0.072; - float iso = gamma.neutralHadronIso; - iso = max(iso - rho*ea, (float)0.); + float iso = gamma.neutralHadronIso; + iso = max(iso - rho*ea, (float)0.); - return iso; + return iso; } float photonIso_corrected(susy::Photon gamma, float rho) { - float eta = fabs(gamma.caloPosition.Eta()); - float ea; + float eta = fabs(gamma.caloPosition.Eta()); + float ea; - if(eta < 1.0) ea = 0.148; - else if(eta < 1.479) ea = 0.130; - else if(eta < 2.0) ea = 0.112; - else if(eta < 2.2) ea = 0.216; - else if(eta < 2.3) ea = 0.262; - else if(eta < 2.4) ea = 0.260; - else ea = 0.266; + if(eta < 1.0) ea = 0.148; + else if(eta < 1.479) ea = 0.130; + else if(eta < 2.0) ea = 0.112; + else if(eta < 2.2) ea = 0.216; + else if(eta < 2.3) ea = 0.262; + else if(eta < 2.4) ea = 0.260; + else ea = 0.266; - float iso = gamma.photonIso; - iso = max(iso - rho*ea, (float)0.); + float iso = gamma.photonIso; + iso = max(iso - rho*ea, (float)0.); - return iso; + return iso; } +float d0correction( susy::Electron electron, susy::Event event ) { + TVector3 beamspot = event.vertices[0].position; + susy::Track track = event.tracks[electron.gsfTrackIndex]; + float d0 = track.d0() - beamspot.X()*sin(track.phi()) + beamspot.Y()*cos(track.phi()); + cout << "return " << d0 << endl; + return d0; +} - +float dZcorrection( susy::Electron electron, susy::Event event ) { + TVector3 beamspot = event.vertices[0].position; + susy::Track track = event.tracks[electron.gsfTrackIndex]; + + if(track.momentum.Pt() == 0.) return 1.e6; + float dz = (track.vertex.Z() - beamspot.Z()) - ((track.vertex.X() - beamspot.X())*track.momentum.Px() + (track.vertex.Y() - beamspot.Y())*track.momentum.Py()) / track.momentum.Pt() * (track.momentum.Pz() / track.momentum.Pt()); + return dz; +} float TreeWriter::getPtFromMatchedJet( susy::Photon myPhoton, susy::Event myEvent ) { /** @@ -132,7 +169,7 @@ float TreeWriter::getPtFromMatchedJet( s } else { susy::PFJetCollection& jetColl = pfJets_it->second; for(std::vector::iterator it = jetColl.begin(); - it != jetColl.end(); it++) { + it != jetColl.end(); ++it) { std::map::iterator s_it = it->jecScaleFactors.find("L2L3"); if (s_it == it->jecScaleFactors.end()) { std::cout << "JEC is not available for this jet!!!" << std::endl; @@ -149,14 +186,15 @@ float TreeWriter::getPtFromMatchedJet( s }// if, else if ( nearJets.size() == 0 ) { - //std::cout << "No jet with deltaR < .3 found, do not change photon_pt" << std::endl; + if( loggingVerbosity > 1 ) + std::cout << "No jet with deltaR < .3 found, do not change photon_pt" << std::endl; return myPhoton.momentum.Et(); } float pt = 0; float minPtDifferenz = 1E20; // should be very high for( std::vector::iterator it = nearJets.begin(), jetEnd = nearJets.end(); - it != jetEnd; it++ ) { + it != jetEnd; ++it ) { float ptDiff = fabs(myPhoton.momentum.Et() - it->momentum.Et()); if ( ptDiff < minPtDifferenz ) { minPtDifferenz = ptDiff; @@ -179,6 +217,8 @@ void TreeWriter::Loop() { // get number of events to be proceeded Long64_t nentries = inputTree->GetEntries(); + // store them in histo + eventNumbers->Fill( "Number of generated Events", nentries ); if(processNEvents <= 0 || processNEvents > nentries) processNEvents = nentries; if( loggingVerbosity > 0 ) @@ -190,45 +230,81 @@ void TreeWriter::Loop() { tree->Branch("photon", &photon); tree->Branch("jet", &jet); + tree->Branch("electron", &electron); + tree->Branch("muon", &muon); tree->Branch("met", &met, "met/F"); + tree->Branch("metPhi", &met_phi, "metPhi/F"); + tree->Branch("type1met", &type1met, "type1met/F"); + tree->Branch("type1metPhi", &type1met_phi, "type1metPhi/F"); tree->Branch("ht", &ht, "ht/F"); tree->Branch("nVertex", &nVertex, "nVertex/I"); - tree->Branch("nElectron", &nElectron, "nElectron/I"); + tree->Branch("pu_weight", &pu_weight, "pu_weight/F"); - for (long jentry=0; jentry < processNEvents; jentry++) { + for (long jentry=0; jentry < processNEvents; ++jentry) { if ( loggingVerbosity>0 && jentry%reportEvery==0 ) std::cout << jentry << " / " << processNEvents << std :: endl; inputTree->GetEntry(jentry); photon.clear(); jet.clear(); + electron.clear(); + muon.clear(); ht = 0; + // weights + if (pileupHisto == 0) { + pu_weight = 1.; + } else { + float trueNumInteractions = -1; + for( susy::PUSummaryInfoCollection::const_iterator iBX = event->pu.begin(); + iBX != event->pu.end() && trueNumInteractions < 0; ++iBX) { + if (iBX->BX == 0) + trueNumInteractions = iBX->trueNumInteractions; + } + pu_weight = pileupHisto->GetBinContent( pileupHisto->FindBin( trueNumInteractions ) ); + } + + // photons if( loggingVerbosity > 1 ) std::cout << "Process photons" << std::endl; std::map >::iterator phoMap = event->photons.find("photons"); for(std::vector::iterator it = phoMap->second.begin(); - it != phoMap->second.end() && phoMap != event->photons.end(); it++ ) { - if( ! it->isEB() && skim ) + it != phoMap->second.end() && phoMap != event->photons.end(); ++it ) { + if( !(it->isEB() || it->isEE()) && skim ) continue; thisphoton->pt = getPtFromMatchedJet( *it, *event ); - if( thisphoton->pt < 80 && skim ) - continue; - thisphoton->eta = it->momentum.Eta(); - thisphoton->phi = it->momentum.Phi(); + thisphoton->chargedIso = chargedHadronIso_corrected(*it, event->rho25); thisphoton->neutralIso = neutralHadronIso_corrected(*it, event->rho25); thisphoton->photonIso = photonIso_corrected(*it, event->rho25); - if ( it->r9 >= 1 && skim ) + + bool loose_photon_barrel = thisphoton->pt>20 + && it->isEB() + && it->passelectronveto + && it->hadTowOverEm<0.05 + && it->sigmaIetaIeta<0.012 + && thisphoton->chargedIso<2.6 + && thisphoton->neutralIso<3.5+0.04*thisphoton->pt + && thisphoton->photonIso<1.3+0.005*thisphoton->pt; + bool loose_photon_endcap = thisphoton->pt > 20 + && it->isEE() + && it->passelectronveto + && it->hadTowOverEm<0.05 + && it->sigmaIetaIeta<0.034 + && thisphoton->chargedIso<2.3 + && thisphoton->neutralIso<2.9+0.04*thisphoton->pt; + if(!(loose_photon_endcap || loose_photon_barrel || thisphoton->pt > 75 ) && skim ) continue; + thisphoton->eta = it->momentum.Eta(); + thisphoton->phi = it->momentum.Phi(); thisphoton->r9 = it->r9; thisphoton->sigmaIetaIeta = it->sigmaIetaIeta; thisphoton->hadTowOverEm = it->hadTowOverEm; thisphoton->pixelseed = it->nPixelSeeds; + thisphoton->conversionSafeVeto = it->passelectronveto; photon.push_back( *thisphoton ); - ht += thisphoton->pt; if( loggingVerbosity > 2 ) std::cout << " p_T, gamma = " << thisphoton->pt << std::endl; } @@ -248,7 +324,7 @@ void TreeWriter::Loop() { susy::PFJetCollection& jetColl = pfJets_it->second; for(std::vector::iterator it = jetColl.begin(); - it != jetColl.end(); it++) { + it != jetColl.end(); ++it) { std::map::iterator s_it = it->jecScaleFactors.find("L2L3"); if (s_it == it->jecScaleFactors.end()) { std::cout << "JEC is not available for this jet!!!" << std::endl; @@ -257,11 +333,27 @@ void TreeWriter::Loop() { float scale = s_it->second; TLorentzVector corrP4 = scale * it->momentum; - if(std::abs(corrP4.Eta()) > 3.0 && skim ) continue; + if(std::fabs(corrP4.Eta()) > 3.0 && skim ) continue; if(corrP4.Et() < 30 && skim ) continue; thisjet->pt = corrP4.Et(); thisjet->eta = corrP4.Eta(); thisjet->phi = corrP4.Phi(); + thisjet->bCSV = it->bTagDiscriminators[susy::kCSV]; + // jet composition + thisjet->chargedHadronEnergy = it->chargedHadronEnergy; + thisjet->neutralHadronEnergy = it->neutralHadronEnergy; + thisjet->photonEnergy = it->photonEnergy; + thisjet->electronEnergy = it->electronEnergy; + thisjet->muonEnergy = it->muonEnergy; + thisjet->HFHadronEnergy = it->HFHadronEnergy; + thisjet->HFEMEnergy = it->HFEMEnergy; + thisjet->chargedEmEnergy = it->chargedEmEnergy; + thisjet->chargedMuEnergy = it->chargedMuEnergy; + thisjet->neutralEmEnergy = it->neutralEmEnergy; + + if( loggingVerbosity > 2 ) + std::cout << " p_T, jet = " << thisjet->pt << std::endl; + jet.push_back( *thisjet ); ht += thisjet->pt; }// for jet @@ -269,18 +361,100 @@ void TreeWriter::Loop() { if( jet.size() < 2 && skim ) continue; std::sort( jet.begin(), jet.end(), tree::EtGreater); + if( loggingVerbosity > 1 ) + std::cout << "Found " << jet.size() << " jets" << std::endl; + // met std::map::iterator met_it = event->metMap.find("pfMet"); susy::MET* metobj = &(met_it->second); met = metobj->met(); + met_phi = metobj->mEt.Phi(); + if( loggingVerbosity > 2 ) + std::cout << " met = " << met << std::endl; + + std::map::iterator type1met_it = event->metMap.find("pfType1CorrectedMet"); + susy::MET* type1metobj = &(type1met_it->second); + type1met = type1metobj->met(); + type1met_phi = type1metobj->mEt.Phi(); + if( loggingVerbosity > 2 ) + std::cout << " type1met = " << type1met << std::endl; // electrons - //std::vector eVector = event->electrons["gsfElectrons"]; - //nElectron = eVector.size(); + tree::Particle* thiselectron = new tree::Particle(); + map >::iterator eleMap = event->electrons.find("gsfElectrons"); + if(eleMap == event->electrons.end() && loggingVerbosity > 0) { + cout << "gsfElectrons not found!" << endl; + } else { + cout << "now begin electrons " << endl; + for(vector::iterator it = eleMap->second.begin(); it < eleMap->second.end(); ++it) { + // for cuts see https://twiki.cern.ch/twiki/bin/viewauth/CMS/EgammaCutBasedIdentification for veto electrons + if( it->momentum.Pt() < 1 || it->momentum.Pt() > 1e6 ) + continue; // spike rejection + float iso = ( it->chargedHadronIso + max(it->neutralHadronIso+it->photonIso + - effectiveAreaElectron(it->momentum.Eta())*event->rho25, (Float_t)0. ) + ) / it->momentum.Pt(); + cout << iso << endl; + cout << d0correction( *it, *event ) << endl; + float test = d0correction( *it, *event ); + cout << test << endl; + float dZ = std::fabs( dZcorrection( *it, *event ) ); + cout <<" find e" << endl; + float d0 = 0.; + if ( it->isEB() ){ + if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.007 + || fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.8 + || it->sigmaIetaIeta > 0.01 + || it->hcalOverEcalBc > 0.15 + || d0 > 0.04 + || dZ > 0.2 + || iso > 0.15 ) + continue; + } + else if( it->isEE() ) { + if ( fabs(it->deltaEtaSuperClusterTrackAtVtx) > 0.01 + || fabs(it->deltaPhiSuperClusterTrackAtVtx) > 0.7 + || it->sigmaIetaIeta > 0.03 + || d0 > 0.04 + || dZ > 0.2 + || iso > 0.15 ) + continue; + } + else // not in barrel nor in endcap + continue; + // TODO: conversion rejection information not implemented yet, see twiki for more details + + thiselectron->pt = it->momentum.Pt(); + if( thiselectron->pt < 20 ) + continue; + if( loggingVerbosity > 2 ) + std::cout << " p_T, electron = " << it->momentum.Et() << std::endl; + thiselectron->eta = it->momentum.Eta(); + thiselectron->phi = it->momentum.Phi(); + electron.push_back( *thiselectron ); + } + } + if( loggingVerbosity > 1 ) + std::cout << "Found " << electron.size() << " electrons" << std::endl; + + // muons + std::vector mVector = event->muons["muons"]; + tree::Particle* thismuon = new tree::Particle(); + for( std::vector::iterator it = mVector.begin(); it != mVector.end(); ++it) { + if( !( it->isPFMuon() && ( it->isGlobalMuon() || it->isTrackerMuon() ) ) ) + continue; // see https://twiki.cern.ch/twiki/bin/view/CMSPublic/SWGuideMuonId#Loose_Muon + thismuon->pt = it->momentum.Et(); + if( thismuon->pt < 20 ) + continue; + thismuon->eta = it->momentum.Eta(); + thismuon->phi = it->momentum.Phi(); + muon.push_back( *thismuon ); + } + if( loggingVerbosity > 1 ) + std::cout << "Found " << muon.size() << " muons" << std::endl; - // vertices + // vertices nVertex = event->vertices.size(); if( ht < 450 && skim) @@ -290,9 +464,9 @@ void TreeWriter::Loop() { } // for jentry - - tree->Write(); outFile->cd(); + eventNumbers->Write(); + tree->Write(); outFile->Write(); outFile->Close(); }