File indexing completed on 2025-01-18 09:14:55
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015 #include <DD4hep/Factories.h>
0016 #include <DD4hep/Printout.h>
0017 #include <DDDigi/noise/FalphaNoise.h>
0018
0019
0020 #include <random>
0021 #include <iostream>
0022
0023 #include <TH1.h>
0024 #include <TCanvas.h>
0025 #include <TPad.h>
0026
0027 using namespace dd4hep;
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037 static long test_FalphaNoise(Detector& , int argc, char** argv) {
0038 using namespace dd4hep::detail;
0039 bool draw = false;
0040 size_t poles = 5;
0041 size_t shots = 10;
0042 double alpha = 0.5, variance = 1.0;
0043 for(int i = 0; i < argc && argv[i]; ++i) {
0044 if ( 0 == ::strncmp("-shots",argv[i],3) )
0045 shots = ::atol(argv[++i]);
0046 else if ( 0 == ::strncmp("-alpha",argv[i],3) )
0047 alpha = ::atof(argv[++i]);
0048 else if ( 0 == ::strncmp("-variance",argv[i],3) )
0049 variance = ::atof(argv[++i]);
0050 else if ( 0 == ::strncmp("-poles",argv[i],3) )
0051 poles = ::atol(argv[++i]);
0052 else if ( 0 == ::strncmp("-draw",argv[i],3) )
0053 draw = true;
0054 else {
0055 std::cout <<
0056 "Usage: -plugin DD4hep_FalphaNoise -arg [-arg] \n"
0057 " -shots <value> Number of samples to be generated [default: 10] \n"
0058 " -alpha <value> Parameter for the 1/f**alpha distribution \n"
0059 " default: 0.5 \n"
0060 " -variance <value> Distribution variance [default: 1] \n"
0061 " -poles <value> Number of IRR poles to fill the distribution \n"
0062 " -draw Fill and draw hoistogram with distribution \n"
0063 "\tArguments given: " << arguments(argc,argv) << std::endl << std::flush;
0064 ::exit(EINVAL);
0065 }
0066 }
0067
0068 std::default_random_engine generator;
0069 FalphaNoise noise(poles, alpha, variance);
0070 FalphaNoise noise0(poles, 0, variance);
0071 FalphaNoise noise1(poles, 1, variance);
0072 FalphaNoise noise2(poles, 1.98, variance);
0073 std::stringstream cpara;
0074 cpara << " distribution alpha=" << alpha << " sigma=" << variance;
0075 TH1D* hist11 = new TH1D("D11", ("1/f**alpha"+cpara.str()).c_str(), 50, -5e0*variance, 5e0*variance);
0076 TH1D* hist12 = new TH1D("D12", ("1/f**alpha"+cpara.str()).c_str(), 50, -5e0*variance, 5e0*variance);
0077 TH1D* hist13 = new TH1D("D13", ("1/f**alpha"+cpara.str()).c_str(), 50, -5e0*variance, 5e0*variance);
0078 TH1D* hist14 = new TH1D("D14", ("1/f**alpha"+cpara.str()).c_str(), 50, -5e0*variance, 5e0*variance);
0079 TH1D* hist21 = new TH1D("D21", ("log10(1/f**alpha)"+cpara.str()).c_str(), 300, 1e0, 1e2*variance);
0080 TH1D* hist22 = new TH1D("D22", "log10(1/f**alpha) alpha=0", 300, 1e0, 1e2*variance);
0081 TH1D* hist23 = new TH1D("D23", "log10(1/f**alpha) alpha=1", 300, 1e0, 1e2*variance);
0082 TH1D* hist24 = new TH1D("D24", "log10(1/f**alpha) alpha=2", 300, 1e0, 1e2*variance);
0083
0084 printout(INFO, "FalphaNoise", "Executing %ld shots .... variance=%.3f alpha=%.3f",
0085 shots, noise.variance(), noise.alpha());
0086 for(size_t i=0; i < shots; ++i) {
0087 double val = noise(generator);
0088 hist11->Fill(val, 1.);
0089 hist21->Fill(std::exp(std::abs(val)), 1.);
0090
0091 val = noise0(generator);
0092 hist12->Fill(val, 1.);
0093 hist22->Fill(std::exp(std::abs(val)), 1.);
0094
0095 val = noise1(generator);
0096 hist13->Fill(val, 1.);
0097 hist23->Fill(std::exp(std::abs(val)), 1.);
0098
0099 val = noise2(generator);
0100 hist14->Fill(val, 1.);
0101 hist24->Fill(std::exp(std::abs(val)), 1.);
0102 }
0103
0104 printout(INFO, "FalphaNoise", "Distribution Mean %10.5f", hist11->GetMean());
0105 printout(INFO, "FalphaNoise", "Distribution Mean Uncertainty %10.5f", hist11->GetMeanError());
0106 printout(INFO, "FalphaNoise", "Distribution RMS %10.5f", hist11->GetRMS());
0107 printout(INFO, "FalphaNoise", "Distribution RMS Uncertainty %10.5f", hist11->GetRMSError());
0108 hist11->GetXaxis()->SetTitle("Energy [arb.units]");
0109 hist11->GetYaxis()->SetTitle("Counts");
0110 hist21->GetXaxis()->SetTitle("exp(Energy) [arb.units]");
0111 hist21->GetYaxis()->SetTitle("Counts");
0112
0113 hist21->Scale(1.0 / hist21->GetBinContent(1));
0114 hist22->Scale(1.0 / hist22->GetBinContent(1));
0115 hist23->Scale(1.0 / hist23->GetBinContent(1));
0116 hist24->Scale(1.0 / hist24->GetBinContent(1));
0117 if ( draw ) {
0118 TCanvas *c = new TCanvas("c1","multigraph",700,500);
0119 c->Divide(1,2);
0120 c->cd(1);
0121
0122 gPad->SetGrid();
0123 hist12->SetLineColor(kBlue);
0124 hist12->Draw("L");
0125
0126 hist11->SetLineWidth(1);
0127 hist11->SetLineColor(kRed);
0128 hist11->Draw("SAME");
0129
0130 hist11->SetLineWidth(1);
0131 hist11->SetLineColor(kRed);
0132 hist11->Draw("LSAME");
0133
0134 hist13->SetLineColor(kMagenta);
0135 hist13->Draw("LSAME");
0136
0137 hist14->SetLineColor(kGreen);
0138 hist14->Draw("LSAME");
0139
0140 c->cd(2);
0141 gPad->SetGrid();
0142 gPad->SetLogx();
0143 gPad->SetLogy();
0144 hist21->Draw("");
0145 hist21->SetLineWidth(1);
0146 hist21->SetLineColor(kRed);
0147 hist21->Draw("LPSAME");
0148
0149 hist24->SetLineWidth(1);
0150 hist24->SetLineColor(kGreen);
0151 hist24->Draw("LPSAME");
0152
0153 hist22->SetLineWidth(1);
0154 hist22->SetLineColor(kBlue);
0155 hist22->Draw("LPSAME");
0156 hist23->SetLineWidth(1);
0157 hist23->SetLineColor(kMagenta);
0158 hist23->Draw("LPSAME");
0159 gPad->Update();
0160 gPad->Modified();
0161 }
0162 return 1;
0163 }
0164 DECLARE_APPLY(DD4hep_FalphaNoise,test_FalphaNoise)
