File indexing completed on 2024-04-06 12:32:02
0001
0002 double crystalball(double *x, double *par) {
0003
0004
0005
0006
0007
0008
0009
0010 double cb = 0.0;
0011 double exponent = 0.0;
0012
0013
0014
0015
0016
0017
0018
0019
0020 if (x[0] > par[0] - par[2]*par[1]) {
0021 exponent = (x[0] - par[0])/par[1];
0022 cb = exp(-exponent*exponent/2.);
0023 }
0024 else {
0025 double nenner = pow(par[3]/par[2], par[3])*exp(-par[2]*par[2]/2.);
0026 double zaehler = (par[0] - x[0])/par[1] + par[3]/par[2] - par[2];
0027 zaehler = pow(zaehler, par[3]);
0028 cb = nenner/zaehler;
0029 }
0030 if (par[4] > 0.) { cb *= par[4]; }
0031
0032
0033
0034 return cb;
0035 }
0036
0037 void photonContainmentAnalysis() {
0038
0039 gStyle->SetOptFit(1111);
0040
0041 TFile *infile = new TFile("contCorrAnalyzer.root");
0042 infile->ls();
0043
0044 TH1F *theHistos[2];
0045 theHistos[0] = (TH1F*)infile->Get("contCorrAnalyzer/EB_e25EtrueReference");
0046 theHistos[1] = (TH1F*)infile->Get("contCorrAnalyzer/EE_e25EtrueReference");
0047
0048
0049 TF1 *gausa;
0050 TF1 *cb_p;
0051 for(int myH=0; myH<2; myH++) {
0052
0053
0054 int peakBin = theHistos[myH]->GetMaximumBin();
0055 double h_norm = theHistos[myH]->GetMaximum();
0056 double h_rms = theHistos[myH]->GetRMS();
0057 double h_mean = theHistos[myH]->GetMean();
0058 double h_peak = theHistos[myH]->GetBinCenter(peakBin);
0059
0060
0061 gausa = new TF1 ("gausa","[0]*exp(-1*(x-[1])*(x-[1])/2/[2]/[2])",h_peak-10*h_rms,h_peak+10*h_rms);
0062 gausa ->SetParameters(h_norm,h_peak,h_rms);
0063 theHistos[myH]->Fit("gausa","","",h_peak-3*h_rms,h_peak+3*h_rms);
0064 double gausNorm = gausa->GetParameter(0);
0065 double gausMean = gausa->GetParameter(1);
0066 double gausSigma = fabs(gausa->GetParameter(2));
0067 double gausChi2 = gausa->GetChisquare()/gausa->GetNDF();
0068 if (gausChi2>100){ gausMean = h_peak; gausSigma = h_rms; }
0069
0070
0071 double myXmin = gausMean - 7.*gausSigma;
0072 double myXmax = gausMean + 5.*gausSigma;
0073
0074
0075 cb_p = new TF1 ("cb_p",crystalball,myXmin,myXmax, 5) ;
0076 cb_p->SetParNames ("Mean","Sigma","alpha","n","Norm","Constant");
0077 cb_p->SetParameter(0, gausMean);
0078 cb_p->SetParameter(1, gausSigma);
0079 cb_p->SetParameter(2, 1.);
0080 cb_p->SetParLimits(2, 0.1, 5.);
0081 cb_p->FixParameter(3, 5.);
0082 cb_p->SetParameter(4, gausNorm);
0083 theHistos[myH]->Fit("cb_p","lR","",myXmin,myXmax);
0084 theHistos[myH]->GetXaxis()->SetRangeUser(0.95,1.05);
0085 double matrix_gmean = cb_p->GetParameter(0);
0086 double matrix_gmean_err = cb_p->GetParError(0);
0087
0088 c1->Update();
0089 if(myH == 0) {
0090 c1->SaveAs("e25EtrueReference_EB.eps");
0091 std::cout << "E25 barrel containment: " << matrix_gmean << " +/- " << matrix_gmean_err << std::endl;
0092 }
0093 if(myH == 1) {
0094 c1->SaveAs("e25EtrueReference_EE.eps");
0095 std::cout << "E25 endcap containment: " << matrix_gmean << " +/- " << matrix_gmean_err << std::endl;
0096 }
0097 }
0098 }