From 7df26ad47487f707943d2524c67e07281be555e2 Mon Sep 17 00:00:00 2001 From: Federico Preiato Date: Tue, 6 Jun 2017 19:18:24 +0200 Subject: [PATCH] Analyzer with run and eta dependent HLTReco --- ...ijetCaloScoutingSelectionRunEtaDependent.C | 1043 +++++++++++++++++ 1 file changed, 1043 insertions(+) create mode 100644 src/analysisClass_mainDijetCaloScoutingSelectionRunEtaDependent.C diff --git a/src/analysisClass_mainDijetCaloScoutingSelectionRunEtaDependent.C b/src/analysisClass_mainDijetCaloScoutingSelectionRunEtaDependent.C new file mode 100644 index 000000000..53868c9e0 --- /dev/null +++ b/src/analysisClass_mainDijetCaloScoutingSelectionRunEtaDependent.C @@ -0,0 +1,1043 @@ +#define analysisClass_cxx +#include "analysisClass.h" +#include +#include +#include +#include +#include +#include +#include + +analysisClass::analysisClass(string * inputList, string * cutFile, string * treeName, string * outputFileName, string * cutEfficFile) + :baseClass(inputList, cutFile, treeName, outputFileName, cutEfficFile) +{ + std::cout << "analysisClass::analysisClass(): begins " << std::endl; + + std::string jetAlgo = getPreCutString1("jetAlgo"); + double rParam = getPreCutValue1("DeltaR"); + + if( jetAlgo == "AntiKt" ) + fjJetDefinition = JetDefPtr( new fastjet::JetDefinition(fastjet::antikt_algorithm, rParam) ); + else if( jetAlgo == "Kt" ) + fjJetDefinition = JetDefPtr( new fastjet::JetDefinition(fastjet::kt_algorithm, rParam) ); + else + fjJetDefinition = JetDefPtr( new fastjet::JetDefinition(fastjet::cambridge_algorithm, rParam) ); + + // For JECs + if( int(getPreCutValue1("useJECs"))==1 ) + { + std::cout << "Reapplying JECs on the fly" << std::endl; + std::cout << "Using IOV implementation for periodic JEC..." << std::endl; + + std::string L1Path = "data/Summer15_25nsV7_MC/Summer15_25nsV7_MC_L1FastJet_AK4PFchs.txt"; + std::string L2Path = "data/Summer15_25nsV7_MC/Summer15_25nsV7_MC_L2Relative_AK4PFchs.txt"; + std::string L3Path = "data/Summer15_25nsV7_MC/Summer15_25nsV7_MC_L3Absolute_AK4PFchs.txt"; + + L1Par = new JetCorrectorParameters(L1Path); + L2Par = new JetCorrectorParameters(L2Path); + L3Par = new JetCorrectorParameters(L3Path); + + std::vector vPar; + vPar.push_back(*L1Par); + vPar.push_back(*L2Par); + vPar.push_back(*L3Par); + + JetCorrector = new FactorizedJetCorrector(vPar); + + // procedure for 2015 CaloScouting data: + //std::string L1DATAPath = "data/74X_HLT_mcRun2_asymptotic_fromSpring15DR_v0_MC/74X_HLT_mcRun2_asymptotic_fromSpring15DR_v0_L1FastJet_AK4CaloHLT.txt"; + //std::string L2DATAPath = "data/74X_HLT_mcRun2_asymptotic_fromSpring15DR_v0_MC/74X_HLT_mcRun2_asymptotic_fromSpring15DR_v0_L2Relative_AK4CaloHLT.txt"; + //std::string L3DATAPath = "data/74X_HLT_mcRun2_asymptotic_fromSpring15DR_v0_MC/74X_HLT_mcRun2_asymptotic_fromSpring15DR_v0_L3Absolute_AK4CaloHLT.txt"; + //std::string L2L3ResidualPath = "data/Summer15_25nsV7_DATA/Summer15_25nsV7_DATA_L2L3Residual_AK4PF.txt" ; + + // procedure for 2016 CaloScouting data: + // 76X 2015 data + //std::string L2L3ResidualPath = "data/Fall15_25nsV2_DATA/Fall15_25nsV2_DATA_L2L3Residual_AK4PF.txt" ; + // 2016 data for ICHEP + //std::string L2L3ResidualPath = "data/Spring16_25nsV6_DATA/Spring16_25nsV6_DATA_L2L3Residual_AK4PF.txt"; + // 2016 data for Moriond + // std::string L2L3ResidualPath = "data/Spring16_V8_DATA/Spring16_25nsV8p2_DATA_L2L3Residual_AK4PFchs.txt"; + + // Latest 2016 JEC for PromptReco + /* + std::string L1DATAPath = "data/80X_dataRun2_HLT_frozen_v12/80X_dataRun2_HLT_frozen_v12_L1FastJet_AK4CaloHLT.txt"; + std::string L2DATAPath = "data/80X_dataRun2_HLT_frozen_v12/80X_dataRun2_HLT_frozen_v12_L2Relative_AK4CaloHLT.txt"; + std::string L3DATAPath = "data/80X_dataRun2_HLT_frozen_v12/80X_dataRun2_HLT_frozen_v12_L3Absolute_AK4CaloHLT.txt"; + std::string L2L3ResidualPath = "data/Spring16_25nsV10p2_DATA/Spring16_25nsV10p2_DATA_L2L3Residual_AK4PFchs.txt"; + + L1DATAPar = new JetCorrectorParameters(L1DATAPath); + L2DATAPar = new JetCorrectorParameters(L2DATAPath); + L3DATAPar = new JetCorrectorParameters(L3DATAPath); + L2L3Residual = new JetCorrectorParameters(L2L3ResidualPath); + + std::vector vPar_data; + //residuals are applied only to data + vPar_data.push_back(*L1DATAPar); + vPar_data.push_back(*L2DATAPar); + vPar_data.push_back(*L3DATAPar); + vPar_data.push_back(*L2L3Residual); + + JetCorrector_data = new FactorizedJetCorrector(vPar_data); + */ + + // Using IOV JEC now. Now that's 2016. See include/IOV.h for changing versions. + // Note that the IOV implementation can be used with just one JEC if wanted. Juska. + + /* + Ranges for 2016 run periods from DAS as of 4 Nov 16. + Periods with no certified luminosity are omitted. + B-v2: 273150 - 275376 + C-v2: 276282 - 276279 + D-v2: 276315 - 276653 + E-v2: 276824 - 277420 + F-v1: 277816 - 278808 + G-v1: 278816 - 280385 + H-v2: 282807 - 283885 + List obtained with commands like this: + [juska@lxplus069 workdir]$ das_client.py --query='run dataset=/JetHT/Run2016H-PromptReco-v2/MINIAOD' --limit 0 + */ + + // data + iov = new jec::IOV("AK4CaloHLT"); + iov->add("BCD",1,276823,true); // Using start-1 of run E instead as upper + iov->add("E",276824,277815,true); // Use start-1 of run F instead as upper + iov->add("F",277816,278801,true); // Note the division before F dataset ends as instructed by JEC group + iov->add("p2",278802,999999,true); // same as V8p2 for 29Sept ReReco (2016B --> 2016G) and for the last period of PromptReco (2016H) + + JetCorrector_data = new FactorizedJetCorrector(); // Will be filled later + + //uncertainty + //unc = new JetCorrectionUncertainty("data/Summer15_50nsV5_DATA/Summer15_50nsV5_DATA_Uncertainty_AK4PFchs.txt"); + //unc = new JetCorrectionUncertainty("data/Summer15_25nsV5_DATA/Summer15_25nsV5_DATA_Uncertainty_AK4PFchs.txt"); + //unc = new JetCorrectionUncertainty("data/Summer15_25nsV6_DATA/Summer15_25nsV6_DATA_Uncertainty_AK4PFchs.txt"); + // for 2015 CaloScouting + //unc = new JetCorrectionUncertainty("data/Summer15_25nsV7_DATA/Summer15_25nsV7_DATA_Uncertainty_AK4PFchs.txt"); + // for ICHEP 2016 CaloScouting + //unc = new JetCorrectionUncertainty("data/Spring16_25nsV6_DATA/Spring16_25nsV6_DATA_Uncertainty_AK4PFchs.txt"); + // for 2016 CaloScouting + unc = new JetCorrectionUncertainty("data/Spring16_25nsV10_DATA/Spring16_25nsV10p2_DATA_Uncertainty_AK4PFchs.txt"); + + } + + std::cout << "analysisClass::analysisClass(): ends " << std::endl; +} + +analysisClass::~analysisClass() +{ + std::cout << "analysisClass::~analysisClass(): begins " << std::endl; + + std::cout << "analysisClass::~analysisClass(): ends " << std::endl; +} + +void analysisClass::Loop() +{ + std::cout << "analysisClass::Loop() begins" <Sumw2(); + // TH1F *h_nVtx = new TH1F ("h_nVtx","",30,0,30); + // h_nVtx->Sumw2(); + // TH1F *h_trueVtx = new TH1F ("h_trueVtx","",40,0,40); + // h_trueVtx->Sumw2(); + // TH1F *h_pT1stJet = new TH1F ("h_pT1stJet","",100,0,3000); + // h_pT1stJet->Sumw2(); + // TH1F *h_pT2ndJet = new TH1F ("h_pT2ndJet","",100,0,3000); + // h_pT2ndJet->Sumw2(); + // TH1F *h_eta1stJet = new TH1F ("h_eta1stJet","",5,-2.5,2.5); + // h_eta1stJet->Sumw2(); + // TH1F *h_eta2ndJet = new TH1F ("h_eta2ndJet","",5,-2.5,2.5); + // h_eta2ndJet->Sumw2(); + // TH1F *h_DijetMass = new TH1F ("h_DijetMass","",600,0,6000); + // h_DijetMass->Sumw2(); + // TH1F *h_DeltaETAjj = new TH1F ("h_DeltaETAjj","",120,0,3.); + // h_DeltaETAjj->Sumw2(); + + // variable binning for mjj trigger efficiency plots + const int nMassBins = 103; + + double massBoundaries[nMassBins+1] = {1, 3, 6, 10, 16, 23, 31, 40, 50, 61, 74, 88, 103, 119, 137, 156, 176, 197, 220, 244, 270, 296, 325, + 354, 386, 419, 453, 489, 526, 565, 606, 649, 693, 740, 788, 838, 890, 944, 1000, 1058, 1118, 1181, 1246, 1313, 1383, 1455, 1530, 1607, 1687, + 1770, 1856, 1945, 2037, 2132, 2231, 2332, 2438, 2546, 2659, 2775, 2895, 3019, 3147, 3279, 3416, 3558, 3704, 3854, 4010, 4171, 4337, 4509, + 4686, 4869, 5058, 5253, 5455, 5663, 5877, 6099, 6328, 6564, 6808, 7060, 7320, 7589, 7866, 8152, 8447, 8752, 9067, 9391, 9726, 10072, 10430, + 10798, 11179, 11571, 11977, 12395, 12827, 13272, 13732, 14000}; + + // char* HLTname[50] = {"noTrig","PFHT475","PFHT800","PFHT650MJJ900","PFHT800_OR_PFHT650MJJ900","PFHT800_noPFHT475", + // "Mu45Eta2p1", "PFHT800AndMu45Eta2p1"}; + // TH1F* h_mjj_HLTpass[8]; + // char name_histoHLT[50]; + // for (int i=0; i<8; i++){ + // sprintf(name_histoHLT,"h_mjj_HLTpass_%s",HLTname[i]); + // h_mjj_HLTpass[i]= new TH1F(name_histoHLT,"",103,massBoundaries); + // } + + //For trigger efficiency measurements + //No trigger selection applied (full offline selection applied) + TH1F* h_mjj_NoTrigger_1GeVbin = new TH1F("h_mjj_NoTrigger_1GeVbin","",14000,0,14000); + TH1F* h_mjj_NoTrigger = new TH1F("h_mjj_NoTrigger","",103,massBoundaries); + //HLT + TH1F* h_mjj_HLTpass_CaloJet40_CaloScouting_PFScouting = new TH1F("h_mjj_HLTpass_CaloJet40_CaloScouting_PFScouting","",103,massBoundaries); + TH1F* h_mjj_HLTpass_L1HTT_CaloScouting_PFScouting = new TH1F("h_mjj_HLTpass_L1HTT_CaloScouting_PFScouting","",103,massBoundaries); + TH1F* h_mjj_HLTpass_CaloScoutingHT250 = new TH1F("h_mjj_HLTpass_CaloScoutingHT250","",103,massBoundaries); + + /////////initialize variables + + Long64_t nentries = fChain->GetEntriesFast(); + std::cout << "analysisClass::Loop(): nentries = " << nentries << std::endl; + + ////// The following ~7 lines have been taken from rootNtupleClass->Loop() ///// + ////// If the root version is updated and rootNtupleClass regenerated, ///// + ////// these lines may need to be updated. ///// + Long64_t nbytes = 0, nb = 0; + for (Long64_t jentry=0; jentryGetEntry(jentry); nbytes += nb; + if(jentry < 10 || jentry%1000 == 0) std::cout << "analysisClass::Loop(): jentry = " << jentry << std::endl; + // if (Cut(ientry) < 0) continue; + + ////////////////////// User's code starts here /////////////////////// + + ///Stuff to be done for every event + + size_t no_jets_ak4=jetPtAK4->size(); + + vector widejets; + vector widejets_noCorr; + TLorentzVector wj1, wj2, wdijet; + TLorentzVector wj1_noCorr, wj2_noCorr, wdijet_noCorr; + TLorentzVector wj1_shift, wj2_shift, wdijet_shift; + + vector AK4jets; + TLorentzVector ak4j1, ak4j2, ak4dijet; + + resetCuts(); + + // //find intime BX + // int idx_InTimeBX=-1; + // for(size_t j=0; jsize(); ++j) + // { + // //cout << PileupOriginBX->at(j) << endl; + // if(PileupOriginBX->at(j)==0) + // { + // idx_InTimeBX = j; + // //cout << "idx_InTimeBX: " << idx_InTimeBX << endl; + // } + // } + + std::vector jecFactors; + std::vector jecUncertainty; + std::vector idCaloJet; // CaloJet ID + // new JECs could change the jet pT ordering. the vector below + // holds sorted jet indices after the new JECs had been applied + std::vector sortedJetIdx; + + if( int(getPreCutValue1("useJECs"))==1 ) + { + // sort jets by increasing pT + std::multimap sortedJets; + for(size_t j=0; jsetJetEta(jetEtaAK4->at(j)); + JetCorrector->setJetPt(jetPtAK4->at(j)/jetJecAK4->at(j)); //pTraw + JetCorrector->setJetA(jetAreaAK4->at(j)); + JetCorrector->setRho(rho); + + JetCorrector_data = iov->get(runNo); // Get IOV dependent JEC + + JetCorrector_data->setJetEta(jetEtaAK4->at(j)); + JetCorrector_data->setJetPt(jetPtAK4->at(j)/jetJecAK4->at(j)); //pTraw + JetCorrector_data->setJetA(jetAreaAK4->at(j)); + JetCorrector_data->setRho(rho); + + + //nominal value of JECs + double correction;//, old_correction, nominal_correction; + //if( int(getPreCutValue1("shiftJECs"))==0 ){ + if (isData == 1) correction = JetCorrector_data->getCorrection(); + else correction = JetCorrector->getCorrection(); + //nominal_correction=correction; + //old_correction = jetJecAK4->at(j); + //} + + // std::cout << "Correction applied = "<< correction << std::endl; + + //JEC uncertainties + unc->setJetEta(jetEtaAK4->at(j)); + unc->setJetPt(jetPtAK4->at(j)/jetJecAK4->at(j)*correction); + double uncertainty = unc->getUncertainty(true); + jecUncertainty.push_back(uncertainty); + + // std::cout << "run:" << runNo << " lumi:" << lumi << " event:" << evtNo << " jet pt:" << jetPtAK4->at(j)/jetJecAK4->at(j)*correction << " correction:" << correction << " uncertainty:" << uncertainty << " nominal correction:" << nominal_correction << " old correction: " << old_correction << std::endl; + //use "shifted" JECs for study of systematic uncertainties + if( int(getPreCutValue1("shiftJECs"))==1 ){ + //flat shift + //if (isData == 1) correction = JetCorrector_data->getCorrection() * getPreCutValue2("shiftJECs"); + //else correction = JetCorrector->getCorrection() * getPreCutValue2("shiftJECs"); + //shift of the corresponding unc + correction = correction + getPreCutValue2("shiftJECs")*uncertainty*correction; + // std::cout << "run:" << runNo << " lumi:" << lumi << " event:" << evtNo << " jet pt:" << jetPtAK3->at(j)/jetJecAK4->at(j)*correction << " correction:" << correction << " uncertainty:" << uncertainty << std::endl << std::endl; + + } + + jecFactors.push_back(correction); + + bool idval = false; + if( jetHadfAK4->at(j)at(j)at(j)/jetJecAK4->at(j))*correction, j)); + + } + + // get jet indices in decreasing pT order + for(std::multimap::const_reverse_iterator it = sortedJets.rbegin(); it != sortedJets.rend(); ++it) + sortedJetIdx.push_back(it->second); + + } + else if( int(getPreCutValue1("noJECs"))==1 ) + { + // sort jets by increasing pT + std::multimap sortedJets; + for(size_t j=0; jat(j)at(j)at(j)/jetJecAK4->at(j)), j)); //raw + } + // get jet indices in decreasing pT order + for(std::multimap::const_reverse_iterator it = sortedJets.rbegin(); it != sortedJets.rend(); ++it) + sortedJetIdx.push_back(it->second); + } + else + { + for(size_t j=0; jat(j)); + jecUncertainty.push_back(0.); + + bool idval = false; + if( jetHadfAK4->at(j)at(j)=2){ + // if( !( fabs(jetEtaAK4->at(0)) < getPreCutValue1("jetFidRegion") && idTAK4->at(0) == getPreCutValue1("tightJetID"))){ + // std::cout << " JET 0 FAIL " << jetEtaAK4->at(0) << " JET 0 ID " << idTAK4->at(0) << std::endl; + // } + // if( !( fabs(jetEtaAK4->at(1)) < getPreCutValue1("jetFidRegion") && idTAK4->at(1) == getPreCutValue1("tightJetID"))){ + // std::cout << " JET 1 FAIL " << jetEtaAK4->at(1) << " JET 1 ID " << idTAK4->at(1) << std::endl; + // } + // } + + //count ak4 jets passing pt threshold and id criteria + int Nak4 = 0; + double HTak4 = 0; + + for(size_t ijet=0; ijetat(sortedJetIdx[ijet]) == getPreCutValue1("tightJetID") // figure out ARTUR + && idCaloJet[sortedJetIdx[ijet]] == getPreCutValue1("tightJetID") + && (jecFactors[sortedJetIdx[ijet]]/jetJecAK4->at(sortedJetIdx[ijet]))*jetPtAK4->at(sortedJetIdx[ijet]) > getPreCutValue1("ptCut")) + { + Nak4 += 1; + HTak4 += (jecFactors[sortedJetIdx[ijet]]/jetJecAK4->at(sortedJetIdx[ijet]))*jetPtAK4->at(sortedJetIdx[ijet]); + } + } + + + if( int(getPreCutValue1("useFastJet"))==1 ) + { + // vector of ak4 jets used for wide jet clustering + std::vector fjInputs, fjInputs_shift; + + for(size_t j=0; jat(sortedJetIdx[j])) < getPreCutValue1("jetFidRegion") + // && idTAK4->at(sortedJetIdx[j]) == getPreCutValue1("tightJetID") + && idCaloJet[sortedJetIdx[j]] == getPreCutValue1("tightJetID") + ) ) continue; + + double rescale = (jecFactors[sortedJetIdx[j]]/jetJecAK4->at(sortedJetIdx[j])); + + if( j==0 && !( rescale*jetPtAK4->at(sortedJetIdx[j]) > getPreCutValue1("pt0Cut")) ) continue; + else if( j==1 && !( rescale*jetPtAK4->at(sortedJetIdx[j]) > getPreCutValue1("pt1Cut")) ) continue; + else if( !( rescale*jetPtAK4->at(sortedJetIdx[j]) > getPreCutValue1("ptCut")) ) continue; + + TLorentzVector tempJet, tempJet_shift; + + tempJet.SetPtEtaPhiM( rescale*jetPtAK4->at(sortedJetIdx[j]) , jetEtaAK4->at(sortedJetIdx[j]) , jetPhiAK4->at(sortedJetIdx[j]) , rescale*jetMassAK4->at(sortedJetIdx[j])); + tempJet_shift.SetPtEtaPhiM( (1+jecUncertainty[sortedJetIdx[j]])* rescale*jetPtAK4->at(sortedJetIdx[j]) , jetEtaAK4->at(sortedJetIdx[j]) , jetPhiAK4->at(sortedJetIdx[j]) , (1+jecUncertainty[sortedJetIdx[j]])* rescale*jetMassAK4->at(sortedJetIdx[j])); + + fjInputs.push_back(fastjet::PseudoJet(tempJet.Px(),tempJet.Py(),tempJet.Pz(),tempJet.E())); + fjInputs_shift.push_back(fastjet::PseudoJet(tempJet_shift.Px(),tempJet_shift.Py(),tempJet_shift.Pz(),tempJet_shift.E())); + } + + fjClusterSeq = ClusterSequencePtr( new fastjet::ClusterSequence( fjInputs, *fjJetDefinition ) ); + fjClusterSeq_shift = ClusterSequencePtr( new fastjet::ClusterSequence( fjInputs_shift, *fjJetDefinition ) ); + + std::vector inclusiveWideJets = fastjet::sorted_by_pt( fjClusterSeq->inclusive_jets(0.) ); + std::vector inclusiveWideJets_shift = fastjet::sorted_by_pt( fjClusterSeq_shift->inclusive_jets(0.) ); + + if( inclusiveWideJets.size()>1 ) + { + wj1.SetPxPyPzE(inclusiveWideJets.at(0).px(), inclusiveWideJets.at(0).py(), inclusiveWideJets.at(0).pz(), inclusiveWideJets.at(0).e()); + wj2.SetPxPyPzE(inclusiveWideJets.at(1).px(), inclusiveWideJets.at(1).py(), inclusiveWideJets.at(1).pz(), inclusiveWideJets.at(1).e()); + wj1_shift.SetPxPyPzE(inclusiveWideJets_shift.at(0).px(), inclusiveWideJets_shift.at(0).py(), inclusiveWideJets_shift.at(0).pz(), inclusiveWideJets_shift.at(0).e()); + wj2_shift.SetPxPyPzE(inclusiveWideJets_shift.at(1).px(), inclusiveWideJets_shift.at(1).py(), inclusiveWideJets_shift.at(1).pz(), inclusiveWideJets_shift.at(1).e()); + } + } + else + { + TLorentzVector wj1_tmp, wj2_tmp; + TLorentzVector wj1_shift_tmp, wj2_shift_tmp; + double wideJetDeltaR_ = getPreCutValue1("DeltaR"); + + if(no_jets_ak4>=2) + { + if( fabs(jetEtaAK4->at(sortedJetIdx[0])) < getPreCutValue1("jetFidRegion") + && (jecFactors[sortedJetIdx[0]]/jetJecAK4->at(sortedJetIdx[0]))*jetPtAK4->at(sortedJetIdx[0]) > getPreCutValue1("pt0Cut")) + { + if( fabs(jetEtaAK4->at(sortedJetIdx[1])) < getPreCutValue1("jetFidRegion") + && (jecFactors[sortedJetIdx[1]]/jetJecAK4->at(sortedJetIdx[1]))*jetPtAK4->at(sortedJetIdx[1]) > getPreCutValue1("pt1Cut")) + { + TLorentzVector jet1, jet2, jet1_shift, jet2_shift; + jet1.SetPtEtaPhiM( (jecFactors[sortedJetIdx[0]]/jetJecAK4->at(sortedJetIdx[0])) *jetPtAK4->at(sortedJetIdx[0]) + ,jetEtaAK4->at(sortedJetIdx[0]),jetPhiAK4->at(sortedJetIdx[0]) + , (jecFactors[sortedJetIdx[0]]/jetJecAK4->at(sortedJetIdx[0])) * jetMassAK4->at(sortedJetIdx[0])); + jet2.SetPtEtaPhiM( (jecFactors[sortedJetIdx[1]]/jetJecAK4->at(sortedJetIdx[1])) *jetPtAK4->at(sortedJetIdx[1]) + ,jetEtaAK4->at(sortedJetIdx[1]),jetPhiAK4->at(sortedJetIdx[1]) + , (jecFactors[sortedJetIdx[1]]/jetJecAK4->at(sortedJetIdx[1])) * jetMassAK4->at(sortedJetIdx[1])); + jet1_shift.SetPtEtaPhiM( (1+jecUncertainty[sortedJetIdx[0]])*(jecFactors[sortedJetIdx[0]]/jetJecAK4->at(sortedJetIdx[0])) *jetPtAK4->at(sortedJetIdx[0]) + ,jetEtaAK4->at(sortedJetIdx[0]),jetPhiAK4->at(sortedJetIdx[0]) + , (1+jecUncertainty[sortedJetIdx[0]])*(jecFactors[sortedJetIdx[0]]/jetJecAK4->at(sortedJetIdx[0])) * jetMassAK4->at(sortedJetIdx[0])); + jet2_shift.SetPtEtaPhiM( (1+jecUncertainty[sortedJetIdx[1]])* (jecFactors[sortedJetIdx[1]]/jetJecAK4->at(sortedJetIdx[1])) *jetPtAK4->at(sortedJetIdx[1]) + ,jetEtaAK4->at(sortedJetIdx[1]),jetPhiAK4->at(sortedJetIdx[1]) + , (1+jecUncertainty[sortedJetIdx[1]])*(jecFactors[sortedJetIdx[1]]/jetJecAK4->at(sortedJetIdx[1])) * jetMassAK4->at(sortedJetIdx[1])); + + for(Long64_t ijet=0; ijetat(sortedJetIdx[ijet])) < getPreCutValue1("jetFidRegion") + && idCaloJet[sortedJetIdx[ijet]] == getPreCutValue1("tightJetID") + && (jecFactors[sortedJetIdx[ijet]]/jetJecAK4->at(sortedJetIdx[ijet]))*jetPtAK4->at(sortedJetIdx[ijet]) > getPreCutValue1("ptCut")) + { + TLorentzVector currentJet, currentJet_shift; + currentJet.SetPtEtaPhiM( (jecFactors[sortedJetIdx[ijet]]/jetJecAK4->at(sortedJetIdx[ijet])) *jetPtAK4->at(sortedJetIdx[ijet]) + ,jetEtaAK4->at(sortedJetIdx[ijet]),jetPhiAK4->at(sortedJetIdx[ijet]) + , (jecFactors[sortedJetIdx[ijet]]/jetJecAK4->at(sortedJetIdx[ijet])) *jetMassAK4->at(sortedJetIdx[ijet])); + currentJet_shift.SetPtEtaPhiM( (1+jecUncertainty[sortedJetIdx[ijet]])*(jecFactors[sortedJetIdx[ijet]]/jetJecAK4->at(sortedJetIdx[ijet])) *jetPtAK4->at(sortedJetIdx[ijet]) + ,jetEtaAK4->at(sortedJetIdx[ijet]),jetPhiAK4->at(sortedJetIdx[ijet]) + , (1+jecUncertainty[sortedJetIdx[ijet]])*(jecFactors[sortedJetIdx[ijet]]/jetJecAK4->at(sortedJetIdx[ijet])) *jetMassAK4->at(sortedJetIdx[ijet])); + + double DeltaR1 = currentJet.DeltaR(jet1); + double DeltaR2 = currentJet.DeltaR(jet2); + + if(DeltaR1 < DeltaR2 && DeltaR1 < wideJetDeltaR_) + { + wj1_tmp += currentJet; + wj1_shift_tmp += currentJet_shift; + } + else if(DeltaR2 < wideJetDeltaR_) + { + wj2_tmp += currentJet; + wj2_shift_tmp += currentJet_shift; + } + } // if AK4 jet passes fid and jetid. + } //end of ak4 jet loop + + // if(wj1_tmp.Pt()==0 && wj2_tmp.Pt() ==0) + // std::cout << " wj1_tmp.Pt() IN " < wj2_tmp.Pt()) + { + wj1 = wj1_tmp; + wj2 = wj2_tmp; + wj1_shift = wj1_shift_tmp; + wj2_shift = wj2_shift_tmp; + } + else + { + wj1 = wj2_tmp; + wj2 = wj1_tmp; + wj1_shift = wj2_shift_tmp; + wj2_shift = wj1_shift_tmp; + } + } + + + double MJJWide = 0; + double MJJWideNoCorr = 0; + double DeltaEtaJJWide = 0; + double DeltaPhiJJWide = 0; + double MJJWide_shift = 0; + float corr1 = 1.; + float corr2 = 1.; + if( wj1.Pt()>0 && wj2.Pt()>0 ) + { + /* + // 2015 bias correction from Mikko/Federico + float p0 = - 0.63; + float p1 = 0.636; + float p2 = - 1.0; + + float f1 = p0 + p1 * pow( 0.01 * wj1.Pt() , p2); + float f2 = p0 + p1 * pow( 0.01 * wj2.Pt() , p2); + */ + + /* + // old 2016 bias correction from Mikko/Federico + // (page 4 https://indico.cern.ch/event/546408/contributions/2217944/attachments/1298388/1936977/Giugno-24-2016_-_CaloScouting.pdf) + float p0 = 0.419; + float p1 = -5; + float p2 = -0.188; + float p3 = -0.5; + float p4 = 1.15; + float p5 = 399; + float p6 = 252; + + float f1 = p0 + p1 * pow( wj1.Pt() , p2) + p3/wj1.Pt() + p4 * exp( -0.5 * ((wj1.Pt() - p5)/p6) * ((wj1.Pt() - p5)/p6) ); + float f2 = p0 + p1 * pow( wj2.Pt() , p2) + p3/wj2.Pt() + p4 * exp( -0.5 * ((wj2.Pt() - p5)/p6) * ((wj2.Pt() - p5)/p6) ); + */ + + /* + // new 2016 bias correction from Federico + // (page 10 https://www.dropbox.com/s/7sporqeim01675d/Luglio_20_2016_CaloScouting.pdf?dl=1) + // flattened above 993.264 (point of zero slope) + float p0 = -31.7198; + float p1 = 8.58611; + float p2 = -0.622092; + + float f1 = 0; + float f2 = 0; + if (wj1.Pt() >= 993.264) + f1 = p0 + p1 * log( 993.264 ) + p2 * log( 993.264 ) * log( 993.264 ) ; + else + f1 = p0 + p1 * log( wj1.Pt() ) + p2 * log( wj1.Pt() ) * log( wj1.Pt() ) ; + + if (wj2.Pt() >= 993.264) + f2 = p0 + p1 * log( 993.264 ) + p2 * log( 993.264 ) * log( 993.264 ) ; + else + f2 = p0 + p1 * log( wj2.Pt() ) + p2 * log( wj2.Pt() ) * log( wj2.Pt() ) ; + */ + + // new 2016 bias correction -- full 2016 dataset + // TF1 *f3 = new TF1("f3","[0]*pow(x,[1]) + [2]*pow(max(0.,log(x/700.)),1.2)",180,2500) + + /* + // HLT-Reco correction for Moriond17 + + float p0 = -1.55133e+02; + float p1 = -7.09058e-01; + float p2 = -3.11990e-01; + + float f1 = 0; + float f2 = 0; + + f1 = p0 * pow(wj1.Pt(),p1) + p2 * pow(max(0.,log( wj1.Pt()/700. )),1.2) ; + f2 = p0 * pow(wj2.Pt(),p1) + p2 * pow(max(0.,log( wj2.Pt()/700. )),1.2) ; + */ + + /* + // HLT-Reco correction after Moriond17 - IOV system also for CaloJets + float p0 = -26.0058; + float p1 = 7.32690; + float p2 = -0.552001; + + float f1 = 0; + float f2 = 0; + + f1 = p0 + p1 * log( wj1.Pt() ) + p2 * log( wj1.Pt() ) * log( wj1.Pt() ) ; + f2 = p0 + p1 * log( wj2.Pt() ) + p2 * log( wj2.Pt() ) * log( wj2.Pt() ) ; + */ + + /* + // new 2016 bias correction with IOV system for Res and HLT-Reco correction + float p0; + float p1; + float p2; + + // "BCD": 1,276823 + // "E": 276824,277815 + // "F": 277816,278801 + // "p2": 278802,999999 + + if( runNo >= 1 && runNo <= 276823 ){ // Run2016BCD + // std::cout<< "Run2016BCD" << std::endl; + p0 = -3.03601e+01; + p1 = 8.54708e+00; + p2 = -6.41985e-01; + } + else if( runNo >= 276824 && runNo <= 277815 ){ // Run2016E + // std::cout<< "Run2016E" << std::endl; + p0 = -3.14290e+01; + p1 = 9.58863e+00; + p2 = -7.47892e-01; + } + else if( runNo >= 277816 && runNo <= 278801 ){ // Run2016F + // std::cout<< "Run2016F" << std::endl; + p0 = -3.21730e+01; + p1 = 1.01856e+01; + p2 = -8.26144e-01; + } + else if( runNo >= 278802 && runNo <= 999999 ){ // Run2016G + // std::cout<< "Run2016G" << std::endl; + p0 = -2.82339e+01; + p1 = 7.95804e+00; + p2 = -5.99060e-01; + }else{ + std::cout<< "runNo OUT OF RANGE" << std::endl; + } + + float f1 = 0; + float f2 = 0; + + f1 = p0 + p1 * log( wj1.Pt() ) + p2 * log( wj1.Pt() ) * log( wj1.Pt() ) ; + f2 = p0 + p1 * log( wj2.Pt() ) + p2 * log( wj2.Pt() ) * log( wj2.Pt() ) ; + */ + + // new 2016 bias correction with IOV system for Res and HLT-Reco correction + eta-dependent + + // DEBUGGING + // std::cout<< "BEFORE CORRECTION"<< std::endl; + // std::cout<< "wj1"<< std::endl; + // std::cout<< "Pt: "<< wj1.Pt()<< " Eta: "< f_HLTReco; + // loop on widejets + for(int ijet=0; ijet<2; ijet++){ + + TLorentzVector wj_i; + if( ijet == 0 ) wj_i = wj1; + if( ijet == 1 ) wj_i = wj2; + + // std::cout<< "wj_i"<< std::endl; + // std::cout<< "Pt: "<< wj_i.Pt()<< " Eta: "<= 1 && runNo <= 276823 ){ // Run2016BCD + // std::cout<< "Run2016BCD" << std::endl; + if( fabs(wj_i.Eta())<=1.3 ){ + // std::cout<< "|eta| < 1.3" << std::endl; + p0 = -2.97322e+01; + p1 = 8.56657e+00; + p2 = -6.56580e-01; + }else if( fabs(wj_i.Eta())>1.3 && fabs(wj_i.Eta()) <=2.5 ){ + // std::cout<< "1.3 < |eta| < 2.5" << std::endl; + p0 = -3.52634e+01; + p1 = 1.00689e+01; + p2 = -7.82283e-01; + } + } + else if( runNo >= 276824 && runNo <= 277815 ){ // Run2016E + // std::cout<< "Run2016E" << std::endl; + if( fabs(wj_i.Eta()) <= 1.3 ){ + // std::cout<< "|eta| < 1.3" << std::endl; + p0 = -3.00851e+01; + p1 = 9.50250e+00; + p2 = -7.62598e-01; + }else if( fabs(wj_i.Eta()) > 1.3 && fabs(wj_i.Eta()) <= 2.5 ){ + // std::cout<< "1.3 < |eta| < 2.5" << std::endl; + p0 = -3.91216e+01; + p1 = 1.18063e+01; + p2 = -9.33997e-01; + } + } + else if( runNo >= 277816 && runNo <= 278801 ){ // Run2016F + // std::cout<< "Run2016F" << std::endl; + if( fabs(wj_i.Eta()) <= 1.3 ){ + // std::cout<< "|eta| < 1.3" << std::endl; + p0 = -2.67199e+01; + p1 = 8.66466e+00; + p2 = -7.19131e-01; + }else if( fabs(wj_i.Eta()) > 1.3 && fabs(wj_i.Eta()) <= 2.5 ){ + // std::cout<< "1.3 < |eta| < 2.5" << std::endl; + p0 = -4.83545e+01; + p1 = 1.51752e+01; + p2 = -1.22605e+00; + } + } + else if( runNo >= 278802 && runNo <= 999999 ){ // Run2016G + // std::cout<< "Run2016G" << std::endl; + if( fabs(wj_i.Eta()) <= 1.3 ){ + // std::cout<< "|eta| < 1.3" << std::endl; + p0 = -2.73704e+01; + p1 = 7.81992e+00; + p2 = -5.96448e-01; + }else if( fabs(wj_i.Eta()) > 1.3 && fabs(wj_i.Eta()) <= 2.5 ){ + // std::cout<< "1.3 < |eta| < 2.5" << std::endl; + p0 = -3.20948e+01; + p1 = 9.11822e+00 ; + p2 = -6.98530e-01; + } + }else{ //just a check + std::cout<< "runNo OUT OF RANGE" << std::endl; + } + + float f_i = 0; + + f_i = p0 + p1 * log( wj_i.Pt() ) + p2 * log( wj_i.Pt() ) * log( wj_i.Pt() ) ; + + f_HLTReco.push_back( f_i ); + } + + // std::cout<< "f1: " << f_HLTReco[0] << std::endl; + // std::cout<< "f2: " << f_HLTReco[1] << std::endl; + + + ///////////////////////////////////////// + + // corr1 = 1. / (1. + 0.01*f1); + // corr2 = 1. / (1. + 0.01*f2); + corr1 = 1. / (1. + 0.01*f_HLTReco[0]); + corr2 = 1. / (1. + 0.01*f_HLTReco[1]); + + // Get MJJ before corrections + wj1_noCorr = TLorentzVector(wj1); + wj2_noCorr = TLorentzVector(wj2); + wdijet_noCorr = wj1_noCorr + wj2_noCorr; + MJJWideNoCorr = wdijet_noCorr.M(); + + // Put widejets in the container + widejets_noCorr.push_back( wj1_noCorr ); + widejets_noCorr.push_back( wj2_noCorr ); + + // Apply corrections + wj1 = wj1*corr1; + wj2 = wj2*corr2; + + // std::cout<< "AFTER CORRECTION"<< std::endl; + // std::cout<< "wj1"<< std::endl; + // std::cout<< "Pt: "<< wj1.Pt()<< " Eta: "<=2) + //cout << "eta j1 " << jetEtaAK4->at(sortedJetIdx[0]) << endl; + //cout << "pt j1 " << (jecFactors[sortedJetIdx[0]]/jetJecAK4->at(sortedJetIdx[0])) *jetPtAK4->at(sortedJetIdx[0]) << endl; + { + if( fabs(jetEtaAK4->at(sortedJetIdx[0])) < getPreCutValue1("jetFidRegion") + && (jecFactors[sortedJetIdx[0]]/jetJecAK4->at(sortedJetIdx[0]))*jetPtAK4->at(sortedJetIdx[0]) > getPreCutValue1("pt0Cut")) + { + if( fabs(jetEtaAK4->at(sortedJetIdx[1])) < getPreCutValue1("jetFidRegion") + && (jecFactors[sortedJetIdx[1]]/jetJecAK4->at(sortedJetIdx[1]))*jetPtAK4->at(sortedJetIdx[1]) > getPreCutValue1("pt1Cut")) + { + //cout << "filling ak4j1 and ak4j2" << endl; + //cout << "pt ak4 j1 = " << (jecFactors[sortedJetIdx[0]]/jetJecAK4->at(sortedJetIdx[0])) *jetPtAK4->at(sortedJetIdx[0]) << endl; + ak4j1.SetPtEtaPhiM( (jecFactors[sortedJetIdx[0]]/jetJecAK4->at(sortedJetIdx[0])) *jetPtAK4->at(sortedJetIdx[0]) + ,jetEtaAK4->at(sortedJetIdx[0]) + ,jetPhiAK4->at(sortedJetIdx[0]) + , (jecFactors[sortedJetIdx[0]]/jetJecAK4->at(sortedJetIdx[0])) *jetMassAK4->at(sortedJetIdx[0])); + ak4j2.SetPtEtaPhiM( (jecFactors[sortedJetIdx[1]]/jetJecAK4->at(sortedJetIdx[1])) *jetPtAK4->at(sortedJetIdx[1]) + ,jetEtaAK4->at(sortedJetIdx[1]) + ,jetPhiAK4->at(sortedJetIdx[1]) + , (jecFactors[sortedJetIdx[1]]/jetJecAK4->at(sortedJetIdx[1])) *jetMassAK4->at(sortedJetIdx[1])); + } + } + } + + double MJJAK4 = 0; + double DeltaEtaJJAK4 = 0; + double DeltaPhiJJAK4 = 0; + + //std::cout << "ak4j1.Pt()=" << ak4j1.Pt() << " ak4j2.Pt()=" << ak4j2.Pt() << std::endl; + if( ak4j1.Pt()>0 && ak4j2.Pt()>0 ) + { + // Create dijet system + ak4dijet = ak4j1 + ak4j2; + MJJAK4 = ak4dijet.M(); + DeltaEtaJJAK4 = fabs(ak4j1.Eta()-ak4j2.Eta()); + DeltaPhiJJAK4 = fabs(ak4j1.DeltaPhi(ak4j2)); + + // Put widejets in the container + AK4jets.push_back( ak4j1 ); + AK4jets.push_back( ak4j2 ); + } + + //== Fill Variables == + fillVariableWithValue("isData",isData); + fillVariableWithValue("run",runNo); + fillVariableWithValue("event",evtNo); + fillVariableWithValue("lumi",lumi); + fillVariableWithValue("nVtx",nvtx); + fillVariableWithValue("nJet",widejets.size()); + fillVariableWithValue("Nak4",Nak4); + fillVariableWithValue ( "PassJSON", passJSON (runNo, lumi, isData)); + + //directly taken from big root tree (i.e. jec not reapplied) + fillVariableWithValue("htAK4",htAK4); // summing all jets with minimum pT cut and no jetid cut (jec not reapplied) + fillVariableWithValue("mhtAK4",mhtAK4); //summing all jets with minimum pT cut and no jetid cut (jec not reapplied) + fillVariableWithValue("mhtAK4Sig",mhtAK4Sig); // mhtAK4/htAK4 summing all jets with minimum pT cut and no jetid cut (jec not reapplied) + fillVariableWithValue("met",met); //directly taken from event + + + if( AK4jets.size() >=1 ){ + fillVariableWithValue( "IdTight_j1",idCaloJet[sortedJetIdx[0]]); + fillVariableWithValue( "pTAK4_j1", AK4jets[0].Pt()); + fillVariableWithValue( "etaAK4_j1", AK4jets[0].Eta()); + fillVariableWithValue( "phiAK4_j1", AK4jets[0].Phi()); + fillVariableWithValue( "jetJecAK4_j1", jecFactors[sortedJetIdx[0]] ); + fillVariableWithValue( "jetJecUncAK4_j1", jecUncertainty[sortedJetIdx[0]] ); + //jetID + fillVariableWithValue( "HadEnFrac_j1", jetHadfAK4->at(sortedJetIdx[0])); + fillVariableWithValue( "EmEnFrac_j1", jetEmfAK4->at(sortedJetIdx[0])); + fillVariableWithValue( "jetCSVAK4_j1", jetCSVAK4->at(sortedJetIdx[0]) ); + } + if( AK4jets.size() >=2 ){ + //cout << "IdTight_j2 : " << idTAK4->at(sortedJetIdx[1]) << endl << endl; + fillVariableWithValue( "IdTight_j2",idCaloJet[sortedJetIdx[1]]); + fillVariableWithValue( "pTAK4_j2", AK4jets[1].Pt() ); + fillVariableWithValue( "etaAK4_j2", AK4jets[1].Eta()); + fillVariableWithValue( "phiAK4_j2", AK4jets[1].Phi()); + fillVariableWithValue( "jetJecAK4_j2", jecFactors[sortedJetIdx[1]]); + fillVariableWithValue( "jetJecUncAK4_j2", jecUncertainty[sortedJetIdx[1]] ); + //jetID + fillVariableWithValue( "HadEnFrac_j2", jetHadfAK4->at(sortedJetIdx[1])); + fillVariableWithValue( "EmEnFrac_j2", jetEmfAK4->at(sortedJetIdx[1])); + fillVariableWithValue( "jetCSVAK4_j2", jetCSVAK4->at(sortedJetIdx[1]) ); + //dijet + fillVariableWithValue( "Dijet_MassAK4", MJJAK4) ; + fillVariableWithValue( "CosThetaStarAK4", TMath::TanH( (AK4jets[0].Eta()-AK4jets[1].Eta())/2 )); + fillVariableWithValue( "deltaETAjjAK4", DeltaEtaJJAK4 ) ; + fillVariableWithValue( "deltaPHIjjAK4", DeltaPhiJJAK4 ) ; + } + + if( widejets.size() >= 1 ){ + fillVariableWithValue( "pTWJ_j1", widejets[0].Pt() ); + fillVariableWithValue( "pTWJ_j1_noCorr", widejets_noCorr[0].Pt() ); + fillVariableWithValue( "etaWJ_j1", widejets[0].Eta()); + //no cuts on these variables, just to store in output + fillVariableWithValue( "massWJ_j1", widejets[0].M()); + fillVariableWithValue( "phiWJ_j1", widejets[0].Phi()); + fillVariableWithValue( "corr1_WJ1", corr1); + } + + if( widejets.size() >= 2 ){ + fillVariableWithValue( "pTWJ_j2", widejets[1].Pt() ); + fillVariableWithValue( "pTWJ_j2_noCorr", widejets_noCorr[1].Pt() ); + fillVariableWithValue( "etaWJ_j2", widejets[1].Eta()); + fillVariableWithValue( "deltaETAjj", DeltaEtaJJWide ) ; + fillVariableWithValue( "mjj", MJJWide ) ; + fillVariableWithValue( "mjj_noCorr", MJJWideNoCorr ) ; + fillVariableWithValue( "mjj_shiftJEC", MJJWide_shift ) ; + //no cuts on these variables, just to store in output + fillVariableWithValue( "massWJ_j2", widejets[1].M()); + fillVariableWithValue( "phiWJ_j2", widejets[1].Phi()); + //dijet + fillVariableWithValue( "CosThetaStarWJ", TMath::TanH( (widejets[0].Eta()-widejets[1].Eta())/2 )); + fillVariableWithValue( "deltaPHIjj", DeltaPhiJJWide ) ; + fillVariableWithValue( "corr2_WJ2", corr2); + //fillVariableWithValue( "Dijet_MassAK8", mjjAK8 ) ; + //fillVariableWithValue( "Dijet_MassC", mjjCA8 ) ; + // if(wdijet.M()<1){ + // std::cout << " INV MASS IS " << wdijet.M() << std::endl; + // std::cout << " Delta Eta IS " << DeltaEtaJJWide << " n is " << widejets.size() << std::endl; + // std::cout << " INV MASS FROM NTUPLE AK8 " << mjjAK8 << std::endl; + // //std::cout << " INV MASS FROM NTUPLE CA8 " << mjjCA8 << std::endl; + } + + //no cuts on these variables, just to store in output + // if(!isData) + // fillVariableWithValue("trueVtx",PileupInteractions->at(idx_InTimeBX)); + // else if(isData) + // fillVariableWithValue("trueVtx",999); + + // Trigger + //int NtriggerBits = triggerResult->size(); + if (isData) + { + fillVariableWithValue("passHLT_CaloJet40_CaloScouting_PFScouting",triggerResult->at(0));// CaloJet40_CaloScouting_PFScouting + fillVariableWithValue("passHLT_L1HTT_CaloScouting_PFScouting",triggerResult->at(2));// L1HTT_CaloScouting_PFScouting + fillVariableWithValue("passHLT_CaloScoutingHT250",triggerResult->at(4));// CaloScoutingHT250 + fillVariableWithValue("passHLT_PFScoutingHT450",triggerResult->at(8)); + fillVariableWithValue("passHLT_PFHT900",triggerResult->at(15)); + fillVariableWithValue("passHLT_PFHT800",triggerResult->at(16)); + fillVariableWithValue("passHLT_PFHT650MJJ950",triggerResult->at(25)); + fillVariableWithValue("passHLT_PFHT650MJJ900",triggerResult->at(26)); + fillVariableWithValue("passHLT_PFJET500",triggerResult->at(27)); + fillVariableWithValue("passHLT_PFJET450",triggerResult->at(28)); + fillVariableWithValue("passHLT_Mu45Eta2p1",triggerResult->at(32)); + fillVariableWithValue("passHLT_AK8PFHT700TriMass50",triggerResult->at(33)); + fillVariableWithValue("passHLT_AK8PFJet360TrimMass50",triggerResult->at(34)); + fillVariableWithValue("passHLT_CaloJet500NoJetID",triggerResult->at(35)); + + fillVariableWithValue("passL1T_HTT200",l1Result->at(0));// + fillVariableWithValue("passL1T_HTT240",l1Result->at(1));// + fillVariableWithValue("passL1T_HTT270",l1Result->at(2));// + fillVariableWithValue("passL1T_HTT280",l1Result->at(3));// + fillVariableWithValue("passL1T_HTT300",l1Result->at(4));// + fillVariableWithValue("passL1T_HTT320",l1Result->at(5));// + fillVariableWithValue("passL1T_ZeroBias",l1Result->at(6));// + } + + // Evaluate cuts (but do not apply them) + evaluateCuts(); + + // optional call to fill a skim with the full content of the input roottuple + //if( passedCut("nJetFinal") ) fillSkimTree(); + if( passedCut("PassJSON") + && passedCut("IdTight_j1") + && passedCut("IdTight_j2") + && passedCut("nJet") + && passedCut("pTWJ_j1") + && passedCut("etaWJ_j1") + && passedCut("pTWJ_j2") + && passedCut("etaWJ_j2") + && getVariableValue("deltaETAjj") < getPreCutValue1("DetaJJforTrig") ){ + + h_mjj_NoTrigger_1GeVbin -> Fill(MJJWide); + h_mjj_NoTrigger -> Fill(MJJWide); + + + // if( (getVariableValue("passHLT_L1HTT150_BtagSeq")||getVariableValue("passHLT_L1HTT150")) ) + // h_mjj_HLTpass_L1HTT150 -> Fill(MJJWide); + + if( getVariableValue("passHLT_CaloJet40_CaloScouting_PFScouting") ) + h_mjj_HLTpass_CaloJet40_CaloScouting_PFScouting -> Fill(MJJWide); + if( getVariableValue("passHLT_L1HTT_CaloScouting_PFScouting") ) + h_mjj_HLTpass_L1HTT_CaloScouting_PFScouting -> Fill(MJJWide); + if( getVariableValue("passHLT_CaloScoutingHT250") ) + h_mjj_HLTpass_CaloScoutingHT250 -> Fill(MJJWide); + } + + // optional call to fill a skim with a subset of the variables defined in the cutFile (use flag SAVE) + if( passedAllPreviousCuts("mjj") && passedCut("mjj") ) + { + fillReducedSkimTree(); + + // ===== Take a look at this ===== + // //Example on how to investigate quickly the data + // if(getVariableValue("mjj")>4000) + // { + // //xbfast creation and filling of histograms + // CreateAndFillUserTH1D("h_dphijj_mjjgt4000", 100, 0, 3.15, getVariableValue("deltaPHIjj")); + // CreateAndFillUserTH1D("h_htak4_mjjgt4000", 1000, 0, 10000, getVariableValue("HTAK4")); + // CreateAndFillUserTH1D("h_nvtx_mjjgt4000", 31, -0.5, 30.5, getVariableValue("nVtx")); + // } + + } + + // ===== Example of mjj spectrum after HLT selection ===== + // if( passedAllPreviousCuts("mjj") ) + // { + // if(getVariableValue("passHLT")>0) + // { + // //fast creation and filling of histograms + // CreateAndFillUserTH1D("h_mjj_passHLT", getHistoNBins("mjj"), getHistoMin("mjj"), getHistoMax("mjj"), getVariableValue("mjj")); + // } + // } + + // reject events that did not pass level 0 cuts + //if( !passedCut("0") ) continue; + // ...... + + // reject events that did not pass level 1 cuts + //if( !passedCut("1") ) continue; + // ...... + + // reject events that did not pass the full cut list + //if( !passedCut("all") ) continue; + // ...... + + // if( widejets.size() >= 2) { + // h_nJetFinal->Fill(widejets.size()); + // h_DijetMass->Fill(wdijet.M()); + // h_pT1stJet->Fill(widejets[0].Pt()); + // h_pT2ndJet->Fill(widejets[1].Pt()); + // h_eta1stJet->Fill(widejets[0].Eta()); + // h_eta2ndJet->Fill(widejets[1].Eta()); + // } + ////////////////////// User's code ends here /////////////////////// + + } // End loop over events + + //////////write histos + + h_mjj_NoTrigger_1GeVbin -> Write(); + h_mjj_NoTrigger -> Write(); + h_mjj_HLTpass_CaloJet40_CaloScouting_PFScouting -> Write(); + h_mjj_HLTpass_L1HTT_CaloScouting_PFScouting -> Write(); + h_mjj_HLTpass_CaloScoutingHT250 -> Write(); + + // h_nVtx->Write(); + // h_trueVtx->Write(); + // h_nJetFinal->Write(); + // h_pT1stJet->Write(); + // h_pT2ndJet->Write(); + // h_DijetMass->Write(); + // h_eta1stJet->Write(); + // h_eta2ndJet->Write(); + + // //pT of both jets, to be built using the histograms produced automatically by baseClass + // TH1F * h_pTJets = new TH1F ("h_pTJets","", getHistoNBins("pT1stJet"), getHistoMin("pT1stJet"), getHistoMax("pT1stJet")); + // h_pTJets->Add( & getHisto_noCuts_or_skim("pT1stJet") ); // all histos can be retrieved, see other getHisto_xxxx methods in baseClass.h + // h_pTJets->Add( & getHisto_noCuts_or_skim("pT2ndJet") ); + // //one could also do: *h_pTJets = getHisto_noCuts_or_skim("pT1stJet") + getHisto_noCuts_or_skim("pT2ndJet"); + // h_pTJets->Write(); + // //one could also do: const TH1F& h = getHisto_noCuts_or_skim// and use h + + std::cout << "analysisClass::Loop() ends" <