00001
00002 #include <memory>
00003 #include <string>
00004 #include <iostream>
00005
00006
00007 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00008 #include "Calibration/HcalCalibAlgos/interface/HcalConstantsASCIIWriter.h"
00009 #include "DataFormats/Provenance/interface/Provenance.h"
00010 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00011 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00012 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
00013 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
00014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00015 #include "FWCore/ParameterSet/interface/FileInPath.h"
00016 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
00017 #include "DataFormats/DetId/interface/DetId.h"
00018 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
00019
00020 #include "CondFormats/HcalObjects/interface/HcalRespCorrs.h"
00021 #include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"
00022 #include "HepMC/GenParticle.h"
00023 #include "HepMC/GenVertex.h"
00024 #include "CalibCalorimetry/HcalAlgos/interface/HcalDbASCIIIO.h"
00025 #include "CondFormats/DataRecord/interface/HcalRespCorrsRcd.h"
00026 #include <fstream>
00027 #include <sstream>
00028 #include <map>
00029 #include <vector>
00030
00031 using namespace std;
00032 using namespace reco;
00033
00034
00035
00036 namespace cms{
00037 HcalConstantsASCIIWriter::HcalConstantsASCIIWriter(const edm::ParameterSet& iConfig)
00038 {
00039
00040 file_input="Calibration/HcalCalibAlgos/data/"+iConfig.getParameter <std::string> ("fileInput")+".txt";
00041 file_output="Calibration/HcalCalibAlgos/data/"+iConfig.getParameter <std::string> ("fileOutput")+".txt";
00042 }
00043
00044 HcalConstantsASCIIWriter::~HcalConstantsASCIIWriter()
00045 {
00046
00047
00048
00049
00050 }
00051
00052 void HcalConstantsASCIIWriter::beginJob( const edm::EventSetup& iSetup)
00053 {
00054 edm::FileInPath f1(file_output);
00055 string fDataFile = f1.fullPath();
00056
00057 myout_hcal = new ofstream(fDataFile.c_str());
00058 if(!myout_hcal) cout << " Output file not open!!! "<<endl;
00059
00060 }
00061
00062 void HcalConstantsASCIIWriter::endJob()
00063 {
00064 delete myout_hcal;
00065 }
00066
00067
00068
00069
00070
00071
00072
00073 void
00074 HcalConstantsASCIIWriter::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00075 {
00076
00077 std::cout<<" Start HcalConstantsASCIIWriter::analyze "<<std::endl;
00078
00079 edm::ESHandle<HcalRespCorrs> r;
00080 iSetup.get<HcalRespCorrsRcd>().get(r);
00081 HcalRespCorrs* oldRespCorrs = new HcalRespCorrs(*r.product());
00082
00083
00084 edm::ESHandle<CaloGeometry> pG;
00085 iSetup.get<CaloGeometryRecord>().get(pG);
00086 const CaloGeometry* geo = pG.product();
00087
00088
00089 std::vector<DetId> did = geo->getValidDetIds();
00090
00091
00092 map<HcalDetId,float> corrold;
00093
00094
00095 int mydet,mysubd,depth,ieta,iphi;
00096 float coradd,corerr;
00097
00098 vector<HcalDetId> theVector;
00099 for(std::vector<DetId>::iterator i = did.begin(); i != did.end(); i++)
00100 {
00101 if( (*i).det() == DetId::Hcal ) {
00102 HcalDetId hid = HcalDetId(*i);
00103 theVector.push_back(hid);
00104 corrold[hid] = (oldRespCorrs->getValues(*i))->getValue();
00105 std::cout<<" Old calibration "<<hid.depth()<<" "<<hid.ieta()<<" "<<hid.iphi()<<std::endl;
00106 }
00107 }
00108
00109 std::cout<<" Get old calibration "<<std::endl;
00110
00111
00112 edm::FileInPath f1(file_input);
00113 string fDataFile = f1.fullPath();
00114
00115 std::ifstream in( fDataFile.c_str() );
00116 string line;
00117
00118 double corrnew_p[5][5][45][75];
00119 double corrnew_m[5][5][45][75];
00120 std::cout<<" Start to read txt file "<<fDataFile.c_str()<<std::endl;
00121 while( std::getline( in, line)){
00122
00123
00124
00125 if(!line.size() || line[0]=='#') continue;
00126 istringstream linestream(line);
00127 double par;
00128 int type;
00129 linestream>>mysubd>>depth>>ieta>>iphi>>coradd>>corerr;
00130
00131
00132
00133
00134 int ietak = ieta;
00135 if(ieta<0) ietak = -1*ieta;
00136 if(ieta>0) corrnew_p[mysubd][depth][ietak][iphi] = coradd;
00137 if(ieta<0) corrnew_m[mysubd][depth][ietak][iphi] = coradd;
00138 std::cout<<" Try to initialize mysubd "<<mysubd<<" depth "<<depth<<" ieta "<<ieta<<" "<<ietak<<" iphi "<<iphi<<" "<<coradd<<
00139 std::endl;
00140 }
00141
00142 HcalRespCorrs* mycorrections = new HcalRespCorrs();
00143
00144 for(vector<HcalDetId>::iterator it = theVector.begin(); it != theVector.end(); it++)
00145 {
00146 float cc1 = (*corrold.find(*it)).second;
00147
00148 float cc2 = 0.;
00149 int ietak = (*it).ieta();
00150
00151 if((*it).ieta()<0) ietak=-1*(*it).ieta();
00152
00153 if((*it).ieta()>0) cc2 = corrnew_p[(*it).subdet()][(*it).depth()][ietak][(*it).iphi()];
00154 if((*it).ieta()<0) cc2 = corrnew_m[(*it).subdet()][(*it).depth()][ietak][(*it).iphi()];
00155
00156 float cc = cc1*cc2;
00157 std::cout<<" Multiply "<<(*it).subdet()<<" "<<(*it).depth()<<" "<<(*it).ieta()<<" "<<ietak<<" "<<(*it).iphi()<<" "<<(*it).rawId()<<" "<<cc1<<" "<<cc2<<std::endl;
00158
00159
00160 HcalRespCorr item ((*it).rawId(), cc);
00161 bool rr = mycorrections->addValues(item);
00162 }
00163
00164 HcalRespCorrs mycc = *mycorrections;
00165 HcalDbASCIIIO::dumpObject (*myout_hcal, mycc);
00166
00167 }
00168 }
00169
00170
00171