EcalSimAlgos/test/testCorrNoise.cpp

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 }