PFClusterTools/test/ProducePFCalibrationObject.cc

0001
0002 /*
0003  *  See header file for a description of this class.
0004  *
0005  *  Original author G. Cerminara - CERN
0006  *  Modified by  spandey/bkansal
0007  */
0008
0009 #include "ProducePFCalibrationObject.h"
0010
0011 #include <iostream>
0012 #include <vector>
0013 #include "FWCore/ServiceRegistry/interface/Service.h"
0014 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0015
0016 #include "FWCore/Framework/interface/ESHandle.h"
0017 #include "FWCore/Framework/interface/EventSetup.h"
0018
0019 #include "TF1.h"
0020
0021 using namespace std;
0022 using namespace edm;
0023
0024 ProducePFCalibrationObject::ProducePFCalibrationObject(const edm::ParameterSet& pSet)
0025     : read(pSet.getUntrackedParameter<bool>("read")),
0026       write(pSet.getUntrackedParameter<bool>("write")),
0027       fToWrite(pSet.getParameter<vector<ParameterSet> >("toWrite")),
0028       fToRead(pSet.getUntrackedParameter<vector<string> >("toRead")),
0029       perfToken(esConsumes<edm::Transition::BeginRun>()) {}
0030
0031 ProducePFCalibrationObject::~ProducePFCalibrationObject() {}
0032
0033 void ProducePFCalibrationObject::beginRun(const edm::Run& run, const edm::EventSetup& eSetup) {
0034   cout << "[ProducePFCalibrationObject] beginJob" << endl;
0035
0036   string record = "PFCalibrationRcd";
0037
0038   static map<string, PerformanceResult::ResultType> functType;
0039   functType["PFfa_BARREL"] = PerformanceResult::PFfa_BARREL;
0040   functType["PFfa_ENDCAP"] = PerformanceResult::PFfa_ENDCAP;
0041   functType["PFfb_BARREL"] = PerformanceResult::PFfb_BARREL;
0042   functType["PFfb_ENDCAP"] = PerformanceResult::PFfb_ENDCAP;
0043   functType["PFfc_BARREL"] = PerformanceResult::PFfc_BARREL;
0044   functType["PFfc_ENDCAP"] = PerformanceResult::PFfc_ENDCAP;
0045   //New added functions for H/EH hadrons
0046   functType["PFfaEta_BARRELH"] = PerformanceResult::PFfaEta_BARRELH;
0047   functType["PFfaEta_ENDCAPH"] = PerformanceResult::PFfaEta_ENDCAPH;
0048   functType["PFfbEta_BARRELH"] = PerformanceResult::PFfbEta_BARRELH;
0049   functType["PFfbEta_ENDCAPH"] = PerformanceResult::PFfbEta_ENDCAPH;
0050   functType["PFfaEta_BARRELEH"] = PerformanceResult::PFfaEta_BARRELEH;
0051   functType["PFfaEta_ENDCAPEH"] = PerformanceResult::PFfaEta_ENDCAPEH;
0052   functType["PFfbEta_BARRELEH"] = PerformanceResult::PFfbEta_BARRELEH;
0053   functType["PFfbEta_ENDCAPEH"] = PerformanceResult::PFfbEta_ENDCAPEH;
0054   //Left older functions untouched for backward compatibility
0055   functType["PFfaEta_BARREL"] = PerformanceResult::PFfaEta_BARREL;
0056   functType["PFfaEta_ENDCAP"] = PerformanceResult::PFfaEta_ENDCAP;
0057   functType["PFfbEta_BARREL"] = PerformanceResult::PFfbEta_BARREL;
0058   functType["PFfbEta_ENDCAP"] = PerformanceResult::PFfbEta_ENDCAP;
0059
0060   functType["PFfcEta_BARRELH"] = PerformanceResult::PFfcEta_BARRELH;
0061   functType["PFfcEta_ENDCAPH"] = PerformanceResult::PFfcEta_ENDCAPH;
0062   functType["PFfdEta_ENDCAPH"] = PerformanceResult::PFfdEta_ENDCAPH;
0063   functType["PFfcEta_BARRELEH"] = PerformanceResult::PFfcEta_BARRELEH;
0064   functType["PFfcEta_ENDCAPEH"] = PerformanceResult::PFfcEta_ENDCAPEH;
0065   functType["PFfdEta_ENDCAPEH"] = PerformanceResult::PFfdEta_ENDCAPEH;
0066
0067   // ---------------------------------------------------------------------------------
0068   // Write the payload
0069
0070   if (write) {
0071     vector<pair<float, float> > limitsToWrite;
0072     vector<string> formulasToWrite;
0073     vector<PerformanceResult::ResultType> resToWrite;
0074     vector<BinningVariables::BinningVariablesType> binsToWrite;
0075     binsToWrite.push_back(BinningVariables::JetEt);
0076
0077     // loop over all the pSets for the TF1 that we want to write to DB
0078     for (vector<ParameterSet>::const_iterator fSetup = fToWrite.begin(); fSetup != fToWrite.end(); ++fSetup) {
0079       // read from cfg
0080       string fType = (*fSetup).getUntrackedParameter<string>("fType");
0081       //FIXME: should check that the give type exists
0082       string formula = (*fSetup).getUntrackedParameter<string>("formula");
0083       pair<float, float> limits = make_pair((*fSetup).getUntrackedParameter<vector<double> >("limits")[0],
0084                                             (*fSetup).getUntrackedParameter<vector<double> >("limits")[1]);
0085       vector<double> parameters = (*fSetup).getUntrackedParameter<vector<double> >("parameters");
0086
0087       TF1* function = new TF1(fType.c_str(), formula.c_str(), limits.first, limits.second);
0088       for (unsigned int index = 0; index != parameters.size(); ++index) {
0089         function->SetParameter(index, parameters[index]);
0090       }
0091
0092       // write them in the containers for the storage
0093       limitsToWrite.push_back(limits);
0094       formulasToWrite.push_back(string(function->GetExpFormula("p").Data()));
0095       resToWrite.push_back(functType[fType]);
0096     }
0097
0098     // create the actual object storing the functions
0099     PhysicsTFormulaPayload allTFormulas(limitsToWrite, formulasToWrite);
0100
0101     // put them in the container
0102     PerformancePayloadFromTFormula pfCalibrationFormulas(resToWrite, binsToWrite, allTFormulas);
0103
0104     // actually write to DB
0105     edm::Service<cond::service::PoolDBOutputService> dbOut;
0106     if (dbOut.isAvailable()) {
0107       dbOut->writeOneIOV(pfCalibrationFormulas, 1, record);
0108     }
0109   }
0110
0111   if (read) {
0112     // ---------------------------------------------------------------------------------
0113     // Read the objects
0114     const PerformancePayloadFromTFormula* pfCalibrations =
0115         static_cast<const PerformancePayloadFromTFormula*>(&eSetup.getData(perfToken));
0116
0117     for (vector<string>::const_iterator name = fToRead.begin(); name != fToRead.end(); ++name) {
0118       cout << "Function: " << *name << endl;
0119       PerformanceResult::ResultType fType = functType[*name];
0120       pfCalibrations->printFormula(fType);
0121
0122       // evaluate it @ 10 GeV
0123       float energy = 10.;
0124
0125       BinningPointByMap point;
0126       point.insert(BinningVariables::JetEt, energy);
0127
0128       if (pfCalibrations->isInPayload(fType, point)) {
0129         float value = pfCalibrations->getResult(fType, point);
0130         cout << "   Energy before:: " << energy << " after: " << value << endl;
0131       } else
0132         cout << "outside limits!" << endl;
0133     }
0134   }
0135   ///  if(pfCalibrationFormulas->isInPayload(etaBin, point)) {
0136   //     float value = pfCalibrationFormulas->getResult(etaBin, point);
0137   //     cout << "t: " << t << " eta: " << eta << " CalibObj: " <<
0138   //       value << endl;
0139   //   } else cout <<  "INVALID result!!!" << endl;
0140
0141   //   TF1* faBarrel = new TF1("faBarrel","[0]+([1]+[2]/sqrt(x))*exp(-x/[3])-[4]*exp(-x*x/[5])",1.,1000.);
0142
0143   //   faBarrel->SetParameter(0,1.10772);
0144   //   faBarrel->SetParameter(1,0.186273);
0145   //   faBarrel->SetParameter(2,-0.47812);
0146   //   faBarrel->SetParameter(3,62.5754);
0147   //   faBarrel->SetParameter(4,1.31965);
0148   //   faBarrel->SetParameter(5,35.2559);
0149
0150   //   // faBarrel->GetExpFormula("p").Data()
0151   //   TF1 * faEndcap = new TF1("faEndcap","[0]+([1]+[2]/sqrt(x))*exp(-x/[3])-[4]*exp(-x*x/[5])",1.,1000.);
0152   //   faEndcap->SetParameter(0,1.0877);
0153   //   faEndcap->SetParameter(1,0.28939);
0154   //   faEndcap->SetParameter(2,-0.57635);
0155   //   faEndcap->SetParameter(3,86.5501);
0156   //   faEndcap->SetParameter(4,1.02296);
0157   //   faEndcap->SetParameter(5,64.0116);
0158
0159   //   std::vector< std::pair<float, float> > limits;
0160   //   limits.push_back(make_pair(0,99999999));
0161   //   //  limits.push_back(make_pair(0,1.48 ));
0162
0163   //   limits.push_back(make_pair(0,99999999));
0164   //   //  limits.push_back(make_pair(1.48,5));
0165
0166   //   std::vector<std::string> formulas;
0167   //   formulas.push_back(string(faBarrel->GetExpFormula("p").Data()));
0168   //   formulas.push_back(string(faEndcap->GetExpFormula("p").Data()));
0169
0170   //   std::vector<PerformanceResult::ResultType> res;
0171   //   res.push_back(PerformanceResult::PFfa_BARREL);
0172   //   res.push_back(PerformanceResult::PFfa_ENDCAP);
0173
0174   //   std::vector<BinningVariables::BinningVariablesType> bin;
0175   //   bin.push_back(BinningVariables::JetEt);
0176   //   //  bin.push_back(BinningVariables::JetAbsEta);
0177
0178   //   bin.push_back(BinningVariables::JetEt);
0179   //   //  bin.push_back(BinningVariables::JetAbsEta);
0180
0181   //   PhysicsTFormulaPayload ppl(limits, formulas);
0182   // //   PerformanceResult::PFfa_BARREL
0183
0184   //   PerformancePayloadFromTFormula * pfCalibrationFormulas =
0185   //     new PerformancePayloadFromTFormula(res,
0186   //                       bin,
0187   //                       ppl);
0188
0189   //   double t = 10.;
0190   //   double eta = 2.1;
0191
0192   //   BinningPointByMap point;
0193   //   point.insert(BinningVariables::JetEt, t);
0194   //   //point.insert(BinningVariables::JetAbsEta, eta);
0195
0196   //   PerformanceResult::ResultType etaBin;
0197
0198   //   if(fabs(eta) < 1.48 ) {
0199   //     // this is the barrel
0200   //     etaBin = PerformanceResult::PFfa_BARREL;
0201   //     cout << " f_barrel(a): " << faBarrel->Eval(t) << endl;
0202
0203   //   } else {
0204   //     // this is the endcap
0205   //     etaBin = PerformanceResult::PFfa_BARREL;
0206   //     cout << " f_endcap(a): " << faEndcap->Eval(t) << endl;
0207
0208   // //     if(pfCalibrationFormulas->isInPayload(PerformanceResult::PFfa_ENDCAP, point)){
0209   // //       float value = pfCalibrationFormulas->getResult(PerformanceResult::PFfa_ENDCAP, point);
0210   // //       cout << "t: " << t << " eta: " << eta << " f_endcap(a): " << faEndcap->Eval(t) << " CalibObj: " <<
0211   // //     value << endl;
0212   // //     } else cout <<  "INVALID result!!!" << endl;
0213   //   }
0214
0215   //   if(pfCalibrationFormulas->isInPayload(etaBin, point)) {
0216   //     float value = pfCalibrationFormulas->getResult(etaBin, point);
0217   //     cout << "t: " << t << " eta: " << eta << " CalibObj: " <<
0218   //       value << endl;
0219   //   } else cout <<  "INVALID result!!!" << endl;
0220 }
0221
0222 #include "FWCore/Framework/interface/MakerMacros.h"
0223
0224 DEFINE_FWK_MODULE(ProducePFCalibrationObject);