File indexing completed on 2024-04-06 11:58:55
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
0046
0047
0048
0049
0050
0051
0052
0053 Double_t langaufun(Double_t* x, Double_t* par) {
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066 Double_t invsq2pi = 0.3989422804014;
0067 Double_t mpshift = -0.22278298;
0068
0069
0070 Double_t np = 100.0;
0071 Double_t sc = 5.0;
0072
0073
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
0082 mpc = par[1] - mpshift * par[0] * par[1];
0083 double scale = 1;
0084 double scale2 = anormglb;
0085 if (scale2 < .1)
0086 scale2 = 0.1;
0087
0088
0089 xlow = x[0] - sc * scale * par[3];
0090 xupp = x[0] + sc * scale * par[3];
0091
0092 step = (xupp - xlow) / np;
0093
0094
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);
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);
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]);
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
0193
0194 TFile* fileIn = new TFile(outfilx, "read");
0195
0196
0197
0198 unsigned ievt, hoflag;
0199 int irun, ilumi, nprim, isect, isect2, ndof, nmuon;
0200
0201 float pileup, 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("pileup", &pileup);
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);
0260
0261
0262 const int nsample = 16;
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
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);
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
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);
0444 gStyle->SetPadBottomMargin(0.12);
0445 gStyle->SetPadLeftMargin(0.17);
0446 gStyle->SetPadRightMargin(0.03);
0447
0448 gStyle->SetOptStat(0);
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
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);
0473 gStyle->SetPadGridY(0);
0474 gStyle->SetGridStyle(3);
0475 gStyle->SetNdivisions(101, "XY");
0476 gStyle->SetOptLogy(1);
0477 int iiter = 0;
0478
0479 ps.NewPage();
0480
0481 int xsiz = 900;
0482 int ysiz = 1200;
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
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) {
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);
0524 parall[1] = 0.00;
0525 parall[2] = 0.03;
0526 parall[3] = 0.135;
0527 parall[4] = 0.7 * signallunb[izone]->GetMean();
0528 parall[5] = signallunb[izone]->GetEntries() / 2.;
0529 parall[6] = 0.2238;
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]);
0605 sum1 += landf;
0606 sum2 += xx[0] * landf;
0607 }
0608
0609 sum2 /= TMath::Max(0.1, sum1);
0610
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);
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);
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 {
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 }
0756 }
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;
0790 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);
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 }