|
|
|
#include <HcalCholeskyDecomp.h>
Public Member Functions | |
| HcalCholeskyDecomp (const edm::ParameterSet &) | |
| ~HcalCholeskyDecomp () | |
Private Member Functions | |
| virtual void | analyze (const edm::Event &, const edm::EventSetup &) |
| virtual void | beginJob () |
| virtual void | endJob () |
Private Attributes | |
| std::string | outfile |
Definition at line 42 of file HcalCholeskyDecomp.h.
| HcalCholeskyDecomp::HcalCholeskyDecomp | ( | const edm::ParameterSet & | iConfig | ) | [explicit] |
Definition at line 6 of file HcalCholeskyDecomp.cc.
References edm::ParameterSet::getUntrackedParameter(), and outfile.
{
outfile = iConfig.getUntrackedParameter<std::string>("outFile","CholeskyMatrices.txt");
}
| HcalCholeskyDecomp::~HcalCholeskyDecomp | ( | ) |
Definition at line 12 of file HcalCholeskyDecomp.cc.
{
}
| void HcalCholeskyDecomp::analyze | ( | const edm::Event & | iEvent, |
| const edm::EventSetup & | iSetup | ||
| ) | [private, virtual] |
Implements edm::EDAnalyzer.
Definition at line 17 of file HcalCholeskyDecomp.cc.
References HcalElectronicsMap::allPrecisionId(), trackerHits::c, prof2calltree::count, gather_cfg::cout, CastorDbASCIIIO::dumpObject(), spr::find(), edm::EventSetup::get(), HcalCovarianceMatrix::getValue(), i, j, gen::k, m, outfile, HcalCholeskyMatrix::setValue(), and mathSSE::sqrt().
{
using namespace edm;
using namespace std;
edm::ESHandle<HcalCovarianceMatrices> refCov;
iSetup.get<HcalCovarianceMatricesRcd>().get("reference",refCov);
const HcalCovarianceMatrices* myCov = refCov.product(); //Fill emap from database
double sig[4][10][10];
double c[4][10][10], cikik[4], cikjk[4], ct[4][10][10];
HcalCholeskyMatrices * outMatrices = new HcalCholeskyMatrices();
std::vector<DetId> listChan = myCov->getAllChannels();
std::vector<DetId>::iterator cell;
int HBcount = 0;
int HEcount = 0;
int HFcount = 0;
int HOcount = 0;
double HBmatrix[4][10][10];
double HEmatrix[4][10][10];
double HFmatrix[4][10][10];
double HOmatrix[4][10][10];
double tempmatrix[4][10][10] ;
for(int m= 0; m != 4; m++){
for(int i = 0; i != 10; i++){
for(int j = 0; j!=10; j++){
HBmatrix[m][i][j] = 0;
HEmatrix[m][i][j] = 0;
HFmatrix[m][i][j] = 0;
HOmatrix[m][i][j] = 0;
}
}
}
for (std::vector<DetId>::iterator it = listChan.begin(); it != listChan.end(); it++)
{
for(int m= 0; m != 4; m++){
for(int i = 0; i != 10; i++){
for(int j = 0; j!=10; j++){
sig[m][i][j] = 0;
c[m][i][j] = 0;
ct[m][i][j] = 0;
cikik[m] =0;
cikjk[m] = 0;
tempmatrix[m][i][j] = 0;
}
}
}
const HcalCovarianceMatrix * CMatrix = myCov->getValues(*it);
HcalDetId hcalid(it->rawId());
for(int m = 0; m != 4; m++){
for(int i = 0; i != 10; i++){
for(int j = 0; j != 10; j++) {sig[m][i][j] = CMatrix->getValue(m,i,j);}
}
}
//Method to check 4x10 storage
// for(int m = 0; m != 4; m++){
// for(int i = 0; i != 6; i++){
// for(int j=(i+4); j!=10; j++) {sig[m][i][j] = 0;}
// }
// }
for(int m = 0; m != 4; m++){
for(int i = 0; i != 10; i++){
for(int j = i; j != 10; j++){ sig[m][j][i] = sig[m][i][j];}
}
}
//.......................Cholesky Decomposition.............................
//Step 1a
for(int m = 0; m!=4; m++)
{
c[m][0][0]=sqrt(sig[m][0][0]);
}
for(int m = 0; m!=4; m++)
{
for(int i = 1; i != 10; i++)
{
c[m][i][0] = (sig[m][0][i] / c[m][0][0]);
}
}
//Chelesky Matrices
for(int m = 0; m!=4; m++)
{
c[m][1][1] = sqrt(sig[m][1][1] - (c[m][1][0]*c[m][1][0])); // i) step 2a of chelesky decomp
if (((sig[m][1][1]) - (c[m][1][0]*c[m][1][0]))<=0) continue;
for(int i = 2; i !=10; i++)
{
cikik[m] = 0;
for(int j=1; j!= i; j++)
{
cikjk[m] = 0;
for(int k=0; k != j; k++)
{
cikjk[m] += (c[m][i][k]*c[m][j][k]); // ii) step 2a of chelesky decomp
}
c[m][i][j] = (sig[m][i][j] - cikjk[m])/c[m][j][j]; // step 3a of chelesky decomp
}//end of j loop
for( int k = 0; k != i; k++)
{
cikik[m] += (c[m][i][k]*c[m][i][k]);
}
double test = (sig[m][i][i] - cikik[m]);
if(test > 0 )
{
c[m][i][i] = sqrt(test);
}
else
{
c[m][i][i] = 1000;
}
}//end of i loop
/*
//Cholesky Matrix for rechit (2-5 ts)
for (int i = 0; i!=2; i++)
{
for (int j=0; j!=10; j++)
{
c[m][i][j] = 0;
}
}
for (int i = 6; i!=10; i++)
{
for (int j=0; j!=10; j++)
{
c[m][i][j] = 0;
}
}
*/
}//end of m loop
HcalCholeskyMatrix * item = new HcalCholeskyMatrix(it->rawId());
for(int m = 0; m != 4; m++){
for(int i = 0; i != 10; i++){
for(int j = i; j != 10; j++){
item->setValue(m,j,i, c[m][j][i]);
tempmatrix[m][j][i] = c[m][j][i];
}//sig[m][j][i]
}
}
if(hcalid.subdet()==1)
{
for(int m = 0; m != 4; m++)
for(int i = 0; i != 10; i++)
for(int j = 0; j != 10; j++)
HBmatrix[m][i][j] += tempmatrix[m][i][j];
HBcount++;
}
if(hcalid.subdet()==2)
{
for(int m = 0; m != 4; m++)
for(int i = 0; i != 10; i++)
for(int j = 0; j != 10; j++)
HEmatrix[m][i][j] += tempmatrix[m][i][j];
HEcount++;
}
if(hcalid.subdet()==3)
{
for(int m = 0; m != 4; m++)
for(int i = 0; i != 10; i++)
for(int j = 0; j != 10; j++)
HOmatrix[m][i][j] += tempmatrix[m][i][j];
HOcount++;
}
if(hcalid.subdet()==4)
{
for(int m = 0; m != 4; m++)
for(int i = 0; i != 10; i++)
for(int j = 0; j != 10; j++)
HFmatrix[m][i][j] += tempmatrix[m][i][j];
HFcount++;
}
outMatrices->addValues(*item);
}
for(int m = 0; m != 4; m++)
{
for(int i = 0; i != 10; i++)
{
for(int j = 0; j != 10; j++)
{
HBmatrix [m][i][j] /= HBcount;
HEmatrix [m][i][j] /= HEcount;
HFmatrix [m][i][j] /= HFcount;
HOmatrix [m][i][j] /= HOcount;
}
}
}
std::vector<DetId> listResult = outMatrices->getAllChannels();
int n_avg = 0;
edm::ESHandle<HcalElectronicsMap> refEMap;
iSetup.get<HcalElectronicsMapRcd>().get(refEMap);
// iSetup.get<HcalElectronicsMapRcd>().get("reference",refEMap);
const HcalElectronicsMap* myRefEMap = refEMap.product();
std::vector<HcalGenericDetId> listEMap = myRefEMap->allPrecisionId();
for (std::vector<HcalGenericDetId>::const_iterator it = listEMap.begin(); it != listEMap.end(); it++)
{
DetId mydetid = DetId(it->rawId());
if (std::find(listResult.begin(), listResult.end(), mydetid ) == listResult.end() )
{
// std::cout << "Cholesky Matrix not found for cell " << HcalGenericDetId(it->rawId());
if(it->isHcalDetId())
{
std::cout << "Cholesky Matrix not found for cell " << HcalGenericDetId(it->rawId());
HcalDetId hcalid2(it->rawId());
HcalCholeskyMatrix * item = new HcalCholeskyMatrix(it->rawId());
for(int m = 0; m != 4; m++)
{
for(int i = 0; i != 10; i++)
{
for(int j = 0; j != 10; j++)
{
if(j <= i){
if(hcalid2.subdet()==1) item->setValue(m,i,j,HBmatrix [m][i][j]);
if(hcalid2.subdet()==2) item->setValue(m,i,j,HEmatrix [m][i][j]);
if(hcalid2.subdet()==3) item->setValue(m,i,j,HFmatrix [m][i][j]);
if(hcalid2.subdet()==4) item->setValue(m,i,j,HOmatrix [m][i][j]);
}else{
if(hcalid2.subdet()==1) item->setValue(m,i,j,0);
if(hcalid2.subdet()==2) item->setValue(m,i,j,0);
if(hcalid2.subdet()==3) item->setValue(m,i,j,0);
if(hcalid2.subdet()==4) item->setValue(m,i,j,0);
}
}
}
}
outMatrices->addValues(*item);
std::cout << "... Average Matrix used\n";
n_avg++;
}
}
}
std::vector<DetId> checkResult = outMatrices->getAllChannels();
int count = 0;
for(std::vector<DetId>::const_iterator it = checkResult.begin(); it != checkResult.end(); it++) count ++;
std::cout << "There are " << count << " channels with Cholesky matrices.\n" << "Used " << n_avg << " average values.\n";
ofstream cmatout(outfile.c_str());
HcalDbASCIIIO::dumpObject(cmatout, *outMatrices);
}
| void HcalCholeskyDecomp::beginJob | ( | void | ) | [private, virtual] |
| void HcalCholeskyDecomp::endJob | ( | void | ) | [private, virtual] |
std::string HcalCholeskyDecomp::outfile [private] |
Definition at line 52 of file HcalCholeskyDecomp.h.
Referenced by analyze(), and HcalCholeskyDecomp().
1.7.3