|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
|
/**
*/
#include "Calibration/EcalCalibAlgos/interface/IMACalibBlock.h"
#include "Calibration/Tools/interface/BlockSolver.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/isFinite.h"
#include "TH1F.h"
#include "TFile.h"
#include <cassert>
// -----------------------------------------------------
IMACalibBlock::IMACalibBlock(const int numberOfElements)
: VEcalCalibBlock(numberOfElements), m_kaliVector(m_numberOfElements), m_kaliMatrix(evalX2Size()) {
reset();
}
// -----------------------------------------------------
IMACalibBlock::~IMACalibBlock() {}
// -----------------------------------------------------
void IMACalibBlock::Fill(std::map<int, double>::const_iterator MapBegin,
std::map<int, double>::const_iterator MapEnd,
double pTk,
double pSubtract,
double sigma) {
// std::cerr<<"\n\nfaccio il fill!\n";
double inverror = 1. / sigma;
//LP fist loop over energies
for (std::map<int, double>::const_iterator itMap1 = MapBegin; itMap1 != MapEnd; ++itMap1) {
// std::cerr<<"numero "<<itMap1->first<<" vale "<<itMap1->second<<"\t";
for (std::map<int, double>::const_iterator itMap2 = itMap1; itMap2 != MapEnd; ++itMap2) {
//LP calculate the chi square value
double dummy = itMap1->second * itMap2->second;
dummy *= inverror;
//LP fill the calib matrix
m_kaliMatrix.at(itMap1->first * m_numberOfElements + itMap2->first) += dummy;
} //LP second loop over xtals
//LP calculate the vector value
double dummy = pTk;
dummy -= pSubtract; //LP borders
dummy *= itMap1->second;
dummy *= inverror;
//LP fill the calib vector
m_kaliVector.at(itMap1->first) += dummy;
} //LP first loop over energies
return;
}
//------------------------------------------------------------
void IMACalibBlock::complete() {
int bef;
int aft;
for (unsigned int i = 0; i < m_numberOfElements; ++i) {
for (unsigned int j = i + 1; j < m_numberOfElements; ++j) {
bef = (i * m_numberOfElements + j);
aft = (j * m_numberOfElements + i);
m_kaliMatrix.at(aft) = m_kaliMatrix.at(bef);
} //LP second loop over xtals
} //LP first loop over xtals
return;
}
// ------------------------------------------------------------
int IMACalibBlock::solve(int usingBlockSolver, double min, double max) {
complete();
// TH1F vettore ("vettore","vettore",10,-0.1,9.9);
// TH1F matrice ("matrice","matrice",100,-0.1,99.9);
CLHEP::HepMatrix kaliMatrix(m_numberOfElements, m_numberOfElements);
// for (std::vector<double>::iterator it = m_kaliVector.begin();
// it!= m_kaliVector.end();++it)
// vettore.Fill(it-m_kaliVector.begin(),*it);
riempiMtr(m_kaliMatrix, kaliMatrix);
// for (std::vector<double>::iterator it = m_kaliMatrix.begin();
// it!= m_kaliMatrix.end();++it)
// matrice.Fill(it-m_kaliMatrix.begin(),*it);
// TFile f ("fileInteressante.root","RECREATE");
// vettore.Write();
// matrice.Write();
CLHEP::HepVector kaliVector(m_numberOfElements);
riempiVtr(m_kaliVector, kaliVector);
//PG linear system solution
CLHEP::HepVector result = CLHEP::solve(kaliMatrix, kaliVector);
for (int i = 0; i < kaliVector.num_row(); ++i)
if (result.normsq() < min * kaliVector.num_row() || result.normsq() > max * kaliVector.num_row()) {
if (usingBlockSolver) {
edm::LogWarning("IML") << "using blocSlover " << std::endl;
BlockSolver()(kaliMatrix, kaliVector, result);
} else {
edm::LogWarning("IML") << "coeff out of range " << std::endl;
for (int i = 0; i < kaliVector.num_row(); ++i)
result[i] = 1.;
}
}
fillMap(result);
return 0;
}
//------------------------------------------------------------
inline int IMACalibBlock::evalX2Size() { return m_numberOfElements * m_numberOfElements; }
// ------------------------------------------------------------
void IMACalibBlock::riempiMtr(const std::vector<double>& piena, CLHEP::HepMatrix& vuota) {
unsigned int max = m_numberOfElements;
assert(piena.size() == max * max);
assert(vuota.num_row() == int(max));
assert(vuota.num_col() == int(max));
for (unsigned int i = 0; i < max; ++i)
for (unsigned int j = 0; j < max; ++j)
if (edm::isNotFinite(piena[i * max + j]))
vuota[i][j] = 0.;
else
vuota[i][j] = piena[i * max + j];
return;
}
// ------------------------------------------------------------
void IMACalibBlock::riempiVtr(const std::vector<double>& pieno, CLHEP::HepVector& vuoto) {
int max = m_numberOfElements;
assert(vuoto.num_row() == max);
for (int i = 0; i < max; ++i)
if (edm::isNotFinite(pieno[i]))
vuoto[i] = 0.;
else
vuoto[i] = pieno[i];
return;
}
// ------------------------------------------------------------
void IMACalibBlock::reset() {
for (std::vector<double>::iterator vecIt = m_kaliVector.begin(); vecIt != m_kaliVector.end(); ++vecIt) {
*vecIt = 0.;
}
for (std::vector<double>::iterator vecIt = m_kaliMatrix.begin(); vecIt != m_kaliMatrix.end(); ++vecIt) {
*vecIt = 0.;
}
}
// ------------------------------------------------------------
//LP sta provando a fare i seguenti due metodi come locali. Speriamo di non fare stronzate.
void IMACalibBlock::fillMap(const CLHEP::HepVector& result) {
for (unsigned int i = 0; i < m_numberOfElements; ++i) {
m_coefficients[i] = result[i];
}
return;
}
|