File indexing completed on 2023-10-25 10:03:35
0001 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0002 #include "SimCalorimetry/EcalSimAlgos/interface/EcalCorrelatedNoiseMatrix.h"
0003 #include "CalibFormats/CaloObjects/interface/CaloSamples.h"
0004 #include "SimGeneral/NoiseGenerators/interface/CorrelatedNoisifier.h"
0005 #include "DataFormats/Math/interface/Error.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007
0008 #include "CLHEP/Random/JamesRandom.h"
0009
0010 #include "TROOT.h"
0011 #include "TStyle.h"
0012 #include "TH1F.h"
0013 #include "TCanvas.h"
0014
0015 #include <iostream>
0016 #include <iomanip>
0017 #include <fstream>
0018
0019 int main() {
0020 edm::MessageDrop::instance()->debugEnabled = false;
0021
0022 const long seed = 12345;
0023 CLHEP::HepJamesRandom engine(seed);
0024
0025 const unsigned int readoutFrameSize = CaloSamples::MAXSAMPLES;
0026
0027 EcalCorrMatrix thisMatrix;
0028
0029 thisMatrix(0, 0) = 1.;
0030 thisMatrix(0, 1) = 0.67;
0031 thisMatrix(0, 2) = 0.53;
0032 thisMatrix(0, 3) = 0.44;
0033 thisMatrix(0, 4) = 0.39;
0034 thisMatrix(0, 5) = 0.36;
0035 thisMatrix(0, 6) = 0.38;
0036 thisMatrix(0, 7) = 0.35;
0037 thisMatrix(0, 8) = 0.36;
0038 thisMatrix(0, 9) = 0.32;
0039
0040 thisMatrix(1, 0) = 0.67;
0041 thisMatrix(1, 1) = 1.;
0042 thisMatrix(1, 2) = 0.67;
0043 thisMatrix(1, 3) = 0.53;
0044 thisMatrix(1, 4) = 0.44;
0045 thisMatrix(1, 5) = 0.39;
0046 thisMatrix(1, 6) = 0.36;
0047 thisMatrix(1, 7) = 0.38;
0048 thisMatrix(1, 8) = 0.35;
0049 thisMatrix(1, 9) = 0.36;
0050
0051 thisMatrix(2, 0) = 0.53;
0052 thisMatrix(2, 1) = 0.67;
0053 thisMatrix(2, 2) = 1.;
0054 thisMatrix(2, 3) = 0.67;
0055 thisMatrix(2, 4) = 0.53;
0056 thisMatrix(2, 5) = 0.44;
0057 thisMatrix(2, 6) = 0.39;
0058 thisMatrix(2, 7) = 0.36;
0059 thisMatrix(2, 8) = 0.38;
0060 thisMatrix(2, 9) = 0.35;
0061
0062 thisMatrix(3, 0) = 0.44;
0063 thisMatrix(3, 1) = 0.53;
0064 thisMatrix(3, 2) = 0.67;
0065 thisMatrix(3, 3) = 1.;
0066 thisMatrix(3, 4) = 0.67;
0067 thisMatrix(3, 5) = 0.53;
0068 thisMatrix(3, 6) = 0.44;
0069 thisMatrix(3, 7) = 0.39;
0070 thisMatrix(3, 8) = 0.36;
0071 thisMatrix(3, 9) = 0.38;
0072
0073 thisMatrix(4, 0) = 0.39;
0074 thisMatrix(4, 1) = 0.44;
0075 thisMatrix(4, 2) = 0.53;
0076 thisMatrix(4, 3) = 0.67;
0077 thisMatrix(4, 4) = 1.;
0078 thisMatrix(4, 5) = 0.67;
0079 thisMatrix(4, 6) = 0.53;
0080 thisMatrix(4, 7) = 0.44;
0081 thisMatrix(4, 8) = 0.39;
0082 thisMatrix(4, 9) = 0.36;
0083
0084 thisMatrix(5, 0) = 0.36;
0085 thisMatrix(5, 1) = 0.39;
0086 thisMatrix(5, 2) = 0.44;
0087 thisMatrix(5, 3) = 0.53;
0088 thisMatrix(5, 4) = 0.67;
0089 thisMatrix(5, 5) = 1.;
0090 thisMatrix(5, 6) = 0.67;
0091 thisMatrix(5, 7) = 0.53;
0092 thisMatrix(5, 8) = 0.44;
0093 thisMatrix(5, 9) = 0.39;
0094
0095 thisMatrix(6, 0) = 0.38;
0096 thisMatrix(6, 1) = 0.36;
0097 thisMatrix(6, 2) = 0.39;
0098 thisMatrix(6, 3) = 0.44;
0099 thisMatrix(6, 4) = 0.53;
0100 thisMatrix(6, 5) = 0.67;
0101 thisMatrix(6, 6) = 1.;
0102 thisMatrix(6, 7) = 0.67;
0103 thisMatrix(6, 8) = 0.53;
0104 thisMatrix(6, 9) = 0.44;
0105
0106 thisMatrix(7, 0) = 0.35;
0107 thisMatrix(7, 1) = 0.38;
0108 thisMatrix(7, 2) = 0.36;
0109 thisMatrix(7, 3) = 0.39;
0110 thisMatrix(7, 4) = 0.44;
0111 thisMatrix(7, 5) = 0.53;
0112 thisMatrix(7, 6) = 0.67;
0113 thisMatrix(7, 7) = 1.;
0114 thisMatrix(7, 8) = 0.67;
0115 thisMatrix(7, 9) = 0.53;
0116
0117 thisMatrix(8, 0) = 0.36;
0118 thisMatrix(8, 1) = 0.35;
0119 thisMatrix(8, 2) = 0.38;
0120 thisMatrix(8, 3) = 0.36;
0121 thisMatrix(8, 4) = 0.39;
0122 thisMatrix(8, 5) = 0.44;
0123 thisMatrix(8, 6) = 0.53;
0124 thisMatrix(8, 7) = 0.67;
0125 thisMatrix(8, 8) = 1.;
0126 thisMatrix(8, 9) = 0.67;
0127
0128 thisMatrix(9, 0) = 0.32;
0129 thisMatrix(9, 1) = 0.36;
0130 thisMatrix(9, 2) = 0.35;
0131 thisMatrix(9, 3) = 0.38;
0132 thisMatrix(9, 4) = 0.36;
0133 thisMatrix(9, 5) = 0.39;
0134 thisMatrix(9, 6) = 0.44;
0135 thisMatrix(9, 7) = 0.53;
0136 thisMatrix(9, 8) = 0.67;
0137 thisMatrix(9, 9) = 1.;
0138
0139 CorrelatedNoisifier<EcalCorrMatrix> theCorrNoise(thisMatrix);
0140
0141 EcalCorrMatrix thisTrivialMatrix;
0142 for (unsigned int i = 0; i < readoutFrameSize; i++) {
0143 thisTrivialMatrix(i, i) = 1.;
0144 for (unsigned int j = i + 1; j < readoutFrameSize; j++) {
0145 thisTrivialMatrix(i, j) = 0.;
0146 thisTrivialMatrix(j, i) = 0.;
0147 }
0148 }
0149 CorrelatedNoisifier<EcalCorrMatrix> theUncorrNoise(thisTrivialMatrix);
0150
0151 std::cout << "Using the correlation matrix: " << thisMatrix << std::endl;
0152
0153 std::cout << "\n And the unit matrix: " << thisTrivialMatrix << std::endl;
0154
0155 EBDetId detId(1, 1);
0156
0157 TH1F* uncorr = new TH1F("uncorr", "first 3 samples, uncorrelated distribution", 200, -10., 10.);
0158 TH1F* corr = new TH1F("corr", "first 3 samples, correlated distribution", 200, -10., 10.);
0159
0160 TH1F* uncorrPed = new TH1F("uncorrPed", "first 3 samples, uncorrelated average", 200, -10., 10.);
0161 TH1F* corrPed = new TH1F("corrPed", "first 3 samples, correlated average", 200, -10., 10.);
0162
0163 TH1F* frame[(int)readoutFrameSize];
0164 Char_t histo[200];
0165 for (Int_t i = 0; i < (int)readoutFrameSize; ++i) {
0166 sprintf(histo, "frame %02d", i);
0167 frame[i] = new TH1F(histo, histo, 200, -10., 10.);
0168 }
0169
0170 for (int i = 0; i < 100000; ++i) {
0171 CaloSamples noiseframe(detId, readoutFrameSize);
0172 theCorrNoise.noisify(noiseframe, &engine);
0173 CaloSamples flatframe(detId, readoutFrameSize);
0174 theUncorrNoise.noisify(flatframe, &engine);
0175 for (int j = 0; j < 3; ++j) {
0176 uncorr->Fill(flatframe[j]);
0177 corr->Fill(noiseframe[j]);
0178 }
0179 for (int j = 0; j < (int)readoutFrameSize; ++j) {
0180 frame[j]->Fill(noiseframe[j]);
0181 }
0182 float thisUncorrPed = (flatframe[0] + flatframe[1] + flatframe[2]) / 3.;
0183 float thisCorrPed = (noiseframe[0] + noiseframe[1] + noiseframe[2]) / 3.;
0184 uncorrPed->Fill(thisUncorrPed);
0185 corrPed->Fill(thisCorrPed);
0186 }
0187
0188 uncorr->Fit("gaus", "VE");
0189 corr->Fit("gaus", "VE");
0190
0191 uncorrPed->Fit("gaus", "VE");
0192 corrPed->Fit("gaus", "VE");
0193
0194 const int csize = 500;
0195 TCanvas* showNoise = new TCanvas("showNoise", "showNoise", 2 * csize, 2 * csize);
0196
0197 showNoise->Divide(2, 2);
0198
0199 gStyle->SetOptFit(1111);
0200
0201 showNoise->cd(1);
0202 uncorr->Draw();
0203 showNoise->cd(2);
0204 corr->Draw();
0205 showNoise->cd(3);
0206 uncorrPed->Draw();
0207 showNoise->cd(4);
0208 corrPed->Draw();
0209 showNoise->SaveAs("EcalNoise.jpg");
0210
0211 TCanvas* showNoiseFrame = new TCanvas("showNoiseFrame", "showNoiseFrame", 2 * csize, 2 * csize);
0212 showNoiseFrame->Divide(2, 5);
0213 for (int i = 0; i < (int)readoutFrameSize; ++i) {
0214 showNoiseFrame->cd(i + 1);
0215 frame[i]->Fit("gaus", "Q");
0216 frame[i]->Draw();
0217 }
0218 showNoiseFrame->SaveAs("EcalNoiseFrame.jpg");
0219
0220 delete uncorr;
0221 delete corr;
0222 delete uncorrPed;
0223 delete corrPed;
0224 delete showNoise;
0225 delete showNoiseFrame;
0226
0227 return 0;
0228 }