Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 10:43:02 

0001 #include "TFile.h"
0002 #include "TH1F.h"
0003 #include "TH2F.h"
0004 #include "TTree.h"
0005 #include "TProfile.h"
0006 #include "TPostScript.h"
0007 #include "TCanvas.h"
0008 #include "TF1.h"
0009 #include "TStyle.h"
0010 #include "TLatex.h"
0011 #include "TLegend.h"
0012 #include "TPaveStats.h"
0013 
0014 #include "TRandom.h"
0015 
0016 #include <string>
0017 
0018 #include <iostream>
0019 #include <fstream>
0020 #include <iomanip>
0021 #include <vector>
0022 #include "TMinuit.h"
0023 #include "TMath.h"
0024 
0025 const char* prehist = "hoCalibc/";
0026 
0027 using namespace std;
0028 static unsigned int mypow_2[32];
0029 
0030 int irunold = 1000;
0031 const int nmxbin = 40000;
0032 float yvalft[nmxbin];
0033 float xvalft[nmxbin];
0034 int nvalft = 0;
0035 
0036 float alowx = -1.0;
0037 float ahighx = 29.0;
0038 int auprange = 25;
0039 
0040 double fitchisq;
0041 int fitndof;
0042 double anormglb = -1;
0043 Double_t gausX(Double_t* x, Double_t* par) { return par[0] * (TMath::Gaus(x[0], par[1], par[2], kTRUE)); }
0044 
0045 // Double_t landX(Double_t* x, Double_t* par) {
0046 //   return par[0]*(TMath::Landau(x[0], par[1], par[2]));
0047 // }
0048 
0049 // Double_t completefit(Double_t* x, Double_t* par) {
0050 //   return gausX(x, par) + landX(x, &par[3]);
0051 // }
0052 
0053 Double_t langaufun(Double_t* x, Double_t* par) {
0054   //Fit parameters:
0055   //par[0]*par[1]=Width (scale) parameter of Landau density
0056   //par[1]=Most Probable (MP, location) parameter of Landau density
0057   //par[2]=Total area (integral -inf to inf, normalization constant)
0058   //par[3]=Width (sigma) of convoluted Gaussian function
0059   //
0060   //In the Landau distribution (represented by the CERNLIB approximation),
0061   //the maximum is located at x=-0.22278298 with the location parameter=0.
0062   //This shift is corrected within this function, so that the actual
0063   //maximum is identical to the MP parameter.
0064   // /*
0065   // Numeric constants
0066   Double_t invsq2pi = 0.3989422804014;  // (2 pi)^(-1/2)
0067   Double_t mpshift = -0.22278298;       // Landau maximum location
0068 
0069   // Control constants
0070   Double_t np = 100.0;  // number of convolution steps
0071   Double_t sc = 5.0;    // convolution extends to +-sc Gaussian sigmas
0072 
0073   // Variables
0074   Double_t xx;
0075   Double_t mpc;
0076   Double_t fland;
0077   Double_t sum = 0.0;
0078   Double_t xlow, xupp;
0079   Double_t step;
0080 
0081   // MP shift correction
0082   mpc = par[1] - mpshift * par[0] * par[1];
0083   double scale = 1;          // par[1];
0084   double scale2 = anormglb;  // for notmalisation this is one, otehrwise use the normalisation;
0085   if (scale2 < .1)
0086     scale2 = 0.1;
0087   //  double scale=par[1];
0088   // Range of convolution integral
0089   xlow = x[0] - sc * scale * par[3];
0090   xupp = x[0] + sc * scale * par[3];
0091 
0092   step = (xupp - xlow) / np;
0093 
0094   // Convolution integral of Landau and Gaussian by sum
0095   for (double ij = 1.0; ij <= np / 2; ij++) {
0096     xx = xlow + (ij - .5) * step;
0097     fland = TMath::Landau(xx, mpc, par[0] * par[1], kTRUE);  // / par[0];
0098     if (xx > par[1]) {
0099       fland *= exp(-(xx - par[1]) / par[4]);
0100     }
0101     sum += fland * TMath::Gaus(x[0], xx, scale * par[3]);
0102     xx = xupp - (ij - .5) * step;
0103     fland = TMath::Landau(xx, mpc, par[0] * par[1], kTRUE);  // / par[0];
0104     if (xx > par[1]) {
0105       fland *= exp(-(xx - par[1]) / par[4]);
0106     }
0107     sum += fland * TMath::Gaus(x[0], xx, scale * par[3]);
0108   }
0109   return (par[2] * step * sum * invsq2pi / (scale2 * par[3]));
0110 }
0111 
0112 Double_t totalfunc(Double_t* x, Double_t* par) {
0113   return gausX(x, par) + langaufun(x, &par[3]);  // /max(.001,anormglb);
0114 }
0115 
0116 const int netamx = 30;
0117 const int nphimx = 72;
0118 
0119 int getHOieta(int ij) { return (ij < netamx / 2) ? -netamx / 2 + ij : -netamx / 2 + ij + 1; }
0120 int invert_HOieta(int ieta) { return (ieta < 0) ? netamx / 2 + ieta : netamx / 2 + ieta - 1; }
0121 
0122 int ietafit;
0123 int iphifit;
0124 
0125 void fcnsg(Int_t& npar, Double_t* gin, Double_t& f, Double_t* par, Int_t flag) {
0126   double xval[2];
0127 
0128   anormglb = 1;
0129   double tmpnorm = 0;
0130   for (int ij = 0; ij < 600; ij++) {
0131     xval[0] = -5 + (ij + 0.5) * .1;
0132     tmpnorm += 0.1 * langaufun(xval, &par[3]) / par[5];
0133   }
0134   anormglb = tmpnorm;
0135 
0136   double fval = -(par[0] + par[5]);
0137   for (int ij = 0; ij < nvalft; ij++) {
0138     if (yvalft[ij] < 1 || xvalft[ij] > auprange)
0139       continue;
0140     xval[0] = xvalft[ij];
0141     fval += yvalft[ij] * log(max(1.e-30, totalfunc(xval, par)));
0142   }
0143   f = -fval;
0144 }
0145 
0146 int main() {
0147   int icol = 0;
0148   int ntotal = 0;
0149   int irebin = 1;
0150   int ifile = 0;
0151   for (int ij = 0; ij < 32; ij++) {
0152     mypow_2[ij] = (int)pow(float(2), ij);
0153   }
0154 
0155   int max_nEvents = 300000000;
0156   int ityp = 1;
0157   bool m_cosmic = false;
0158   bool m_histFill = true;
0159   bool m_treeFill = false;
0160   bool m_zeroField = false;
0161   float pival = acos(-1.);
0162 
0163   int nsize;
0164   char outfile[100];
0165   char outfilx[100];
0166   char nametag[100];
0167   char nametag2[100];
0168   char infile[200];
0169   char datafile[100];
0170   char histname[100];
0171   char name[100];
0172   char title[100];
0173 
0174   cout << "Give the value of rebinning: a nonzero positive integer" << endl;
0175   cin >> irebin;
0176   if (irebin < 1) {
0177     irebin = 1;
0178   }
0179   cout << "Give the upper range of signal for fit [5 - 29]" << endl;
0180   cin >> auprange;
0181 
0182   cout << "Give the histname, e.g., 2017b_v2 2018a_v3a" << endl;
0183   cin >> nametag2;
0184 
0185   int len = strlen(nametag2);
0186   nametag2[len] = '\0';
0187 
0188   sprintf(outfilx, "hist_hoprompt_%s.root", nametag2);
0189   len = strlen(outfilx);
0190   outfilx[len] = '\0';
0191 
0192   //  TFile* fileIn = new TFile("histalla_hcalcalho_2016bcde_v2a.root", "read");
0193 
0194   TFile* fileIn = new TFile(outfilx, "read");
0195   //  const char* nametag="2016e_ar1";
0196   //  sprintf(nametag, "2016b_ar%i", irebin);
0197 
0198   unsigned ievt, hoflag;
0199   int irun, ilumi, nprim, isect, isect2, ndof, nmuon;
0200 
0201   float inslumi, trkdr, trkdz, trkvx, trkvy, trkvz, trkmm, trkth, trkph, chisq, therr, pherr, hodx, hody, hoang, htime,
0202       hosig[9], hocorsig[18], hocro, hbhesig[9], caloen[3];
0203   float momatho, tkpt03, ecal03, hcal03;
0204   float tmphoang;
0205 
0206   TTree* Tin;
0207   if (m_treeFill) {
0208     Tin = (TTree*)fileIn->Get("T1");
0209 
0210     Tin->SetBranchAddress("irun", &irun);
0211     Tin->SetBranchAddress("ievt", &ievt);
0212 
0213     Tin->SetBranchAddress("isect", &isect);
0214     Tin->SetBranchAddress("isect2", &isect2);
0215     Tin->SetBranchAddress("ndof", &ndof);
0216     Tin->SetBranchAddress("nmuon", &nmuon);
0217 
0218     Tin->SetBranchAddress("ilumi", &ilumi);
0219     if (!m_cosmic) {
0220       Tin->SetBranchAddress("inslumi", &inslumi);
0221       Tin->SetBranchAddress("nprim", &nprim);
0222       Tin->SetBranchAddress("tkpt03", &tkpt03);
0223       Tin->SetBranchAddress("ecal03", &ecal03);
0224       Tin->SetBranchAddress("hcal03", &hcal03);
0225     }
0226 
0227     Tin->SetBranchAddress("trkdr", &trkdr);
0228     Tin->SetBranchAddress("trkdz", &trkdz);
0229 
0230     Tin->SetBranchAddress("trkvx", &trkvx);
0231     Tin->SetBranchAddress("trkvy", &trkvy);
0232     Tin->SetBranchAddress("trkvz", &trkvz);
0233     Tin->SetBranchAddress("trkmm", &trkmm);
0234     Tin->SetBranchAddress("trkth", &trkth);
0235     Tin->SetBranchAddress("trkph", &trkph);
0236 
0237     Tin->SetBranchAddress("chisq", &chisq);
0238     Tin->SetBranchAddress("therr", &therr);
0239     Tin->SetBranchAddress("pherr", &pherr);
0240     Tin->SetBranchAddress("hodx", &hodx);
0241     Tin->SetBranchAddress("hody", &hody);
0242     Tin->SetBranchAddress("hoang", &hoang);
0243 
0244     Tin->SetBranchAddress("momatho", &momatho);
0245     Tin->SetBranchAddress("hoflag", &hoflag);
0246     Tin->SetBranchAddress("htime", &htime);
0247     Tin->SetBranchAddress("hosig", hosig);
0248     Tin->SetBranchAddress("hocro", &hocro);
0249     Tin->SetBranchAddress("hocorsig", hocorsig);
0250     Tin->SetBranchAddress("caloen", caloen);
0251   }
0252 
0253   sprintf(nametag, "%s_ar%i_float_par8_rng%i", nametag2, irebin, auprange);
0254 
0255   TLatex latex;
0256   latex.SetNDC();
0257   latex.SetTextSize(0.06);
0258   latex.SetTextFont(22);
0259   latex.SetTextAlign(11);  // 11 left; // 21 centre, // (31); // align right, 22, 23, shift bottom
0260 
0261   //Related fitting, not for storing
0262   const int nsample = 16;  //# of signal plots in the .ps page
0263   TF1* pedfun[nsample] = {0};
0264   TF1* sigfun[nsample] = {0};
0265   TF1* signalx[nsample] = {0};
0266 
0267   TH1F* signall[nsample] = {0};
0268 
0269   TH1F* signallunb[nsample] = {0};
0270 
0271   const int nbgpr = 3;
0272   const int nsgpr = 8;
0273   double fitprm[nsgpr][netamx];
0274   double xmn = -1.0;
0275   double xmx = 29.0;
0276 
0277   sprintf(outfilx, "fit_%s", nametag);
0278   len = strlen(outfilx);
0279   outfilx[len] = '\0';
0280 
0281   sprintf(outfile, "%s.txt", outfilx);
0282   ofstream file_out(outfile);
0283 
0284   sprintf(outfile, "%s.root", outfilx);
0285 
0286   TFile* fileOut = new TFile(outfile, "recreate");
0287 
0288   TTree* Tout;
0289 
0290   TH1F* sigrsg[netamx][nphimx];
0291 
0292   TH1F* fit_chi;
0293   TH1F* sig_evt;
0294   TH1F* fit_sigevt;
0295   TH1F* fit_bkgevt;
0296   TH1F* sig_mean;
0297   TH1F* sig_diff;
0298   TH1F* sig_width;
0299   TH1F* sig_sigma;
0300   TH1F* sig_expo;
0301   TH1F* sig_meanerr;
0302   TH1F* sig_meanerrp;
0303   TH1F* sig_signf;
0304 
0305   TH1F* sig_statmean;
0306   TH1F* sig_rms;
0307 
0308   TH2F* ped2d_evt;
0309   TH2F* ped2d_mean;
0310   TH2F* ped2d_width;
0311   TH2F* fit2d_chi;
0312   TH2F* sig2d_evt;
0313   TH2F* fit2d_sigevt;
0314   TH2F* fit2d_bkgevt;
0315   TH2F* sig2d_mean;
0316   TH2F* sig2d_diff;
0317   TH2F* sig2d_width;
0318   TH2F* sig2d_sigma;
0319   TH2F* sig2d_expo;
0320   TH2F* sig2d_meanerr;
0321   TH2F* sig2d_meanerrp;
0322   TH2F* sig2d_signf;
0323   TH2F* sig2d_rms;
0324   TH2F* sig2d_statmean;
0325 
0326   fit_chi = new TH1F("fit_chi", "fit_chi", netamx * nphimx, -0.5, netamx * nphimx - 0.5);
0327   sig_evt = new TH1F("sig_evt", "sig_evt", netamx * nphimx, -0.5, netamx * nphimx - 0.5);
0328   fit_sigevt = new TH1F("fit_sigevt", "fit_sigevt", netamx * nphimx, -0.5, netamx * nphimx - 0.5);
0329   fit_bkgevt = new TH1F("fit_bkgevt", "fit_bkgevt", netamx * nphimx, -0.5, netamx * nphimx - 0.5);
0330   sig_mean = new TH1F("sig_mean", "sig_mean", netamx * nphimx, -0.5, netamx * nphimx - 0.5);
0331   sig_diff = new TH1F("sig_diff", "sig_diff", netamx * nphimx, -0.5, netamx * nphimx - 0.5);
0332   sig_width = new TH1F("sig_width", "sig_width", netamx * nphimx, -0.5, netamx * nphimx - 0.5);
0333   sig_sigma = new TH1F("sig_sigma", "sig_sigma", netamx * nphimx, -0.5, netamx * nphimx - 0.5);
0334   sig_expo = new TH1F("sig_expo", "sig_expo", netamx * nphimx, -0.5, netamx * nphimx - 0.5);
0335 
0336   sig_meanerr = new TH1F("sig_meanerr", "sig_meanerr", netamx * nphimx, -0.5, netamx * nphimx - 0.5);
0337   sig_meanerrp = new TH1F("sig_meanerrp", "sig_meanerrp", netamx * nphimx, -0.5, netamx * nphimx - 0.5);
0338   sig_signf = new TH1F("sig_signf", "sig_signf", netamx * nphimx, -0.5, netamx * nphimx - 0.5);
0339 
0340   sig_statmean = new TH1F("sig_statmean", "sig_statmean", netamx * nphimx, -0.5, netamx * nphimx - 0.5);
0341   sig_rms = new TH1F("sig_rms", "sig_rms", netamx * nphimx, -0.5, netamx * nphimx - 0.5);
0342 
0343   fit2d_chi =
0344       new TH2F("fit2d_chi", "fit2d_chi", netamx + 1, -netamx / 2. - 0.5, netamx / 2. + 0.5, nphimx, 0.5, nphimx + 0.5);
0345   sig2d_evt =
0346       new TH2F("sig2d_evt", "sig2d_evt", netamx + 1, -netamx / 2. - 0.5, netamx / 2. + 0.5, nphimx, 0.5, nphimx + 0.5);
0347   fit2d_sigevt = new TH2F(
0348       "fit2d_sigevt", "fit2d_sigevt", netamx + 1, -netamx / 2. - 0.5, netamx / 2. + 0.5, nphimx, 0.5, nphimx + 0.5);
0349   fit2d_bkgevt = new TH2F(
0350       "fit2d_bkgevt", "fit2d_bkgevt", netamx + 1, -netamx / 2. - 0.5, netamx / 2. + 0.5, nphimx, 0.5, nphimx + 0.5);
0351   sig2d_mean = new TH2F(
0352       "sig2d_mean", "sig2d_mean", netamx + 1, -netamx / 2. - 0.5, netamx / 2. + 0.5, nphimx, 0.5, nphimx + 0.5);
0353   sig2d_diff = new TH2F(
0354       "sig2d_diff", "sig2d_diff", netamx + 1, -netamx / 2. - 0.5, netamx / 2. + 0.5, nphimx, 0.5, nphimx + 0.5);
0355   sig2d_width = new TH2F(
0356       "sig2d_width", "sig2d_width", netamx + 1, -netamx / 2. - 0.5, netamx / 2. + 0.5, nphimx, 0.5, nphimx + 0.5);
0357   sig2d_sigma = new TH2F(
0358       "sig2d_sigma", "sig2d_sigma", netamx + 1, -netamx / 2. - 0.5, netamx / 2. + 0.5, nphimx, 0.5, nphimx + 0.5);
0359   sig2d_expo = new TH2F(
0360       "sig2d_expo", "sig2d_expo", netamx + 1, -netamx / 2. - 0.5, netamx / 2. + 0.5, nphimx, 0.5, nphimx + 0.5);
0361 
0362   sig2d_statmean = new TH2F(
0363       "sig2d_statmean", "sig2d_statmean", netamx + 1, -netamx / 2. - 0.5, netamx / 2. + 0.5, nphimx, 0.5, nphimx + 0.5);
0364   sig2d_rms =
0365       new TH2F("sig2d_rms", "sig2d_rms", netamx + 1, -netamx / 2. - 0.5, netamx / 2. + 0.5, nphimx, 0.5, nphimx + 0.5);
0366 
0367   sig2d_meanerr = new TH2F(
0368       "sig2d_meanerr", "sig2d_meanerr", netamx + 1, -netamx / 2. - 0.5, netamx / 2. + 0.5, nphimx, 0.5, nphimx + 0.5);
0369   sig2d_meanerrp = new TH2F(
0370       "sig2d_meanerrp", "sig2d_meanerrp", netamx + 1, -netamx / 2. - 0.5, netamx / 2. + 0.5, nphimx, 0.5, nphimx + 0.5);
0371   sig2d_signf = new TH2F(
0372       "sig2d_signf", "sig2d_signf", netamx + 1, -netamx / 2. - 0.5, netamx / 2. + 0.5, nphimx, 0.5, nphimx + 0.5);
0373 
0374   if (m_histFill) {
0375     for (int jk = 0; jk < netamx; jk++) {
0376       for (int ij = 0; ij < nphimx; ij++) {
0377         //              sprintf(name, "%sho_indenergy_%i_%i", prehist, jk, ij);
0378         sprintf(name, "%ssig_eta%i_phi%i", prehist, jk, ij);
0379         sigrsg[jk][ij] = (TH1F*)fileIn->Get(name);
0380         cout << "jkij " << jk << " " << ij << " " << name << endl;
0381       }
0382     }
0383   } else if (m_treeFill) {
0384     for (int jk = 0; jk < netamx; jk++) {
0385       int ieta = getHOieta(jk);
0386       for (int ij = 0; ij < nphimx; ij++) {
0387         sprintf(name, "ho_indenergy_%i_%i", jk, ij);
0388         sprintf(title, "signal (i#eta=%i - i#phi=%i", ieta, ij + 1);
0389         sigrsg[jk][ij] = new TH1F(name, name, 120, -1.0, 14.0);  //1200, -1.0, 29.0);
0390       }
0391     }
0392 
0393     int nentries = Tin->GetEntries();
0394     for (int iev = 0; iev < nentries; iev++) {
0395       fileIn->cd();
0396       Tin->GetEntry(iev);
0397       fileOut->cd();
0398       int ieta = int((abs(isect) % 10000) / 100.) - 50;
0399       int iphi = abs(isect) % 100;
0400       int tmpxet = invert_HOieta(ieta);
0401 
0402       //event selection
0403       if (hosig[4] > -90) {
0404         sigrsg[tmpxet][iphi - 1]->Fill(hosig[4]);
0405       }
0406     }
0407   } else {
0408     cout << " You must read either histogramme or tree " << endl;
0409     return 1;
0410   }
0411 
0412   gStyle->SetTitleFontSize(0.075);
0413   gStyle->SetTitleBorderSize(1);
0414   gStyle->SetPadTopMargin(0.12);
0415   gStyle->SetPadBottomMargin(0.10);
0416   gStyle->SetPadLeftMargin(0.15);
0417   gStyle->SetPadRightMargin(0.03);
0418 
0419   gStyle->SetOptStat(1110);
0420   gStyle->SetLabelSize(0.095, "XYZ");
0421   gStyle->SetLabelOffset(-0.01, "XYZ");
0422   gStyle->SetHistLineColor(1);
0423   gStyle->SetHistLineWidth(2);
0424   gStyle->SetPadGridX(0);
0425   gStyle->SetPadGridY(0);
0426   gStyle->SetGridStyle(0);
0427   gStyle->SetOptLogy(1);
0428 
0429   int ips = 111;
0430   sprintf(outfile, "%s.ps", outfilx);
0431   TPostScript ps(outfile, ips);
0432   ps.Range(20, 28);
0433   gStyle->SetOptLogy(0);
0434   gStyle->SetPadLeftMargin(0.17);
0435 
0436   gStyle->SetPadGridX(3);
0437   gStyle->SetPadGridY(3);
0438   gStyle->SetGridStyle(2);
0439   gStyle->SetPadRightMargin(0.17);
0440   gStyle->SetPadLeftMargin(0.10);
0441 
0442   gStyle->SetTitleFontSize(0.045);
0443   gStyle->SetPadTopMargin(0.10);     //.12
0444   gStyle->SetPadBottomMargin(0.12);  //.14
0445   gStyle->SetPadLeftMargin(0.17);
0446   gStyle->SetPadRightMargin(0.03);
0447 
0448   gStyle->SetOptStat(0);  //GMA (1110);
0449 
0450   gStyle->SetOptFit(101);
0451   gStyle->SetCanvasBorderMode(0);
0452   gStyle->SetPadBorderMode(0);
0453   gStyle->SetStatBorderSize(1);
0454   gStyle->SetStatStyle(1001);
0455   gStyle->SetTitleColor(10);
0456   gStyle->SetTitleFontSize(0.09);
0457   gStyle->SetTitleOffset(-0.05);
0458   gStyle->SetTitleFillColor(10);
0459   gStyle->SetTitleBorderSize(1);
0460 
0461   gStyle->SetCanvasColor(10);
0462   gStyle->SetPadColor(10);
0463   gStyle->SetStatColor(10);
0464   //    gStyle->SetStatFontSize(.07);
0465   gStyle->SetStatX(0.99);
0466   gStyle->SetStatY(0.99);
0467   gStyle->SetStatW(0.44);
0468   gStyle->SetStatH(0.16);
0469   gStyle->SetTitleSize(0.065, "XYZ");
0470   gStyle->SetLabelSize(0.075, "XYZ");
0471   gStyle->SetLabelOffset(0.012, "XYZ");
0472   gStyle->SetPadGridX(0);  //(1)
0473   gStyle->SetPadGridY(0);  //(1)
0474   gStyle->SetGridStyle(3);
0475   gStyle->SetNdivisions(101, "XY");
0476   gStyle->SetOptLogy(1);  //0); //GMA 1
0477   int iiter = 0;
0478 
0479   ps.NewPage();
0480 
0481   int xsiz = 900;   //900;
0482   int ysiz = 1200;  //600;
0483   TCanvas* c0 = new TCanvas("c0", " Pedestal vs signal", xsiz, ysiz);
0484   c0->Divide(4, 4);
0485   fileOut->cd();
0486   for (int ij = 0; ij < nphimx; ij++) {
0487     int iphi = ij + 1;
0488     //      for (int jk=0; jk<netamx; jk++) {
0489     for (int jk = 7; jk < 8; jk++) {
0490       int ieta = getHOieta(jk);
0491       int izone = iiter % nsample;
0492 
0493       signall[izone] = (TH1F*)sigrsg[jk][ij]->Clone("hnew");
0494       if (irebin > 1) {
0495         signall[izone]->Rebin(irebin);
0496       }
0497 
0498       signallunb[izone] = (TH1F*)signall[izone]->Clone("hnew");
0499       signall[izone]->Rebin(10 / irebin);
0500 
0501       if (izone == 0) {  //iiter%8 ==0) {
0502         ps.NewPage();
0503       }
0504       c0->cd(izone + 1);
0505       if (signall[izone]->GetEntries() / 2. > 5) {
0506         double binwid = signall[izone]->GetBinWidth(1);
0507 
0508         Double_t parall[nsgpr];
0509         double parserr[nsgpr];
0510         double fitres[nsgpr];
0511         double pedht = 0;
0512 
0513         char temp[20];
0514         xmn = signall[izone]->GetXaxis()->GetXmin();
0515         xmx = signall[izone]->GetXaxis()->GetXmax();
0516         int nbn = signall[izone]->FindBin(0);
0517 
0518         cout << "bincenter ===================================== " << signall[izone]->GetBinCenter(nbn) << endl;
0519         pedht = (signall[izone]->GetBinContent(nbn - 1) + signall[izone]->GetBinContent(nbn) +
0520                  signall[izone]->GetBinContent(nbn + 1)) /
0521                 3.;
0522 
0523         parall[0] = max(1.0, 0.9 * pedht);               //Pedestal peak
0524         parall[1] = 0.00;                                //pedestal mean
0525         parall[2] = 0.03;                                //pedestal width
0526         parall[3] = 0.135;                               //Gaussian smearing of Landau function
0527         parall[4] = 0.7 * signallunb[izone]->GetMean();  //fitprm[4][jk]; //from 2015 cosmic data
0528         parall[5] = signallunb[izone]->GetEntries() / 2.;
0529         parall[6] = 0.2238;  // from 2015 cosmic data
0530         parall[7] = 5.0;
0531 
0532         nvalft = min(nmxbin, signallunb[izone]->GetNbinsX());
0533 
0534         for (int lm = 0; lm < nvalft; lm++) {
0535           xvalft[lm] = signallunb[izone]->GetBinCenter(lm + 1);
0536           yvalft[lm] = signallunb[izone]->GetBinContent(lm + 1);
0537         }
0538 
0539         TMinuit* gMinuit = new TMinuit(nsgpr);
0540         TString namex[nsgpr] = {"const", "mean", "sigma", "Width", "MP", "Area", "GSigma", "Exp"};
0541         double strt[nsgpr] = {parall[0], parall[1], parall[2], parall[3], parall[4], parall[5], parall[6], parall[7]};
0542         double alx = max(0.0, 0.1 * parall[0] - 1.1);
0543         double alowmn[nsgpr] = {
0544             alx, -0.1, -0.1, 0.06, 0.5 * strt[4] - 0.5, 0.1 * strt[5], 0.1 * strt[6], 0.5 * parall[7]};
0545         double ahighmn[nsgpr] = {
0546             5.0 * parall[0] + 10.1, 0.1, 0.1, 0.25, 1.5 * strt[4] + 0.5, 1.5 * strt[5], 3.2 * strt[6], 10.0 * parall[7]};
0547 
0548         double step[nsgpr] = {0.1, 0.01, 0.01, 0.001, 0.001, 1.0, 0.001, 0.01};
0549 
0550         gMinuit->SetFCN(fcnsg);
0551 
0552         double arglist[10];
0553         int ierflg = 0;
0554         arglist[0] = 0.5;
0555         gMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
0556 
0557         for (int lm = 0; lm < nsgpr; lm++) {
0558           gMinuit->mnparm(lm, namex[lm], strt[lm], step[lm], alowmn[lm], ahighmn[lm], ierflg);
0559         }
0560 
0561         arglist[0] = 0;
0562         gMinuit->mnexcm("MINIMIZE", arglist, 0, ierflg);
0563 
0564         arglist[0] = 0;
0565         gMinuit->mnexcm("IMPROVE", arglist, 0, ierflg);
0566 
0567         TString chnam;
0568         double parv, err, xlo, xup, plerr, mierr, eparab, gcc;
0569         int iuit;
0570 
0571         for (int lm = 0; lm < nsgpr; lm++) {
0572           gMinuit->mnpout(lm, chnam, parv, err, xlo, xup, iuit);
0573           gMinuit->mnerrs(lm, plerr, mierr, eparab, gcc);
0574           fitres[lm] = fitprm[lm][jk] = parv;
0575           parserr[lm] = err;
0576         }
0577 
0578         fitres[0] *= binwid;
0579         fitres[5] *= binwid;
0580 
0581         double fedm, errdef;
0582         int nparx, istat;
0583         gMinuit->mnstat(fitchisq, fedm, errdef, fitndof, nparx, istat);
0584 
0585         delete gMinuit;
0586 
0587         anormglb = 1;
0588         double tmpnorm = 0;
0589         for (int ix = 0; ix < 600; ix++) {
0590           double xval[2];
0591           xval[0] = -5 + (ix + 0.5) * .1;
0592           tmpnorm += 0.1 * langaufun(xval, &fitres[3]) / fitres[5];
0593         }
0594         anormglb = tmpnorm;
0595 
0596         double stp = 30 * fitres[3] * fitres[4] / 1000.;
0597         double str = fitres[4] * (1. - 5. * fitres[3]);
0598         double xx[2];
0599         double sum1 = 0;
0600         double sum2 = 0;
0601 
0602         for (int lm = 0; lm < 1000; lm++) {
0603           xx[0] = str + (lm + 0.5) * stp;
0604           double landf = langaufun(xx, &fitres[3]);  //No need of normalisation
0605           sum1 += landf;
0606           sum2 += xx[0] * landf;
0607         }
0608 
0609         sum2 /= TMath::Max(0.1, sum1);
0610         //      signall[izone]->GetXaxis()->SetRangeUser(-0.25, 4.75);
0611         signall[izone]->Draw("hist");
0612 
0613         sprintf(temp, "pedfun_%i", izone);
0614         pedfun[izone] = new TF1(temp, gausX, xmn, xmx, nbgpr);
0615         pedfun[izone]->SetParameters(fitres);
0616         pedfun[izone]->SetLineColor(3);
0617         pedfun[izone]->SetLineWidth(1);
0618         pedfun[izone]->Draw("same");
0619 
0620         sprintf(temp, "signalfun_%i", izone);
0621 
0622         sigfun[izone] = new TF1(temp, langaufun, xmn, xmx, nsgpr - nbgpr);
0623         sigfun[izone]->SetParameters(&fitres[3]);
0624         sigfun[izone]->SetLineWidth(1);
0625         sigfun[izone]->SetLineColor(4);
0626         sigfun[izone]->Draw("same");
0627 
0628         cout << "sum2 " << sum2 << " " << sigfun[izone]->Integral(fitres[4] * (1. - 5. * fitres[3]), str + 1000.5 * stp)
0629              << " " << binwid << endl;
0630 
0631         sprintf(temp, "total_%i", izone);
0632         signalx[izone] = new TF1(temp, totalfunc, xmn, xmx, nsgpr);
0633         signalx[izone]->SetParameters(fitres);
0634         signalx[izone]->SetLineWidth(1);
0635         signalx[izone]->Draw("same");
0636 
0637         latex.DrawLatex(
0638             0.60, 0.83, Form("#mu: %g#pm%g", int(1000 * fitres[4]) / 1000., int(1000 * parserr[4]) / 1000.));
0639         latex.DrawLatex(0.60,
0640                         0.765,
0641                         Form("#Gamma: %g#pm%g",
0642                              int(1000 * fitres[3] * fitres[4]) / 1000.,
0643                              int(1000 * (sqrt((fitres[3] * parserr[4]) * (fitres[3] * parserr[4]) +
0644                                               (fitres[4] * parserr[3]) * (fitres[4] * parserr[3])))) /
0645                                  1000.));
0646         latex.DrawLatex(
0647             0.60, 0.70, Form("#sigma: %g#pm%g", int(1000 * fitres[6]) / 1000., int(1000 * parserr[6]) / 1000.));
0648         latex.DrawLatex(0.65, 0.64, Form("A: %g#pm%g", int(1 * fitres[5]) / 1., int(10 * parserr[5]) / 10.));
0649         latex.DrawLatex(0.67, 0.58, Form("Ex: %g#pm%g", int(10 * fitres[7]) / 10., int(10 * parserr[7]) / 10.));
0650 
0651         latex.DrawLatex(0.67, 0.52, Form("Mean: %g ", int(100 * sum2) / 100.));
0652 
0653         cout << "histinfo fit " << std::setw(3) << ieta << " " << std::setw(3) << ij + 1 << " " << std::setw(5)
0654              << signallunb[izone]->GetEntries() / 2. << " " << std::setw(6) << signallunb[izone]->GetMean() << " "
0655              << std::setw(6) << signallunb[izone]->GetRMS() << " " << std::setw(6) << fitchisq << " " << std::setw(3)
0656              << fitndof << endl;
0657 
0658         file_out << "histinfo fit " << std::setw(3) << ieta << " " << std::setw(3) << ij + 1 << " " << std::setw(5)
0659                  << signallunb[izone]->GetEntries() / 2. << " " << std::setw(6) << signallunb[izone]->GetMean() << " "
0660                  << std::setw(6) << signallunb[izone]->GetRMS() << endl;
0661 
0662         file_out << "fitresx " << ieta << " " << ij + 1 << " " << fitres[0] / binwid << " " << fitres[1] << " "
0663                  << fitres[2] << " " << fitres[3] << " " << fitres[4] << " " << fitres[5] / binwid << " " << fitres[6]
0664                  << " " << signallunb[izone]->GetEntries() / 2. << " " << fitchisq << " " << fitndof << " " << binwid
0665                  << " " << sum2 << " " << fitres[7] << endl;
0666         file_out << "parserr " << ieta << " " << ij + 1 << " " << parserr[0] << " " << parserr[1] << " " << parserr[2]
0667                  << " " << parserr[3] << " " << parserr[4] << " " << parserr[5] << " " << parserr[6] << " "
0668                  << parserr[7] << endl;
0669 
0670         double diff = fitres[4] - fitres[1];
0671         if (diff <= 0)
0672           diff = 0.000001;
0673 
0674         int ifl = nphimx * jk + ij;
0675 
0676         fit_chi->Fill(ifl, fitchisq);  //signal[izone]->GetChisquare());
0677         sig_evt->Fill(ifl, signallunb[izone]->GetEntries() / 2.);
0678         fit_sigevt->Fill(ifl, fitres[5] / binwid);
0679         fit_bkgevt->Fill(ifl, fitres[0] / binwid);
0680 
0681         sig_mean->Fill(ifl, fitres[4]);
0682         sig_diff->Fill(ifl, fitres[4] - fitres[1]);
0683         sig_width->Fill(ifl, fitres[3]);
0684         sig_sigma->Fill(ifl, fitres[6]);
0685         sig_expo->Fill(ifl, fitres[7]);
0686         sig_meanerr->Fill(ifl, parserr[4]);
0687         if (fitres[4] - fitres[1] > 1.e-4)
0688           sig_meanerrp->Fill(ifl, 100 * parserr[4] / (fitres[4] - fitres[1]));
0689         if (fitres[2] > 1.e-4)
0690           sig_signf->Fill(ifl, (fitres[4] - fitres[1]) / fitres[2]);
0691 
0692         sig_statmean->Fill(ifl, signallunb[izone]->GetMean());
0693         sig_rms->Fill(ifl, signallunb[izone]->GetRMS());
0694 
0695         fit2d_chi->Fill(ieta, iphi, fitchisq);  //signal[izone]->GetChisquare());
0696         sig2d_evt->Fill(ieta, iphi, signallunb[izone]->GetEntries() / 2.);
0697         fit2d_sigevt->Fill(ieta, iphi, fitres[5] / binwid);
0698         fit2d_bkgevt->Fill(ieta, iphi, fitres[0] / binwid);
0699 
0700         sig2d_mean->Fill(ieta, iphi, fitres[4]);
0701         sig2d_diff->Fill(ieta, iphi, fitres[4] - fitres[1]);
0702         sig2d_width->Fill(ieta, iphi, fitres[3]);
0703         sig2d_sigma->Fill(ieta, iphi, fitres[6]);
0704         sig2d_expo->Fill(ieta, iphi, fitres[7]);
0705 
0706         sig2d_meanerr->Fill(ieta, iphi, parserr[4]);
0707         if (fitres[4] - fitres[1] > 1.e-4)
0708           sig2d_meanerrp->Fill(ieta, iphi, 100 * parserr[4] / (fitres[4] - fitres[1]));
0709         if (fitres[2] > 1.e-4)
0710           sig2d_signf->Fill(ieta, iphi, (fitres[4] - fitres[1]) / fitres[2]);
0711 
0712         sig2d_statmean->Fill(ieta, iphi, signallunb[izone]->GetMean());
0713         sig2d_rms->Fill(ieta, iphi, signallunb[izone]->GetRMS());
0714       } else {  //if (signallunb[izone]->GetEntries()/2. >10) {
0715         signall[izone]->Draw();
0716         float varx = 0.000;
0717 
0718         file_out << "histinfo nof " << std::setw(3) << ieta << " " << std::setw(3) << ij + 1 << " " << std::setw(5)
0719                  << signallunb[izone]->GetEntries() / 2. << " " << std::setw(6) << signallunb[izone]->GetMean() << " "
0720                  << std::setw(6) << signallunb[izone]->GetRMS() << endl;
0721 
0722         file_out << "fitresx " << ieta << " " << ij + 1 << " " << varx << " " << varx << " " << varx << " " << varx
0723                  << " " << varx << " " << varx << " " << varx << " " << varx << " " << varx << " " << varx << " "
0724                  << varx << " " << varx << " " << varx << endl;
0725         file_out << "parserr " << ieta << " " << ij + 1 << " " << varx << " " << varx << " " << varx << " " << varx
0726                  << " " << varx << " " << varx << " " << varx << " " << varx << " " << varx << endl;
0727       }
0728       iiter++;
0729       if (iiter % nsample == 0) {
0730         c0->Update();
0731 
0732         for (int lm = 0; lm < nsample; lm++) {
0733           if (pedfun[lm]) {
0734             delete pedfun[lm];
0735             pedfun[lm] = 0;
0736           }
0737           if (sigfun[lm]) {
0738             delete sigfun[lm];
0739             sigfun[lm] = 0;
0740           }
0741           if (signalx[lm]) {
0742             delete signalx[lm];
0743             signalx[lm] = 0;
0744           }
0745           if (signall[lm]) {
0746             delete signall[lm];
0747             signall[lm] = 0;
0748           }
0749           if (signallunb[lm]) {
0750             delete signallunb[lm];
0751             signallunb[lm] = 0;
0752           }
0753         }
0754       }
0755     }  //for (int jk=0; jk<netamx; jk++) {
0756   }    //for (int ij=0; ij<nphimx; ij++) {
0757 
0758   if (iiter % nsample != 0) {
0759     c0->Update();
0760     for (int lm = 0; lm < nsample; lm++) {
0761       if (pedfun[lm]) {
0762         delete pedfun[lm];
0763         pedfun[lm] = 0;
0764       }
0765       if (sigfun[lm]) {
0766         delete sigfun[lm];
0767         sigfun[lm] = 0;
0768       }
0769       if (signalx[lm]) {
0770         delete signalx[lm];
0771         signalx[lm] = 0;
0772       }
0773       if (signall[lm]) {
0774         delete signall[lm];
0775         signall[lm] = 0;
0776       }
0777       if (signallunb[lm]) {
0778         delete signallunb[lm];
0779         signallunb[lm] = 0;
0780       }
0781     }
0782   }
0783 
0784   if (c0) {
0785     delete c0;
0786     c0 = 0;
0787   }
0788 
0789   xsiz = 600;  //int xsiz = 600;
0790   ysiz = 800;  //int ysiz = 800;
0791 
0792   ps.NewPage();
0793   gStyle->SetOptLogy(0);
0794   gStyle->SetTitleFontSize(0.05);
0795   gStyle->SetTitleSize(0.025, "XYZ");
0796   gStyle->SetLabelSize(0.025, "XYZ");
0797   gStyle->SetStatFontSize(.045);
0798 
0799   gStyle->SetPadGridX(1);  //bool input yes/no, must give an input
0800   gStyle->SetPadGridY(1);
0801   gStyle->SetGridStyle(3);
0802 
0803   gStyle->SetOptStat(0);
0804   gStyle->SetPadTopMargin(.07);
0805   gStyle->SetPadLeftMargin(0.07);
0806 
0807   ps.NewPage();
0808   TCanvas* c1 = new TCanvas("c1", " Pedestal vs signal", xsiz, ysiz);
0809   sig_evt->Draw("hist");
0810   c1->Update();
0811 
0812   ps.NewPage();
0813   sig_statmean->Draw("hist");
0814   c1->Update();
0815 
0816   ps.NewPage();
0817   sig_rms->Draw("hist");
0818   c1->Update();
0819 
0820   ps.NewPage();
0821   fit_chi->Draw("hist");
0822   c1->Update();
0823 
0824   ps.NewPage();
0825   fit_sigevt->Draw("hist");
0826   c1->Update();
0827 
0828   ps.NewPage();
0829   fit_bkgevt->Draw("hist");
0830   c1->Update();
0831 
0832   ps.NewPage();
0833   sig_mean->Draw("hist");
0834   c1->Update();
0835 
0836   ps.NewPage();
0837   sig_width->Draw("hist");
0838   c1->Update();
0839 
0840   ps.NewPage();
0841   sig_sigma->Draw("hist");
0842   c1->Update();
0843 
0844   ps.NewPage();
0845   sig_expo->Draw("hist");
0846   c1->Update();
0847 
0848   ps.NewPage();
0849   sig_meanerr->Draw("hist");
0850   c1->Update();
0851 
0852   ps.NewPage();
0853   sig_meanerrp->Draw("hist");
0854   c1->Update();
0855 
0856   ps.NewPage();
0857   sig_signf->Draw("hist");
0858   c1->Update();
0859 
0860   gStyle->SetPadLeftMargin(0.06);
0861   gStyle->SetPadRightMargin(0.15);
0862 
0863   ps.NewPage();
0864   TCanvas* c2y = new TCanvas("c2y", " Pedestal vs Signal", xsiz, ysiz);
0865 
0866   sig2d_evt->Draw("colz");
0867   c2y->Update();
0868 
0869   ps.NewPage();
0870   sig2d_statmean->SetMaximum(min(3.0, sig2d_statmean->GetMaximum()));
0871   sig2d_statmean->SetMinimum(max(-2.0, sig2d_statmean->GetMinimum()));
0872   sig2d_statmean->Draw("colz");
0873   c2y->Update();
0874 
0875   ps.NewPage();
0876   sig2d_rms->SetMaximum(min(3.0, sig2d_rms->GetMaximum()));
0877   sig2d_rms->SetMinimum(max(0.0, sig2d_rms->GetMinimum()));
0878   sig2d_rms->Draw("colz");
0879   c2y->Update();
0880 
0881   ps.NewPage();
0882 
0883   fit2d_chi->Draw("colz");
0884   c2y->Update();
0885 
0886   ps.NewPage();
0887   fit2d_sigevt->Draw("colz");
0888   c2y->Update();
0889 
0890   ps.NewPage();
0891   fit2d_bkgevt->Draw("colz");
0892   c2y->Update();
0893 
0894   ps.NewPage();
0895   sig2d_mean->SetMaximum(min(2.0, sig2d_mean->GetMaximum()));
0896   sig2d_mean->SetMinimum(max(0.1, sig2d_mean->GetMinimum()));
0897   sig2d_mean->Draw("colz");
0898   c2y->Update();
0899 
0900   ps.NewPage();
0901   sig2d_width->SetMaximum(min(0.5, sig2d_width->GetMaximum()));
0902   sig2d_width->SetMinimum(max(0.01, sig2d_width->GetMinimum()));
0903   sig2d_width->Draw("colz");
0904   c2y->Update();
0905 
0906   ps.NewPage();
0907   sig2d_sigma->SetMaximum(min(0.5, sig2d_sigma->GetMaximum()));
0908   sig2d_sigma->SetMinimum(max(0.01, sig2d_sigma->GetMinimum()));
0909   sig2d_sigma->Draw("colz");
0910   c2y->Update();
0911 
0912   ps.NewPage();
0913   sig2d_expo->SetMaximum(min(20., sig2d_expo->GetMaximum()));
0914   sig2d_expo->SetMinimum(max(2.0, sig2d_expo->GetMinimum()));
0915   sig2d_expo->Draw("colz");
0916   c2y->Update();
0917 
0918   ps.NewPage();
0919   sig2d_meanerr->Draw("colz");
0920   c2y->Update();
0921 
0922   ps.NewPage();
0923   sig2d_meanerrp->Draw("colz");
0924   c2y->Update();
0925 
0926   ps.NewPage();
0927   sig2d_signf->Draw("colz");
0928   c2y->Update();
0929 
0930   ps.Close();
0931 
0932   delete c1;
0933   delete c2y;
0934 
0935   file_out.close();
0936 
0937   fileOut->cd();
0938   fileOut->Write();
0939   fileOut->Close();
0940 
0941   fileIn->cd();
0942   fileIn->Close();
0943 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!