EcalClusters/test/photonContainmentAnalysis.C

0001 // fitting function
0002 double crystalball(double *x, double *par) {
0003
0004   // par[0]:  mean
0005   // par[1]:  sigma
0006   // par[2]:  alpha, crossover point
0007   // par[3]:  n, length of tail
0008   // par[4]:  N, normalization
0009
0010   double cb = 0.0;
0011   double exponent = 0.0;
0012
0013 /*   std::cout << " x = " << x[0] << " Par = "  */
0014 /*             << par[0] << " "  */
0015 /*             << par[1] << " "  */
0016 /*             << par[2] << " "  */
0017 /*             << par[3] << " "  */
0018 /*             << par[4] << std::endl; */
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   //  std::cout << "CB = " << std::endl;
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   // fitting the distributions to extract the parameter
0049   TF1 *gausa;
0050   TF1 *cb_p;
0051   for(int myH=0; myH<2; myH++) {
0052
0053     // histos parameters
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     // gaussian fit to initialize
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     // crystalball limits
0071     double myXmin = gausMean - 7.*gausSigma;
0072     double myXmax = gausMean + 5.*gausSigma;
0073
0074     // crystalball fit
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 }