8 #include "TVirtualFitter.h"
10 using namespace pftools;
13 name_(name), barrelLimit_(1.4), transitionLimit_(1.4), endcapLimit_(5.0) {
23 for (std::map<SpaceVoxelPtr, CalibratorPtr>::iterator it =
31 const double etaSeg,
const double phiSeg,
const double energySeg) {
33 <<
": this method has not yet been implemented!\n";
42 <<
": creating default calibration schema.\n";
169 const unsigned nEta,
const double etaMin,
const double etaMax,
170 const unsigned nPhi,
const double phiMin,
const double phiMax,
171 const unsigned nEnergy,
const double energyMin,
const double energyMax)
175 if (nEta == 0|| nPhi ==0|| nEnergy == 0) {
177 me(
"Can't create calibrators with zero values for nEta, nPhi or nEnergy!");
181 double etaSeg = (etaMax - etaMin) / nEta;
182 double phiSeg = (phiMax - phiMin) / nPhi;
183 double energySeg = (energyMax - energyMin) / nEnergy;
185 double eta1, eta2, phi1, phi2, energy1, energy2;
186 for (
unsigned k(0);
k < nEta; ++
k) {
187 for (
unsigned l(0);
l < nPhi; ++
l) {
188 for (
unsigned m(0);
m < nEnergy; ++
m) {
189 eta1 = etaMin +
k * etaSeg;
190 eta2 = eta1 + etaSeg;
192 phi1 = phiMin +
l * phiSeg;
193 phi2 = phi1 + phiSeg;
195 energy1 = energyMin +
m * energySeg;
196 energy2 = energy1 + energySeg;
198 myKnownSpaceVoxels.push_back(sv);
200 myAddressBook[sv] =
c;
204 unsigned nCalibrators = nEta * nPhi * nEnergy;
205 std::cout <<
"Created "<< nCalibrators <<
" calibrators.\n";
206 std::cout <<
"Address book size: \t\t"<< myAddressBook.size() <<
"\n";
207 std::cout <<
"Known space voxels size: \t"<< myKnownSpaceVoxels.size()
209 assert(myAddressBook.size() == myKnownSpaceVoxels.size());
210 makeInverseAddressBook();
219 c.reset(toClone.
clone());
230 const double energy)
const {
232 for (std::vector<SpaceVoxelPtr>::const_iterator cit =
235 if (s->contains(eta, phi, energy)) {
247 const std::map<DetectorElementPtr, double>&
result) {
259 std::vector<SpaceVoxelPtr> region;
277 TH2F hDist(name.c_str(), name.c_str(), 100, 0, 300, 50, 0.0, 2.5);
280 for (std::vector<SpaceVoxelPtr>::iterator
i = region.begin();
i
281 != region.end(); ++
i) {
284 if (s->maxEnergy() > maxE)
285 maxE = s->maxEnergy();
286 if (s->minEnergy() < minE)
287 minE = s->minEnergy();
291 std::vector<ParticleDepositPtr> particles = c->getParticles();
292 for (std::vector<ParticleDepositPtr>::iterator it =
293 particles.begin(); it != particles.end(); ++it) {
295 hDist.Fill((*it)->getTruthEnergy(), coeff);
297 hDist.Fill((*it)->getRecEnergy(), coeff);
306 nameProfile.append(
"_pfx");
309 TH1D* slices = (TH1D*) gDirectory->Get(name.c_str());
311 TProfile* profile = (TProfile*) gDirectory->Get(nameProfile.c_str());
324 return func.Eval(energy);
327 return func.Eval(energy);
336 double midEnergy = (s->minEnergy() + s->maxEnergy())/2.0;
338 double diffEnergy = energy - midEnergy;
341 double interpolatedCoeff = thisCoeff;
342 double adjacentCoeff = 0.0;
343 double adjacentEnergy = 0.0;
344 if (diffEnergy > 0) {
350 adjacentEnergy = (adjS->minEnergy() + adjS->maxEnergy()) / 2.0;
358 adjacentEnergy = (adjS->minEnergy() + adjS->maxEnergy()) / 2.0;
361 if (adjacentCoeff != 0) {
362 interpolatedCoeff = thisCoeff + diffEnergy* (adjacentCoeff - thisCoeff)
363 / (adjacentEnergy - midEnergy);
365 return interpolatedCoeff;
369 stream <<
"Calibration results: \n";
372 stream <<
"WARNING! Haven't sorted space voxels properly!\n";
377 stream << *s <<
"\n";
380 for (std::map<DetectorElementPtr, double>::iterator
b = result.begin();
b
381 != result.end(); ++
b) {
383 stream << *d <<
": ";
384 double ans = (*b).second;
385 stream << ans <<
", ";
395 for (std::map<SpaceVoxelPtr, CalibratorPtr>::iterator it =
U second(std::pair< T, U > const &p)
void clear(CLHEP::HepGenMatrix &m)
Helper function: Reset all elements of a matrix to 0.