CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HOCalibAnalyzer.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: HOCalibAnalyzer
4 // Class: HOCalibAnalyzer
5 //
13 //
14 // Original Author: Gobinda Majumder
15 // Created: Sat Jul 7 09:51:31 CEST 2007
16 //
17 //
18 
19 
20 // system include files
21 #include <memory>
22 
23 // user include files
24 #include "FWCore/Framework/interface/Frameworkfwd.h"
25 #include "FWCore/Framework/interface/EDAnalyzer.h"
26 
27 #include "FWCore/Framework/interface/Event.h"
28 #include "FWCore/Framework/interface/MakerMacros.h"
29 
30 #include "FWCore/ParameterSet/interface/ParameterSet.h"
31 //#include "DataFormats/HOCalibHit/interface/HOCalibVariables.h"
32 #include "DataFormats/HcalCalibObjects/interface/HOCalibVariables.h"
33 
34 #include "FWCore/ServiceRegistry/interface/Service.h"
35 #include "FWCore/Utilities/interface/InputTag.h"
36 #include "CommonTools/UtilAlgos/interface/TFileService.h"
37 
38 #include "TMath.h"
39 #include "TFile.h"
40 #include "TH1F.h"
41 #include "TH2F.h"
42 #include "TTree.h"
43 #include "TProfile.h"
44 #include "TPostScript.h"
45 #include "TCanvas.h"
46 #include "TF1.h"
47 #include "TStyle.h"
48 #include "TMinuit.h"
49 #include "TMath.h"
50 
51 #include <string>
52 
53 #include <iostream>
54 #include <fstream>
55 #include <iomanip>
56 //#include <sstream>
57 
58 
59 using namespace std;
60 using namespace edm;
61 
62 
63  //
64  // Look for nearby pixel through eta, phi informations for pixel cross-talk
65  // 1. Look PIXEL code from (eta,phi)
66  // 2. Go to nearby pixel code
67  // 3. Come back to (eta,phi) from pixel code
68  // Though it works, it is a very ugly/crude way to get cross talk, need better algorithms
69  //
70 
71 static const int mapx1[6][3]={{1,4,8}, {12,7,3}, {5,9,13}, {11,6,2}, {16,15,14}, {19,18,17}};
72 
73 static const int mapx2[6][3]={{1,4,8}, {12,7,3}, {5,9,13}, {11,6,2}, {16,15,14}, {-1,-1,-1}};
74 
75 static const int mapx0p[9][2]={{3,1}, {7,4}, {6,5}, {12,8}, {0,0}, {11,9}, {16,13}, {15,14}, {19,17}};
76 static const int mapx0m[9][2]={{17,19}, {14,15}, {13,16}, {9,11}, {0,0}, {8,12}, {5,6}, {4,7}, {1,3}};
77 
78 static const int etamap[4][21]={{-1, 0,3,1, 0,2,3, 1,0,2, -1, 3,1,2, 4,4,4, -1,-1,-1, -1}, //etamap2
79  {-1, 0,3,1, 0,2,3, 1,0,2, -1, 3,1,2, 4,4,4, 5,5,5, -1}, //etamap1
80  {-1, 0,-1,0, 1,2,2, 1,3,5, -1, 5,3,6, 7,7,6, 8,-1,8, -1}, //etamap0p
81  {-1, 8,-1,8, 7,6,6, 7,5,3, -1, 3,5,2, 1,1,2, 0,-1,0, -1}}; //etamap0m
82 
83 static const int phimap[4][21] ={{-1, 0,2,2, 1,0,1, 1,2,1, -1, 0,0,2, 2,1,0, 2,1,0, -1}, //phimap2
84  {-1, 0,2,2, 1,0,1, 1,2,1, -1, 0,0,2, 2,1,0, 2,1,0, -1}, //phimap1
85  {-1, 1,-1,0, 1,1,0, 0,1,1, -1, 0,0,1, 1,0,0, 1,-1,0, -1}, //phimap0p
86  {-1, 0,-1,1, 0,0,1, 1,0,0, -1, 1,1,0, 0,1,1, 0,-1,1, -1}}; //phimap0m
87 //swapped phi map for R0+/R0- (15/03/07)
88 
89 static const int npixleft[21]={0, 0, 1, 2, 0, 4, 5, 6, 0, 8, 0, 0,11, 0,13,14,15, 0,17,18,0};
90 static const int npixrigh[21]={0, 2, 3, 0, 5, 6, 7, 0, 9, 0, 0,12, 0,14,15,16, 0,18,19, 0,0};
91 static const int npixlebt[21]={0, 0, 0, 0, 0, 1, 2, 3, 0, 4, 0, 6, 7, 8, 9, 0,11,13,14,15,0};
92 static const int npixribt[21]={0, 0, 0, 0, 1, 2, 3, 0, 4, 5, 0, 7, 0, 9, 0,11,12,14,15,16,0};
93 static const int npixleup[21]={0, 4, 5, 6, 8, 9, 0,11, 0,13, 0,15,16, 0,17,18,19, 0, 0, 0,0};
94 static const int npixriup[21]={0, 5, 6, 7, 9, 0,11,12,13,14, 0,16, 0,17,18,19, 0, 0, 0, 0,0};
95 
96 static const int netamx=30;
97 static const int nphimx =72;
98  const int nbgpr = 3;
99  const int nsgpr = 7;
100 
101 int ietafit;
102 int iphifit;
103 vector<float>sig_reg[netamx][nphimx+1];
104 vector<float>cro_ssg[netamx][nphimx+1];
105 
106 
107 //#define CORREL
108 //
109 // class decleration
110 //
111 
112 Double_t gausX(Double_t* x, Double_t* par){
113  return par[0]*(TMath::Gaus(x[0], par[1], par[2], kTRUE));
114 }
115 
116 Double_t langaufun(Double_t *x, Double_t *par) {
117 
118  //Fit parameters:
119  //par[0]*par[1]=Width (scale) parameter of Landau density
120  //par[1]=Most Probable (MP, location) parameter of Landau density
121  //par[2]=Total area (integral -inf to inf, normalization constant)
122  //par[3]=Width (sigma) of convoluted Gaussian function
123  //
124  //In the Landau distribution (represented by the CERNLIB approximation),
125  //the maximum is located at x=-0.22278298 with the location parameter=0.
126  //This shift is corrected within this function, so that the actual
127  //maximum is identical to the MP parameter.
128  // /*
129  // Numeric constants
130  Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
131  Double_t mpshift = -0.22278298; // Landau maximum location
132 
133  // Control constants
134  Double_t np = 100.0; // number of convolution steps
135  Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas
136 
137  // Variables
138  Double_t xx;
139  Double_t mpc;
140  Double_t fland;
141  Double_t sum = 0.0;
142  Double_t xlow,xupp;
143  Double_t step;
144  Double_t i;
145 
146 
147  // MP shift correction
148  mpc = par[1] - mpshift * par[0]*par[1];
149 
150  // Range of convolution integral
151  xlow = x[0] - sc * par[3];
152  xupp = x[0] + sc * par[3];
153 
154  step = (xupp-xlow) / np;
155 
156  // Convolution integral of Landau and Gaussian by sum
157  for(i=1.0; i<=np/2; i++) {
158  xx = xlow + (i-.5) * step;
159  fland = TMath::Landau(xx,mpc,par[0]*par[1], kTRUE); // / par[0];
160  sum += fland * TMath::Gaus(x[0],xx,par[3]);
161  xx = xupp - (i-.5) * step;
162  fland = TMath::Landau(xx,mpc,par[0]*par[1], kTRUE); // / par[0];
163  sum += fland * TMath::Gaus(x[0],xx,par[3]);
164  }
165 
166  return (par[2] * step * sum * invsq2pi / par[3]);
167 }
168 
169 Double_t totalfunc(Double_t* x, Double_t* par){
170  return gausX(x, par) + langaufun(x, &par[3]);
171 }
172 
173 void fcnbg(Int_t &npar, Double_t* gin, Double_t &f, Double_t* par, Int_t flag) {
174 
175  double fval = -par[0];
176  for (unsigned i=0; i<cro_ssg[ietafit][iphifit].size(); i++) {
177  double xval = (double)cro_ssg[ietafit][iphifit][i];
178  fval +=log(max(1.e-30,par[0]*TMath::Gaus(xval, par[1], par[2], 1)));
179  // fval +=log(par[0]*TMath::Gaus(xval, par[1], par[2], 1));
180  }
181  f = -fval;
182 }
183 
184 void fcnsg(Int_t &npar, Double_t* gin, Double_t &f, Double_t* par, Int_t flag) {
185 
186  double xval[2];
187  double fval = -(par[0]+par[5]);
188  for (unsigned i=0; i<sig_reg[ietafit][iphifit].size(); i++) {
189  xval[0] = (double)sig_reg[ietafit][iphifit][i];
190  fval +=log(totalfunc(xval, par));
191  }
192  f = -fval;
193 }
194 
195 void set_mean(double& x, bool mdigi) {
196  if(mdigi) {
197  x = min(x, 0.5);
198  x = max(x, -0.5);
199  } else {
200  x = min(x, 0.1);
201  x = max(x, -0.1);
202  }
203 }
204 
205 void set_sigma(double& x, bool mdigi) {
206  if(mdigi) {
207  x = min(x, 1.2);
208  x = max(x, -1.2);
209  } else {
210  x = min(x, 0.24);
211  x = max(x, -0.24);
212  }
213 }
214 
215 
216 
217 class HOCalibAnalyzer : public edm::EDAnalyzer {
218  public:
219  explicit HOCalibAnalyzer(const edm::ParameterSet&);
220  ~HOCalibAnalyzer();
221 
222 
223  private:
224 
225  virtual void beginJob() override ;
226  virtual void analyze(const edm::Event&, const edm::EventSetup&) override;
227  virtual void endJob() override ;
228 
229  TFile* theFile;
230  std::string theRootFileName;
231  std::string theoutputtxtFile;
232  std::string theoutputpsFile;
233 
234  bool m_hotime;
235  bool m_hbtime;
236  bool m_correl;
237  bool m_checkmap;
238  bool m_hbinfo;
239  bool m_combined;
240  bool m_constant;
241  bool m_figure;
242  bool m_digiInput;
243  bool m_cosmic;
244  bool m_histfit;
245  bool m_pedsuppr;
246  double m_sigma;
247 
248  static const int ncut = 13;
249  static const int mypow_2_ncut = 8192; // 2^13, should be changed to match ncut
250 
251  int ipass;
252 
253  TTree* T1;
254 
255  TH1F* muonnm;
256  TH1F* muonmm;
257  TH1F* muonth;
258  TH1F* muonph;
259  TH1F* muonch;
260 
261  TH1F* sel_muonnm;
262  TH1F* sel_muonmm;
263  TH1F* sel_muonth;
264  TH1F* sel_muonph;
265  TH1F* sel_muonch;
266 
267  // if (m_hotime) { // #ifdef HOTIME
268  TProfile* hotime[netamx][nphimx];
269  TProfile* hopedtime[netamx][nphimx];
270  // } // if (m_hotime) #endif
271  // if (m_hbtime) { // #ifdef HBTIME
272  TProfile* hbtime[netamx][nphimx];
273  // } // m_hbtime #endif
274 
275  // double pedval[netamx][nphimx]; // ={0};
276 
277  // if (m_correl) { // #ifdef CORREL
278  TH1F* corrsglb[netamx][nphimx];
279  TH1F* corrsgrb[netamx][nphimx];
280  TH1F* corrsglu[netamx][nphimx];
281  TH1F* corrsgru[netamx][nphimx];
282  TH1F* corrsgall[netamx][nphimx];
283 
284  TH1F* corrsgl[netamx][nphimx];
285  TH1F* corrsgr[netamx][nphimx];
286 
287  TH1F* mncorrsglb;
288  TH1F* mncorrsgrb;
289  TH1F* mncorrsglu;
290  TH1F* mncorrsgru;
291  TH1F* mncorrsgall;
292 
293  TH1F* mncorrsgl;
294  TH1F* mncorrsgr;
295 
296  TH1F* rmscorrsglb;
297  TH1F* rmscorrsgrb;
298  TH1F* rmscorrsglu;
299  TH1F* rmscorrsgru;
300  TH1F* rmscorrsgall;
301 
302  TH1F* rmscorrsgl;
303  TH1F* rmscorrsgr;
304 
305  TH1F* nevcorrsglb;
306  TH1F* nevcorrsgrb;
307  TH1F* nevcorrsglu;
308  TH1F* nevcorrsgru;
309  TH1F* nevcorrsgall;
310 
311  TH1F* nevcorrsgl;
312  TH1F* nevcorrsgr;
313 
314  // } //m_correl #endif
315  // if (m_checkmap) { // #ifdef CHECKMAP
316  TH1F* corrsgc[netamx][nphimx];
317  TH1F* mncorrsgc;
318  TH1F* rmscorrsgc;
319  TH1F* nevcorrsgc;
320 
321  // } //m_checkmap #endif
322 
323  TH1F* sigrsg[netamx][nphimx+1];
324  TH1F* crossg[netamx][nphimx+1];
325  float invang[netamx][nphimx+1];
326 
327  TH1F* mnsigrsg;
328  TH1F* mncrossg;
329 
330  TH1F* rmssigrsg;
331  TH1F* rmscrossg;
332 
333  TH1F* nevsigrsg;
334  TH1F* nevcrossg;
335 
336  TH1F* ho_sig2p[9];
337  TH1F* ho_sig1p[9];
338  TH1F* ho_sig00[9];
339  TH1F* ho_sig1m[9];
340  TH1F* ho_sig2m[9];
341 
342  // if (m_hbinfo) { // #ifdef HBINFO
343  TH1F* hbhe_sig[9];
344  // } // m_hbinfo #endif
345 
346 
347  // if (m_combined) { //#ifdef COMBINED
348  static const int ringmx=5;
349  static const int sectmx=12;
350  static const int routmx=36;
351  static const int rout12mx=24;
352  static const int neffip=6;
353 
354  // if (m_hotime) { // #ifdef HOTIME
355  TProfile* com_hotime[ringmx][sectmx];
356  TProfile* com_hopedtime[ringmx][sectmx];
357  // } //m_hotime #endif
358  // if (m_hbtime) { #ifdef HBTIME
359  TProfile* com_hbtime[ringmx][sectmx];
360  // } // m_hbtime #endif
361  // if (m_correl) { //#ifdef CORREL
362  TH1F* com_corrsglb[ringmx][sectmx];
363  TH1F* com_corrsgrb[ringmx][sectmx];
364  TH1F* com_corrsglu[ringmx][sectmx];
365  TH1F* com_corrsgru[ringmx][sectmx];
366  TH1F* com_corrsgall[ringmx][sectmx];
367 
368  TH1F* com_corrsgl[ringmx][sectmx];
369  TH1F* com_corrsgr[ringmx][sectmx];
370  // } //m_correl #endif
371  // if (m_checkmap) { #ifdef CHECKMAP
372  TH1F* com_corrsgc[ringmx][sectmx];
373  // } // m_checkmap #endif
374 
375  TH1F* com_sigrsg[ringmx][routmx+1];
376  TH1F* com_crossg[ringmx][routmx+1];
377  float com_invang[ringmx][routmx+1];
378 
379 
380  // } //m_combined #endif
381 
382  TH1F* ped_evt;
383  TH1F* ped_mean;
384  TH1F* ped_width;
385  TH1F* fit_chi;
386  TH1F* sig_evt;
387  TH1F* fit_sigevt;
388  TH1F* fit_bkgevt;
389  TH1F* sig_mean;
390  TH1F* sig_diff;
391  TH1F* sig_width;
392  TH1F* sig_sigma;
393  TH1F* sig_meanerr;
394  TH1F* sig_meanerrp;
395  TH1F* sig_signf;
396 
397  TH1F* ped_statmean;
398  TH1F* sig_statmean;
399  TH1F* ped_rms;
400  TH1F* sig_rms;
401 
402  TH2F* const_eta_phi;
403 
404  TH1F* const_eta[netamx];
405  TH1F* stat_eta[netamx];
406  TH1F* statmn_eta[netamx];
407  TH1F* peak_eta[netamx];
408 
409  TH1F* const_hpdrm[ringmx];
410  // TH1F* stat_hpdrm[ringmx];
411  // TH1F* statmn_hpdrm[ringmx];
412  TH1F* peak_hpdrm[ringmx];
413 
414  TH1F* mean_eta_ave;
415  TH1F* mean_phi_ave;
416  TH1F* mean_phi_hst;
417 
418  TH2F* sig_effi[neffip];
419  TH2F* mean_energy;
420 
421  double fitprm[nsgpr][netamx];
422 
423 
424  TProfile* sigvsevt[15][ncut];
425 
426 
427  // int irun, ievt, itrg1, itrg2, isect, nrecht, nfound, nlost, ndof, nmuon;
428  int irun, ievt, itrg1, itrg2, isect, ndof, nmuon;
429  float trkdr, trkdz, trkvx, trkvy, trkvz, trkmm, trkth, trkph, chisq, therr, pherr, hodx, hody,
430  hoang, htime, hosig[9], hocorsig[18], hocro, hbhesig[9], caloen[3];
431 
432  int Nevents;
433  int nbn;
434  float alow;
435  float ahigh;
436  float binwid;
437  int irunold;
438 
439  edm::InputTag hoCalibVariableCollectionTag;
440  edm::EDGetTokenT<HOCalibVariableCollection> tok_ho_;
441  // ----------member data ---------------------------
442 
443 };
444 
445 const int HOCalibAnalyzer::ringmx;
446 const int HOCalibAnalyzer::sectmx;
447 const int HOCalibAnalyzer::routmx;
448 const int HOCalibAnalyzer::rout12mx;
449 const int HOCalibAnalyzer::neffip;
450 
451 //
452 // constants, enums and typedefs
453 //
454 
455 //
456 // static data member definitions
457 //
458 
459 //
460 // constructors and destructor
461 //
462 
463 HOCalibAnalyzer::HOCalibAnalyzer(const edm::ParameterSet& iConfig)
464  // It is very likely you want the following in your configuration
465  // hoCalibVariableCollectionTag = cms.InputTag('hoCalibProducer', 'HOCalibVariableCollection')
466 {
467  tok_ho_ = consumes<HOCalibVariableCollection>(iConfig.getParameter<edm::InputTag>("hoCalibVariableCollectionTag"));
468  //now do what ever initialization is needed
469  ipass = 0;
470  Nevents = 0;
471 
472  theRootFileName = iConfig.getUntrackedParameter<string>("RootFileName", "test.root");
473  theoutputtxtFile = iConfig.getUntrackedParameter<string>("txtFileName", "test.txt");
474  theoutputpsFile = iConfig.getUntrackedParameter<string>("psFileName", "test.ps");
475 
476  m_hbinfo = iConfig.getUntrackedParameter<bool>("hbinfo", false);
477  m_hbtime = iConfig.getUntrackedParameter<bool>("hbtime", false);
478  m_hotime = iConfig.getUntrackedParameter<bool>("hotime", false);
479  m_correl = iConfig.getUntrackedParameter<bool>("correl", false);
480  m_checkmap = iConfig.getUntrackedParameter<bool>("checkmap", false);
481  m_combined = iConfig.getUntrackedParameter<bool>("combined", false);
482  m_constant = iConfig.getUntrackedParameter<bool>("get_constant", false);
483  m_figure = iConfig.getUntrackedParameter<bool>("get_figure", true);
484  m_digiInput = iConfig.getUntrackedParameter<bool>("digiInput", true);
485  m_histfit = iConfig.getUntrackedParameter<bool>("histFit", true);
486  m_pedsuppr = iConfig.getUntrackedParameter<bool>("pedSuppr", true);
487  m_cosmic = iConfig.getUntrackedParameter<bool>("cosmic", true);
488  m_sigma = iConfig.getUntrackedParameter<double>("sigma", 1.0);
489 
490  edm::Service<TFileService> fs;
491 
492  theFile = new TFile(theRootFileName.c_str(), "RECREATE");
493  theFile->cd();
494 
495  T1 = new TTree("T1", "DT+CSC+HO");
496 
497  T1->Branch("irun",&irun,"irun/I");
498  T1->Branch("ievt",&ievt,"ievt/I");
499  T1->Branch("itrg1",&itrg1,"itrg1/I");
500  T1->Branch("itrg2",&itrg2,"itrg2/I");
501 
502  T1->Branch("isect",&isect,"isect/I");
503  // T1->Branch("nrecht",&nrecht,"nrecht/I");
504  T1->Branch("ndof",&ndof,"ndof/I");
505  T1->Branch("nmuon",&nmuon,"nmuon/I");
506 
507  T1->Branch("trkdr",&trkdr,"trkdr/F");
508  T1->Branch("trkdz",&trkdz,"trkdz/F");
509 
510  T1->Branch("trkvx",&trkvx,"trkvx/F");
511  T1->Branch("trkvy",&trkvy,"trkvy/F");
512  T1->Branch("trkvz",&trkvz,"trkvz/F");
513  T1->Branch("trkmm",&trkmm,"trkmm/F");
514  T1->Branch("trkth",&trkth,"trkth/F");
515  T1->Branch("trkph",&trkph,"trkph/F");
516 
517  T1->Branch("chisq",&chisq,"chisq/F");
518  T1->Branch("therr",&therr,"therr/F");
519  T1->Branch("pherr",&pherr,"pherr/F");
520  T1->Branch("hodx",&hodx,"hodx/F");
521  T1->Branch("hody",&hody,"hody/F");
522  T1->Branch("hoang",&hoang,"hoang/F");
523 
524  T1->Branch("htime",&htime,"htime/F");
525  T1->Branch("hosig",hosig,"hosig[9]/F");
526  T1->Branch("hocro",&hocro,"hocro/F");
527  T1->Branch("hocorsig",hocorsig,"hocorsig[18]/F");
528  T1->Branch("caloen",caloen,"caloen[3]/F");
529 
530  if (m_hbinfo) { // #ifdef HBINFO
531  T1->Branch("hbhesig",hbhesig,"hbhesig[9]/F");
532  } //m_hbinfo #endif
533 
534  muonnm = fs->make<TH1F>("muonnm", "No of muon", 10, -0.5, 9.5);
535  muonmm = fs->make<TH1F>("muonmm", "P_{mu}", 200, -100., 100.);
536  muonth = fs->make<TH1F>("muonth", "{Theta}_{mu}", 180, 0., 180.);
537  muonph = fs->make<TH1F>("muonph", "{Phi}_{mu}", 180, -180., 180.);
538  muonch = fs->make<TH1F>("muonch", "{chi^2}/ndf", 100, 0., 1000.);
539 
540  sel_muonnm = fs->make<TH1F>("sel_muonnm", "No of muon(sel)", 10, -0.5, 9.5);
541  sel_muonmm = fs->make<TH1F>("sel_muonmm", "P_{mu}(sel)", 200, -100., 100.);
542  sel_muonth = fs->make<TH1F>("sel_muonth", "{Theta}_{mu}(sel)", 180, 0., 180.);
543  sel_muonph = fs->make<TH1F>("sel_muonph", "{Phi}_{mu}(sel)", 180, -180., 180.);
544  sel_muonch = fs->make<TH1F>("sel_muonch", "{chi^2}/ndf(sel)", 100, 0., 1000.);
545 
546 
547  int nbin = 50; //40;// 45; //50; //55; //60; //55; //45; //40; //50;
548  alow = -2.0;// -1.85; //-1.90; // -1.95; // -2.0;
549  ahigh = 8.0;// 8.15; // 8.10; // 8.05; // 8.0;
550 
551  if (m_digiInput) {alow = -10.0; ahigh = 40.;}
552 
553  float tmpwid = (ahigh-alow)/nbin;
554  nbn = int(-alow/tmpwid)+1;
555  if (nbn <0) nbn = 0;
556  if (nbn>nbin) nbn = nbin;
557 
558  char name[200];
559  char title[200];
560 
561  cout <<"nbin "<< nbin<<" "<<alow<<" "<<ahigh<<" "<<tmpwid<<" "<<nbn<<endl;
562 
563  for (int i=0; i<15; i++) {
564 
565  sprintf(title, "sigvsndof_ring%i", i+1);
566  sigvsevt[i][0] = fs->make<TProfile>(title, title, 50, 0., 50.,-9., 20.);
567 
568  sprintf(title, "sigvschisq_ring%i", i+1);
569  sigvsevt[i][1] = fs->make<TProfile>(title, title, 50, 0., 30.,-9., 20.);
570 
571  sprintf(title, "sigvsth_ring%i", i+1);
572  sigvsevt[i][2] = fs->make<TProfile>(title, title, 50, .7, 2.4,-9., 20.);
573 
574  sprintf(title, "sigvsph_ring%i", i+1);
575  sigvsevt[i][3] = fs->make<TProfile>(title, title, 50, -2.4, -0.7,-9., 20.);
576 
577  sprintf(title, "sigvstherr_ring%i", i+1);
578  sigvsevt[i][4] = fs->make<TProfile>(title, title, 50, 0., 0.2,-9., 20.);
579 
580  sprintf(title, "sigvspherr_ring%i", i+1);
581  sigvsevt[i][5] = fs->make<TProfile>(title, title, 50, 0., 0.2,-9., 20.);
582 
583  sprintf(title, "sigvsdircos_ring%i", i+1);
584  sigvsevt[i][6] = fs->make<TProfile>(title, title, 50, 0.5, 1.,-9., 20.);
585 
586  sprintf(title, "sigvstrkmm_ring%i", i+1);
587  sigvsevt[i][7] = fs->make<TProfile>(title, title, 50, 0., 50.,-9., 20.);
588 
589  sprintf(title, "sigvsnmuon_ring%i", i+1);
590  sigvsevt[i][8] = fs->make<TProfile>(title, title, 5, 0.5, 5.5,-9., 20.);
591 
592  sprintf(title, "sigvserr_ring%i", i+1);
593  sigvsevt[i][9] = fs->make<TProfile>(title, title, 50, 0., .3, -9., 20.);
594 
595  sprintf(title, "sigvsaccx_ring%i", i+1);
596  sigvsevt[i][10] = fs->make<TProfile>(title, title, 100, -25., 25., -9., 20.);
597 
598  sprintf(title, "sigvsaccy_ring%i", i+1);
599  sigvsevt[i][11] = fs->make<TProfile>(title, title, 100, -25., 25., -9., 20.);
600 
601  sprintf(title, "sigvscalo_ring%i", i+1);
602  sigvsevt[i][12] = fs->make<TProfile>(title, title, 100, 0., 15., -9., 20.);
603  }
604 
605  for (int j=0; j<netamx; j++) {
606  int ieta = (j<15) ? j+1 : 14-j;
607  for (int i=0;i<nphimx+1;i++) {
608  if (i==nphimx) {
609  sprintf(title, "sig_eta%i_allphi", ieta);
610  } else {
611  sprintf(title, "sig_eta%i_phi%i", ieta,i+1);
612  }
613  sigrsg[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
614  if (i==nphimx) {
615  sprintf(title, "ped_eta%i_allphi", ieta);
616  } else {
617  sprintf(title, "ped_eta%i_phi%i", ieta,i+1);
618  }
619  crossg[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
620  }
621 
622  for (int i=0;i<nphimx;i++) {
623  if (m_hotime) { //#ifdef HOTIME
624  sprintf(title, "hotime_eta%i_phi%i", (j<=14) ? j+1 : 14-j, i+1);
625  hotime[j][i] = fs->make<TProfile>(title, title, 10, -0.5, 9.5, -1.0, 30.0);
626 
627  sprintf(title, "hopedtime_eta%i_phi%i", (j<=14) ? j+1 : 14-j, i+1);
628  hopedtime[j][i] = fs->make<TProfile>(title, title, 10, -0.5, 9.5, -1.0, 30.0);
629 
630  } //m_hotime #endif
631  if (m_hbtime) { //#ifdef HBTIME
632  sprintf(title, "hbtime_eta%i_phi%i", (j<=15) ? j+1 : 15-j, i+1);
633  hbtime[j][i] = fs->make<TProfile>(title, title, 10, -0.5, 9.5, -1.0, 30.0);
634  } //m_hbtime #endif
635 
636  if (m_correl) { //#ifdef CORREL
637  sprintf(title, "corrsg_eta%i_phi%i_leftbottom", ieta,i+1);
638  corrsglb[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
639 
640  sprintf(title, "corrsg_eta%i_phi%i_rightbottom", ieta,i+1);
641  corrsgrb[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
642 
643  sprintf(title, "corrsg_eta%i_phi%i_leftup", ieta,i+1);
644  corrsglu[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
645 
646  sprintf(title, "corrsg_eta%i_phi%i_rightup", ieta,i+1);
647  corrsgru[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
648 
649  sprintf(title, "corrsg_eta%i_phi%i_all", ieta,i+1);
650  corrsgall[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
651 
652  sprintf(title, "corrsg_eta%i_phi%i_left", ieta,i+1);
653  corrsgl[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
654 
655  sprintf(title, "corrsg_eta%i_phi%i_right", ieta,i+1);
656  corrsgr[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
657  } //m_correl #endif
658  if (m_checkmap) {// #ifdef CHECKMAP
659  sprintf(title, "corrsg_eta%i_phi%i_centrl", ieta,i+1);
660  corrsgc[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
661  } //m_checkmap #endif
662  }
663  }
664 
665  mnsigrsg = fs->make<TH1F>("mnsigrsg","mnsigrsg", netamx*nphimx+ringmx*routmx, -0.5, netamx*nphimx+ringmx*routmx -0.5);
666  rmssigrsg = fs->make<TH1F>("rmssigrsg","rmssigrsg", netamx*nphimx+ringmx*routmx, -0.5, netamx*nphimx+ringmx*routmx -0.5);
667  nevsigrsg = fs->make<TH1F>("nevsigrsg","nevsigrsg", netamx*nphimx+ringmx*routmx, -0.5, netamx*nphimx+ringmx*routmx -0.5);
668 
669  mncrossg = fs->make<TH1F>("mncrossg","mncrossg", netamx*nphimx+ringmx*routmx, -0.5, netamx*nphimx+ringmx*routmx -0.5);
670  rmscrossg = fs->make<TH1F>("rmscrossg","rmscrossg", netamx*nphimx+ringmx*routmx, -0.5, netamx*nphimx+ringmx*routmx -0.5);
671  nevcrossg = fs->make<TH1F>("nevcrossg","nevcrossg", netamx*nphimx+ringmx*routmx, -0.5, netamx*nphimx+ringmx*routmx -0.5);
672 
673 
674  for (int i=0; i<neffip; i++) {
675  if (i==0) {
676  sprintf(title, "Total projected muon in tower");
677  sprintf(name, "total_evt");
678  } else {
679  sprintf(title, "Efficiency with sig >%i #sigma", i);
680  sprintf(name, "Effi_with_gt%i_sig", i);
681  }
682  sig_effi[i] = fs->make<TH2F>(name, title, netamx+1, -netamx/2-0.5, netamx/2+0.5, nphimx, 0.5, nphimx+0.5);
683  }
684 
685  sprintf(title, "Mean Energy of all towers");
686  sprintf(name, "mean_energy");
687  mean_energy = fs->make<TH2F>(name, title, netamx+1, -netamx/2-0.5, netamx/2+0.5, nphimx, 0.5, nphimx+0.5);
688 
689  if (m_correl) { //#ifdef CORREL
690  mncorrsglb = fs->make<TH1F>("mncorrsglb","mncorrsglb", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
691  rmscorrsglb = fs->make<TH1F>("rmscorrsglb","rmscorrsglb", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
692  nevcorrsglb = fs->make<TH1F>("nevcorrsglb","nevcorrsglb", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
693 
694  mncorrsgrb = fs->make<TH1F>("mncorrsgrb","mncorrsgrb", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
695  rmscorrsgrb = fs->make<TH1F>("rmscorrsgrb","rmscorrsgrb", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
696  nevcorrsgrb = fs->make<TH1F>("nevcorrsgrb","nevcorrsgrb", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
697 
698  mncorrsglu = fs->make<TH1F>("mncorrsglu","mncorrsglu", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
699  rmscorrsglu = fs->make<TH1F>("rmscorrsglu","rmscorrsglu", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
700  nevcorrsglu = fs->make<TH1F>("nevcorrsglu","nevcorrsglu", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
701 
702  mncorrsgru = fs->make<TH1F>("mncorrsgru","mncorrsgru", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
703  rmscorrsgru = fs->make<TH1F>("rmscorrsgru","rmscorrsgru", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
704  nevcorrsgru = fs->make<TH1F>("nevcorrsgru","nevcorrsgru", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
705 
706  mncorrsgall = fs->make<TH1F>("mncorrsgall","mncorrsgall", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
707  rmscorrsgall = fs->make<TH1F>("rmscorrsgall","rmscorrsgall", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
708  nevcorrsgall = fs->make<TH1F>("nevcorrsgall","nevcorrsgall", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
709 
710  mncorrsgl = fs->make<TH1F>("mncorrsgl","mncorrsgl", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
711  rmscorrsgl = fs->make<TH1F>("rmscorrsgl","rmscorrsgl", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
712  nevcorrsgl = fs->make<TH1F>("nevcorrsgl","nevcorrsgl", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
713 
714  mncorrsgr = fs->make<TH1F>("mncorrsgr","mncorrsgr", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
715  rmscorrsgr = fs->make<TH1F>("rmscorrsgr","rmscorrsgr", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
716  nevcorrsgr = fs->make<TH1F>("nevcorrsgr","nevcorrsgr", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
717  } //m_correl #endif
718 
719  if (m_checkmap) { //#ifdef CHECKMAP
720  mncorrsgc = fs->make<TH1F>("mncorrsgc","mncorrsgc", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
721  rmscorrsgc = fs->make<TH1F>("rmscorrsgc","rmscorrsgc", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
722  nevcorrsgc = fs->make<TH1F>("nevcorrsgc","nevcorrsgc", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
723  } //m_checkmap #endif
724 
725  if (m_combined) { //#ifdef COMBINED
726  for (int j=0; j<ringmx; j++) {
727 
728  for (int i=0;i<routmx+1;i++) {
729  if (j!=2 && i>rout12mx) continue;
730  int phmn = 3*i-1;
731  int phmx = 3*i+1;
732  if (j==2) {phmn = 2*i-1; phmx=2*i;}
733  if (phmn <=0) phmn = nphimx+phmn;
734  if (phmx <=0) phmx = nphimx+phmx;
735 
736  if ((j==2 && i==routmx) || (j!=2 && i==rout12mx)) {
737  sprintf(title, "sig_ring%i_allrm", j-2);
738  sprintf(name, "sig_ring%i_allrm", j-2);
739  } else {
740  sprintf(title, "sig_ring%i_phi%i-%i", j-2,phmn,phmx);
741  sprintf(name, "sig_ring%i_rout%i", j-2,i+1);
742  }
743  com_sigrsg[j][i] = fs->make<TH1F>(name, title, nbin, alow, ahigh);
744  if ((j==2 && i==routmx) || (j!=2 && i==rout12mx)) {
745  sprintf(title, "ped_ring%i_allrm", j-2);
746  sprintf(name, "ped_ring%i_allrm", j-2);
747  } else {
748  sprintf(title, "ped_ring%i_phi%i-%i", j-2,phmn, phmx);
749  sprintf(name, "ped_ring%i_rout%i", j-2,i+1);
750  }
751  com_crossg[j][i] = fs->make<TH1F>(name, title, nbin, alow, ahigh);
752  }
753 
754  for (int i=0;i<sectmx;i++) {
755  if (m_hotime) { //#ifdef HOTIME
756  sprintf(title, "com_hotime_ring%i_sect%i", j-2, i+1);
757  com_hotime[j][i] = fs->make<TProfile>(title, title, 10, -0.5, 9.5, -1.0, 30.0);
758 
759  sprintf(title, "com_hopedtime_ring%i_sect%i", j-2, i+1);
760  com_hopedtime[j][i] = fs->make<TProfile>(title, title, 10, -0.5, 9.5, -1.0, 30.0);
761  } //m_hotime #endif
762  if (m_hbtime){ //#ifdef HBTIME
763  sprintf(title, "_com_hbtime_ring%i_serrct%i", j-2, i+1);
764  com_hbtime[j][i] = fs->make<TProfile>(title, title, 10, -0.5, 9.5, -1.0, 30.0);
765  } //m_hbtime #endif
766 
767  if (m_correl) { //#ifdef CORREL
768  sprintf(title, "com_corrsg_ring%i_sect%i_leftbottom", j-2,i+1);
769  com_corrsglb[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
770 
771  sprintf(title, "com_corrsg_ring%i_sect%i_rightbottom", j-2,i+1);
772  com_corrsgrb[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
773 
774  sprintf(title, "com_corrsg_ring%i_sect%i_leftup", j-2,i+1);
775  com_corrsglu[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
776 
777  sprintf(title, "com_corrsg_ring%i_sect%i_rightup", j-2,i+1);
778  com_corrsgru[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
779 
780  sprintf(title, "com_corrsg_ring%i_sect%i_all", j-2,i+1);
781  com_corrsgall[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
782 
783  sprintf(title, "com_corrsg_ring%i_sect%i_left", j-2,i+1);
784  com_corrsgl[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
785 
786  sprintf(title, "com_corrsg_ring%i_sect%i_right", j-2,i+1);
787  com_corrsgr[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
788  } //m_correl #endif
789 
790  if (m_checkmap) { // #ifdef CHECKMAP
791  sprintf(title, "com_corrsg_ring%i_sect%i_centrl", j-2,i+1);
792  com_corrsgc[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
793  } //m_checkmap #endif
794  }
795  }
796  } //m_combined #endif
797 
798  for (int i=-1; i<=1; i++) {
799  for (int j=-1; j<=1; j++) {
800  int k = 3*(i+1)+j+1;
801 
802  sprintf(title, "hosct2p_eta%i_phi%i", i, j);
803  ho_sig2p[k] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
804 
805  sprintf(title, "hosct1p_eta%i_phi%i", i, j);
806  ho_sig1p[k] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
807 
808  sprintf(title, "hosct00_eta%i_phi%i", i, j);
809  ho_sig00[k] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
810 
811  sprintf(title, "hosct1m_eta%i_phi%i", i, j);
812  ho_sig1m[k] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
813 
814  sprintf(title, "hosct2m_eta%i_phi%i", i, j);
815  ho_sig2m[k] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
816 
817  if (m_hbinfo) { // #ifdef HBINFO
818  sprintf(title, "hbhesig_eta%i_phi%i", i, j);
819  hbhe_sig[k] = fs->make<TH1F>(title, title, 51, -10.5, 40.5);
820  } //m_hbinfo #endif
821  }
822  }
823 
824  if (m_constant) {
825  ped_evt = fs->make<TH1F>("ped_evt", "ped_evt", netamx*nphimx, -0.5, netamx*nphimx-0.5);
826  ped_mean = fs->make<TH1F>("ped_mean", "ped_mean", netamx*nphimx, -0.5, netamx*nphimx-0.5);
827  ped_width = fs->make<TH1F>("ped_width", "ped_width", netamx*nphimx, -0.5, netamx*nphimx-0.5);
828 
829  fit_chi = fs->make<TH1F>("fit_chi", "fit_chi", netamx*nphimx, -0.5, netamx*nphimx-0.5);
830  sig_evt = fs->make<TH1F>("sig_evt", "sig_evt", netamx*nphimx, -0.5, netamx*nphimx-0.5);
831  fit_sigevt = fs->make<TH1F>("fit_sigevt", "fit_sigevt", netamx*nphimx, -0.5, netamx*nphimx-0.5);
832  fit_bkgevt = fs->make<TH1F>("fit_bkgevt", "fit_bkgevt", netamx*nphimx, -0.5, netamx*nphimx-0.5);
833  sig_mean = fs->make<TH1F>("sig_mean", "sig_mean", netamx*nphimx, -0.5, netamx*nphimx-0.5);
834  sig_diff = fs->make<TH1F>("sig_diff", "sig_diff", netamx*nphimx, -0.5, netamx*nphimx-0.5);
835  sig_width = fs->make<TH1F>("sig_width", "sig_width", netamx*nphimx, -0.5, netamx*nphimx-0.5);
836  sig_sigma = fs->make<TH1F>("sig_sigma", "sig_sigma", netamx*nphimx, -0.5, netamx*nphimx-0.5);
837  sig_meanerr = fs->make<TH1F>("sig_meanerr", "sig_meanerr", netamx*nphimx, -0.5, netamx*nphimx-0.5);
838  sig_meanerrp = fs->make<TH1F>("sig_meanerrp", "sig_meanerrp", netamx*nphimx, -0.5, netamx*nphimx-0.5);
839  sig_signf = fs->make<TH1F>("sig_signf", "sig_signf", netamx*nphimx, -0.5, netamx*nphimx-0.5);
840 
841  ped_statmean = fs->make<TH1F>("ped_statmean", "ped_statmean", netamx*nphimx, -0.5, netamx*nphimx-0.5);
842  sig_statmean = fs->make<TH1F>("sig_statmean", "sig_statmean", netamx*nphimx, -0.5, netamx*nphimx-0.5);
843  ped_rms = fs->make<TH1F>("ped_rms", "ped_rms", netamx*nphimx, -0.5, netamx*nphimx-0.5);
844  sig_rms = fs->make<TH1F>("sig_rms", "sig_rms", netamx*nphimx, -0.5, netamx*nphimx-0.5);
845 
846  const_eta_phi = fs->make<TH2F>("const_eta_phi", "const_eta_phi", netamx+1, -(netamx+1)/2., (netamx+1)/2., nphimx, 0.5, nphimx+0.5);
847 
848  for (int ij=0; ij<netamx; ij++) {
849  int ieta = (ij<15) ? ij+1 : 14-ij;
850  sprintf(title, "Cont_Eta_%i", ieta);
851  const_eta[ij] = fs->make<TH1F>(title, title, nphimx, 0.5, nphimx+0.5);
852 
853  sprintf(title, "Peak_Eta_%i", ieta);
854  peak_eta[ij] = fs->make<TH1F>(title, title, nphimx, 0.5, nphimx+0.5);
855  }
856 
857  for (int ij=0; ij<ringmx; ij++) {
858  int iring = ij-2;
859  int iread = (ij==2) ? routmx : rout12mx;
860  sprintf(title, "Cont_hpdrm_%i", iring);
861  const_hpdrm[ij] = fs->make<TH1F>(title, title, iread, 0.5, iread+0.5);
862 
863  sprintf(title, "Peak_hpdrm_%i", iring);
864  peak_hpdrm[ij] = fs->make<TH1F>(title, title, iread, 0.5, iread+0.5);
865  }
866 
867  mean_phi_hst = fs->make<TH1F>("mean_phi_hst", "mean_phi_hst", netamx+1, -(netamx+1)/2., (netamx+1)/2.);
868  mean_phi_ave = fs->make<TH1F>("mean_phi_ave", "mean_phi_ave", netamx+1, -(netamx+1)/2., (netamx+1)/2.);
869 
870  mean_eta_ave = fs->make<TH1F>("mean_eta_ave", "mean_eta_ave", nphimx, 0.5, nphimx+0.5);
871 
872  } //m_constant
873 
874  for (int ij=0; ij<netamx; ij++) {
875  int ieta = (ij<15) ? ij+1 : 14-ij;
876 
877  sprintf(title, "Stat_Eta_%i", ieta);
878  stat_eta[ij] = fs->make<TH1F>(title, title, nphimx, 0.5, nphimx+0.5);
879 
880  sprintf(title, "#mu(stat)_Eta_%i", ieta);
881  statmn_eta[ij] = fs->make<TH1F>(title, title, nphimx, 0.5, nphimx+0.5);
882  }
883 
884  for (int j=0; j<netamx; j++) {
885  for (int i=0; i<nphimx; i++) {
886  invang[j][i] = 0.0;
887  }
888  }
889  for (int j=0; j<ringmx; j++) {
890  for (int i=0;i<routmx+1;i++) {
891  com_invang[j][i] = 0.0;
892  }
893  }
894 
895 }
896 
897 
898 HOCalibAnalyzer::~HOCalibAnalyzer()
899 {
900 
901  // do anything here that needs to be done at desctruction time
902  // (e.g. close files, deallocate resources etc.)
903 
904  theFile->cd();
905  theFile->Write();
906  theFile->Close();
907  cout <<" Selected events # is "<<ipass<<endl;
908 }
909 
910 
911 //
912 // member functions
913 //
914 
915 // ------------ method called to for each event ------------
916 void
917 HOCalibAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
918 {
919 
920  // calcualte these once (and avoid the pow(int,int) ambiguities for c++)
921  int mypow_2_0 = 0; // 2^0
922  int mypow_2_1 = 2; // 2^1
923  int mypow_2_2 = 4; // 2^2
924 
925  int mypow_2_3 = 8; // 2^3
926  int mypow_2_4 = 16; // 2^4
927  int mypow_2_5 = 32; // 2^5
928  int mypow_2_6 = 64; // 2^6
929  int mypow_2_7 = 128; // 2^7
930  int mypow_2_8 = 256; // 2^8
931  int mypow_2_9 = 512; // 2^9
932  int mypow_2_10 = 1024; // 2^10
933  int mypow_2_11 = 2048; // 2^11
934  int mypow_2_12 = 4096; // 2^12
935 
936  /*
937  //FIXGM Put this is initialiser
938  int mapx1[6][3]={{1,4,8}, {12,7,3}, {5,9,13}, {11,6,2}, {16,15,14}, {19,18,17}};
939  // int etamap1[21]={-1, 0,3,1, 0,2,3, 1,0,2, -1, 3,1,2, 4,4,4, 5,5,5, -1};
940  // int phimap1[21]={-1, 0,2,2, 1,0,1, 1,2,1, -1, 0,0,2, 2,1,0, 2,1,0,-1};
941 
942  int mapx2[6][3]={{1,4,8}, {12,7,3}, {5,9,13}, {11,6,2}, {16,15,14}, {-1,-1,-1}};
943  // int etamap2[21]={-1, 0,3,1, 0,2,3, 1,0,2, -1, 3,1,2, 4,4,4, -1,-1,-1, -1};
944  // int phimap2[21]={-1, 0,2,2, 1,0,1, 1,2,1, -1, 0,0,2, 2,1,0, 2, 1, 0, -1};
945 
946  int mapx0p[9][2]={{3,1}, {7,4}, {6,5}, {12,8}, {0,0}, {11,9}, {16,13}, {15,14}, {19,17}};
947  int mapx0m[9][2]={{17,19}, {14,15}, {13,16}, {9,11}, {0,0}, {8,12}, {5,6}, {4,7}, {1,3}};
948 
949  // int etamap0p[21]={-1, 0,-1,0, 1,2,2, 1,3,5, -1, 5,3,6, 7,7,6, 8,-1,8, -1};
950  // int phimap0p[21]={-1, 1,-1,0, 1,1,0, 0,1,1, -1, 0,0,1, 1,0,0, 1,-1,0, -1};
951 
952  // int etamap0m[21]={-1, 8,-1,8, 7,6,6, 7,5,3, -1, 3,5,2, 1,1,2, 0,-1,0, -1};
953  // int phimap0m[21]={-1, 0,-1,1, 0,0,1, 1,0,0, -1, 1,1,0, 0,1,1, 0,-1,1, -1};
954 
955  int etamap[4][21]={{-1, 0,3,1, 0,2,3, 1,0,2, -1, 3,1,2, 4,4,4, -1,-1,-1, -1}, //etamap2
956  {-1, 0,3,1, 0,2,3, 1,0,2, -1, 3,1,2, 4,4,4, 5,5,5, -1}, //etamap1
957  {-1, 0,-1,0, 1,2,2, 1,3,5, -1, 5,3,6, 7,7,6, 8,-1,8, -1}, //etamap0p
958  {-1, 8,-1,8, 7,6,6, 7,5,3, -1, 3,5,2, 1,1,2, 0,-1,0, -1}}; //etamap0m
959 
960  int phimap[4][21] ={{-1, 0,2,2, 1,0,1, 1,2,1, -1, 0,0,2, 2,1,0, 2,1,0, -1}, //phimap2
961  {-1, 0,2,2, 1,0,1, 1,2,1, -1, 0,0,2, 2,1,0, 2,1,0, -1}, //phimap1
962  {-1, 1,-1,0, 1,1,0, 0,1,1, -1, 0,0,1, 1,0,0, 1,-1,0, -1}, //phimap0p
963  {-1, 0,-1,1, 0,0,1, 1,0,0, -1, 1,1,0, 0,1,1, 0,-1,1, -1}}; //phimap0m
964  //swapped phi map for R0+/R0- (15/03/07)
965  for (int i=0; i<4; i++) {
966  for (int j=0; j<21; j++) {
967  cout <<"ieta "<<i<<" "<<j<<" "<<etamap[i][j]<<endl;
968  }
969  }
970 
971  // Character convention for R+/-1/2
972  // int npixleft[21] = {-1, F, Q,-1, M, D, J,-1, T,-1, C,-1, R, P, H,-1, N, G};
973  // int npixrigh[21] = { Q, S,-1, D, J, L,-1, K,-1, E,-1, P, H, B,-1, G, A,-1};
974 
975  // int npixlb1[21]={-1,-1,-1,-1, F, Q, S,-1, M, J, L, T, K,-1, C, R, P, H};
976  // int npixrb1[21]={-1,-1,-1, F, Q, S,-1, M, D, L,-1, K,-1, C, E, P, H, B};
977  // int npixlu1[21]={ M, D, J, T, K,-1, C,-1, R, H, B,-1, N, G, A,-1,-1,-1};
978  // int npixru1[21]={ D, J, L, K,-1, C, E, R, P, B,-1, N, G, A,-1,-1,-1,-1};
979 
980  int npixleft[21]={0, 0, 1, 2, 0, 4, 5, 6, 0, 8, 0, 0,11, 0,13,14,15, 0,17,18,0};
981  int npixrigh[21]={0, 2, 3, 0, 5, 6, 7, 0, 9, 0, 0,12, 0,14,15,16, 0,18,19, 0,0};
982  int npixlebt[21]={0, 0, 0, 0, 0, 1, 2, 3, 0, 4, 0, 6, 7, 8, 9, 0,11,13,14,15,0};
983  int npixribt[21]={0, 0, 0, 0, 1, 2, 3, 0, 4, 5, 0, 7, 0, 9, 0,11,12,14,15,16,0};
984  int npixleup[21]={0, 4, 5, 6, 8, 9, 0,11, 0,13, 0,15,16, 0,17,18,19, 0, 0, 0,0};
985  int npixriup[21]={0, 5, 6, 7, 9, 0,11,12,13,14, 0,16, 0,17,18,19, 0, 0, 0, 0,0};
986  */
987 
988  int iaxxx = 0;
989  int ibxxx = 0;
990 
991  Nevents++;
992 
993  using namespace edm;
994 
995  float pival = acos(-1.);
996  irunold = irun = iEvent.id().run();
997  ievt = iEvent.id().event();
998 
999  // cout <<"Nevents "<<Nevents<<endl;
1000 
1001  edm::Handle<HOCalibVariableCollection>HOCalib;
1002  bool isCosMu = true;
1003  try {
1004  iEvent.getByToken(tok_ho_, HOCalib);
1005 
1006  } catch ( cms::Exception &iEvent ) { isCosMu = false; }
1007 
1008  if (isCosMu && (*HOCalib).size() >0 ) {
1009  nmuon = (*HOCalib).size();
1010  if (Nevents%5000==1) cout <<"nmuon event # "<<Nevents<<" Run # "<<iEvent.id().run()<<" Evt # "<<iEvent.id().event()<<" "<<nmuon<<endl;
1011 
1012  for (HOCalibVariableCollection::const_iterator hoC=(*HOCalib).begin(); hoC!=(*HOCalib).end(); hoC++){
1013  itrg1 = (*hoC).trig1;
1014  itrg2 = (*hoC).trig2;
1015 
1016  trkdr = (*hoC).trkdr;
1017  trkdz = (*hoC).trkdz;
1018 
1019  trkvx = (*hoC).trkvx;
1020  trkvy = (*hoC).trkvy;
1021  trkvz = (*hoC).trkvz;
1022 
1023  trkmm = (*hoC).trkmm;
1024  trkth = (*hoC).trkth;
1025  trkph = (*hoC).trkph;
1026 
1027  ndof = (int)(*hoC).ndof;
1028  // nrecht = (int)(*hoC).nrecht;
1029  chisq = (*hoC).chisq;
1030 
1031  therr = (*hoC).therr;
1032  pherr = (*hoC).pherr;
1033  trkph = (*hoC).trkph;
1034 
1035  isect = (*hoC).isect;
1036  hodx = (*hoC).hodx;
1037  hody = (*hoC).hody;
1038  hoang = (*hoC).hoang;
1039  htime = (*hoC).htime;
1040  for (int i=0; i<9; i++) { hosig[i] = (*hoC).hosig[i];} //cout<<"hosig "<<i<<" "<<hosig[i]<<endl;}
1041  for (int i=0; i<18; i++) { hocorsig[i] = (*hoC).hocorsig[i];} // cout<<"hocorsig "<<i<<" "<<hocorsig[i]<<endl;}
1042  hocro = (*hoC).hocro;
1043  for (int i=0; i<3; i++) { caloen[i] = (*hoC).caloen[i];}
1044 
1045  if (m_hbinfo) { for (int i=0; i<9; i++) { hbhesig[i] = (*hoC).hbhesig[i];}} // cout<<"hbhesig "<<i<<" "<<hbhesig[i]<<endl;}}
1046 
1047 
1048  int ipsall=0;
1049  int ips0=0;
1050  int ips1=0;
1051  int ips2=0;
1052  int ips3=0;
1053  int ips4=0;
1054  int ips5=0;
1055  int ips6=0;
1056  int ips7=0;
1057  int ips8=0;
1058  int ips9=0;
1059  int ips10 =0;
1060  int ips11 =0;
1061  int ips12 = 0;
1062 
1063  int iselect3 = 0;
1064  if (ndof >=15 && chisq <30) iselect3 = 1;
1065 
1066  if (isect <0) continue; //FIXGM Is it proper place ?
1067  if (fabs(trkth-pival/2)<0.000001) continue; //22OCT07
1068 
1069  int ieta = int((abs(isect)%10000)/100.)-30; //an offset to acodate -ve eta values
1070  if (abs(ieta)>=16) continue;
1071  int iphi = abs(isect)%100;
1072 
1073  int tmpsect = int((iphi + 1)/6.) + 1;
1074  if (tmpsect >12) tmpsect = 1;
1075 
1076  int iring = 0;
1077  int tmpeta = ieta + 4; //For pixel mapping
1078  if (ieta >=-15 && ieta <=-11) {iring = -2; tmpeta =-11-ieta; } //abs(ieta)-11;}
1079  if (ieta >=-10 && ieta <=-5) {iring = -1; tmpeta =-5-ieta; } // abs(ieta)-5;}
1080  if (ieta >= 5 && ieta <= 10) {iring = 1; tmpeta =ieta-5; }
1081  if (ieta >= 11 && ieta <= 15) {iring = 2; tmpeta =ieta-11;}
1082 
1083  int iring2 = iring + 2;
1084 
1085  int tmprout = int((iphi + 1)/3.) + 1;
1086  int tmproutmx =routmx;
1087  if (iring==0) {
1088  tmprout = int((iphi + 1)/2.) + 1;
1089  if (tmprout >routmx) tmprout = 1;
1090  } else {
1091  if (tmprout >rout12mx) tmprout = 1;
1092  tmproutmx =rout12mx;
1093  }
1094 
1095  // CRUZET1
1096  if (m_cosmic) {
1097  /* GMA temoparily change to increase event size at 3 & 9 O'clock position */
1098  if (abs(ndof) >=20 && abs(ndof) <40) {ips0 = (int)mypow_2_0; ipsall += ips0;}
1099  if (chisq >0 && chisq<15) {ips1 = (int)mypow_2_1; ipsall +=ips1;} //18Jan2008
1100  if (fabs(trkth-pival/2) <21.5) {ips2 = (int)mypow_2_2; ipsall +=ips2;} //No nead for pp evt
1101  if (fabs(trkph+pival/2) <21.5) {ips3 = (int)mypow_2_3; ipsall +=ips3;} //No nead for pp evt
1102 
1103  if (therr <0.02) {ips4 = (int)mypow_2_4; ipsall +=ips4;}
1104  if (pherr <0.0002) {ips5 = (int)mypow_2_5; ipsall +=ips5;}
1105  if (fabs(hoang) >0.30) {ips6 = (int)mypow_2_6; ipsall +=ips6;}
1106  if (fabs(trkmm) >0.100) {ips7 = (int)mypow_2_7; ipsall +=ips7;}
1107  // if (nmuon ==1) {ips8 = (int)mypow_2_8; ipsall +=ips8;}
1108  if (nmuon >=1 && nmuon <=4) {ips8 = (int)mypow_2_8; ipsall +=ips8;}
1109 
1110  if (iring2==2) {
1111  if (fabs(hodx)<100 && fabs(hodx)>2 && fabs(hocorsig[8]) <40 && fabs(hocorsig[8]) >2 )
1112  {ips10=(int)mypow_2_10;ipsall +=ips10;}
1113 
1114  if (fabs(hody)<100 && fabs(hody)>2 && fabs(hocorsig[9]) <40 && fabs(hocorsig[9]) >2 )
1115  {ips11=(int)mypow_2_11;ipsall +=ips11;}
1116 
1117  } else {
1118  if (fabs(hodx)<100 && fabs(hodx)>2 )
1119  {ips10=(int)mypow_2_10;ipsall +=ips10;}
1120 
1121  if (fabs(hody)<100 && fabs(hody)>2)
1122  {ips11=(int)mypow_2_11;ipsall +=ips11;}
1123  }
1124  if (caloen[0] ==0) { ips12=(int)mypow_2_12;ipsall +=ips12;}
1125  } else {
1126  //csa08
1127  if (abs(ndof) >=20 && abs(ndof) <40) {ips0 = (int)mypow_2_0; ipsall += ips0;}
1128  if (chisq >0 && chisq<15) {ips1 = (int)mypow_2_1; ipsall +=ips1;} //18Jan2008
1129  if (fabs(trkth-pival/2) <21.5) {ips2 = (int)mypow_2_2; ipsall +=ips2;} //No nead for pp evt
1130  if (fabs(trkph+pival/2) <21.5) {ips3 = (int)mypow_2_3; ipsall +=ips3;} //No nead for pp evt
1131 
1132  if (therr <0.02) {ips4 = (int)mypow_2_4; ipsall +=ips4;}
1133  if (pherr <0.0002) {ips5 = (int)mypow_2_5; ipsall +=ips5;}
1134  if (fabs(hoang) >0.30) {ips6 = (int)mypow_2_6; ipsall +=ips6;}
1135  if (fabs(trkmm) >4.0) {ips7 = (int)mypow_2_7; ipsall +=ips7;}
1136  if (nmuon >=1 && nmuon <=2) {ips8 = (int)mypow_2_8; ipsall +=ips8;}
1137 
1138  if (iring2==2) {
1139  if (fabs(hodx)<100 && fabs(hodx)>2 && fabs(hocorsig[8]) <40 && fabs(hocorsig[8]) >2 )
1140  {ips10=(int)mypow_2_10;ipsall +=ips10;}
1141 
1142  if (fabs(hody)<100 && fabs(hody)>2 && fabs(hocorsig[9]) <40 && fabs(hocorsig[9]) >2 )
1143  {ips11=(int)mypow_2_11;ipsall +=ips11;}
1144 
1145  } else {
1146  if (fabs(hodx)<100 && fabs(hodx)>2 )
1147  {ips10=(int)mypow_2_10;ipsall +=ips10;}
1148 
1149  if (fabs(hody)<100 && fabs(hody)>2)
1150  {ips11=(int)mypow_2_11;ipsall +=ips11;}
1151  }
1152  // if (m_cosmic || (caloen[0] >0.5 && caloen[0]<5.0)) {ips12=(int)pow_2_12;ipsall +=ips12;}
1153  if (ndof >0 && caloen[0]<5.0) {ips12=(int)mypow_2_12;ipsall +=ips12;}
1154  /* */
1155  }
1156 
1157  if (m_digiInput) {
1158  if (htime >2.5 && htime <6.5) {ips9=(int)mypow_2_9;ipsall +=ips9;}
1159  } else {
1160  if (htime >-40 && htime <60) {ips9=(int)mypow_2_9;ipsall +=ips9;}
1161  }
1162 
1163  if (ipsall-ips0==mypow_2_ncut-mypow_2_0-1) sigvsevt[iring2][0]->Fill(abs(ndof), hosig[4]);
1164  if (ipsall-ips1==mypow_2_ncut-mypow_2_1-1) sigvsevt[iring2][1]->Fill(chisq, hosig[4]);
1165  if (ipsall-ips2==mypow_2_ncut-mypow_2_2-1) sigvsevt[iring2][2]->Fill(trkth, hosig[4]);
1166  if (ipsall-ips3==mypow_2_ncut-mypow_2_3-1) sigvsevt[iring2][3]->Fill(trkph, hosig[4]);
1167  if (ipsall-ips4==mypow_2_ncut-mypow_2_4-1) sigvsevt[iring2][4]->Fill(therr, hosig[4]);
1168  if (ipsall-ips5==mypow_2_ncut-mypow_2_5-1) sigvsevt[iring2][5]->Fill(pherr, hosig[4]);
1169  if (ipsall-ips6==mypow_2_ncut-mypow_2_6-1) sigvsevt[iring2][6]->Fill(hoang, hosig[4]);
1170  if (ipsall-ips7==mypow_2_ncut-mypow_2_7-1) sigvsevt[iring2][7]->Fill(fabs(trkmm), hosig[4]);
1171  if (ipsall-ips8==mypow_2_ncut-mypow_2_8-1) sigvsevt[iring2][8]->Fill(nmuon, hosig[4]);
1172  if (ipsall-ips9==mypow_2_ncut-mypow_2_9-1) sigvsevt[iring2][9]->Fill(htime, hosig[4]);
1173  if (ipsall-ips10==mypow_2_ncut-mypow_2_10-1) sigvsevt[iring2][10]->Fill(hodx, hosig[4]);
1174  if (ipsall-ips11==mypow_2_ncut-mypow_2_11-1) sigvsevt[iring2][11]->Fill(hody, hosig[4]);
1175  if (!m_cosmic) {
1176  if (ipsall-ips12==mypow_2_ncut-mypow_2_12-1) sigvsevt[iring2][12]->Fill(caloen[0], hosig[4]);
1177  }
1178 
1179  sigvsevt[iring2+5][0]->Fill(abs(ndof), hosig[4]);
1180  if (ips0 >0) {
1181  sigvsevt[iring2+5][1]->Fill(chisq, hosig[4]);
1182  if (ips1 >0) {
1183  sigvsevt[iring2+5][2]->Fill(trkth, hosig[4]);
1184  if (ips2 >0) {
1185  sigvsevt[iring2+5][3]->Fill(trkph, hosig[4]);
1186  if (ips3 >0) {
1187  sigvsevt[iring2+5][4]->Fill(therr, hosig[4]);
1188  if (ips4 >0) {
1189  sigvsevt[iring2+5][5]->Fill(pherr, hosig[4]);
1190  if (ips5 >0) {
1191  sigvsevt[iring2+5][6]->Fill(hoang, hosig[4]);
1192  if (ips6 >0) {
1193  sigvsevt[iring2+5][7]->Fill(fabs(trkmm), hosig[4]);
1194  if (ips7 >0) {
1195  sigvsevt[iring2+5][8]->Fill(nmuon, hosig[4]);
1196  if (ips8 >0) {
1197  sigvsevt[iring2+5][9]->Fill(htime, hosig[4]);
1198  if (ips9 >0) {
1199  sigvsevt[iring2+5][10]->Fill(hodx, hosig[4]);
1200  if (ips10>0) {
1201  sigvsevt[iring2+5][11]->Fill(hody, hosig[4]);
1202  if (ips11>0) {
1203  if (!m_cosmic) sigvsevt[iring2+5][12]->Fill(caloen[0], hosig[4]);
1204  }
1205  }
1206  }
1207  }
1208  }
1209  }
1210  }
1211  }
1212  }
1213  }
1214  }
1215  }
1216 
1217  sigvsevt[iring2+10][0]->Fill(abs(ndof), hosig[4]);
1218  sigvsevt[iring2+10][1]->Fill(chisq, hosig[4]);
1219  sigvsevt[iring2+10][2]->Fill(trkth, hosig[4]);
1220  sigvsevt[iring2+10][3]->Fill(trkph, hosig[4]);
1221  sigvsevt[iring2+10][4]->Fill(therr, hosig[4]);
1222  sigvsevt[iring2+10][5]->Fill(pherr, hosig[4]);
1223  sigvsevt[iring2+10][6]->Fill(hoang, hosig[4]);
1224  sigvsevt[iring2+10][7]->Fill(fabs(trkmm), hosig[4]);
1225  sigvsevt[iring2+10][8]->Fill(nmuon, hosig[4]);
1226  sigvsevt[iring2+10][9]->Fill(htime, hosig[4]);
1227  sigvsevt[iring2+10][10]->Fill(hodx, hosig[4]);
1228  sigvsevt[iring2+10][11]->Fill(hody, hosig[4]);
1229  if (!m_cosmic) sigvsevt[iring2+10][12]->Fill(caloen[0], hosig[4]);
1230 
1231  int iselect = (ipsall == mypow_2_ncut - 1) ? 1 : 0;
1232 
1233  if (hocro !=-100.0 && hocro <-50.0) hocro +=100.;
1234 
1235  muonnm->Fill(nmuon);
1236  muonmm->Fill(trkmm);
1237  muonth->Fill(trkth*180/pival);
1238  muonph->Fill(trkph*180/pival);
1239  muonch->Fill(chisq);
1240 
1241  if (iselect==1) {
1242  ipass++;
1243  sel_muonnm->Fill(nmuon);
1244  sel_muonmm->Fill(trkmm);
1245  sel_muonth->Fill(trkth*180/pival);
1246  sel_muonph->Fill(trkph*180/pival);
1247  sel_muonch->Fill(chisq);
1248  }
1249 
1250  if (iselect3) T1->Fill();
1251 
1252  int tmpphi = (iphi + 1)%3; //pixel mapping
1253  int npixel = 0;
1254  int itag = -1;
1255  int iflip = 0;
1256  int fact = 2;
1257 
1258  if (iring ==0) {
1259  tmpphi = (iphi+1)%2;
1260  if (tmpsect==2 || tmpsect==3 || tmpsect==6 || tmpsect==7 || tmpsect==10 || tmpsect==11) {
1261  npixel = mapx0p[tmpeta][tmpphi];
1262  itag = 2;
1263  } else {
1264  npixel = mapx0m[tmpeta][tmpphi];
1265  itag = 3;
1266  }
1267  } else {
1268  fact = 3;
1269  if (tmpsect%2==1) iflip =1;
1270  if (abs(iring)==1) {
1271  npixel = mapx1[tmpeta][(iflip==0) ? tmpphi : abs(tmpphi-2)];
1272  itag = 1;
1273  } else {
1274  npixel = mapx2[tmpeta][(iflip==0) ? tmpphi : abs(tmpphi-2)];
1275  itag = 0;
1276  }
1277  }
1278 
1279  int tmpeta1 = (ieta>0) ? ieta -1 : -ieta +14;
1280 
1281  //Histogram filling for noise study: phi shift according to DTChamberAnalysis
1282  int tmpphi1 = iphi -1 ; //(iphi+6 <=nphimx) ? iphi+5 : iphi+5-nphimx;
1283 
1284  int iselect2=0;
1285  if (hosig[4]!=-100) {
1286  if (m_cosmic) {
1287  if (caloen[2]<=0.0) iselect2=1;
1288  } else {
1289  if (caloen[2]<=3.0) iselect2=1;
1290  }
1291  }
1292 
1293  // cout <<"cosmic "<<hosig[4]<<" "<<caloen[3]<<" "<<int(iselect2)<<" "<<int(m_cosmic)<<endl;
1294 
1295 
1296  if (iselect2==1) {
1297  int tmpphi2 = iphi - 1;
1298  if (!m_digiInput) { tmpphi2 = (iphi+6 <=nphimx) ? iphi+5 : iphi+5-nphimx;}
1299 
1300  int tmpsect2 = int((tmpphi2 + 2)/6.) + 1;
1301  if (tmpsect2 >12) tmpsect2 = 1;
1302 
1303  int tmprout2 = int((tmpphi2 + 2)/3.) + 1;
1304  if (iring==0) {
1305  tmprout2 = int((tmpphi2 + 2)/2.) + 1;
1306  if (tmprout2 >routmx) tmprout2 = 1;
1307  } else {
1308  if (tmprout2 >rout12mx) tmprout2 = 1;
1309  }
1310 
1311  if (cro_ssg[tmpeta1][tmpphi2].size()<4000) {
1312  if (hocro>alow && hocro<ahigh) {
1313  if (!m_histfit) cro_ssg[tmpeta1][tmpphi2].push_back(hocro);
1314  crossg[tmpeta1][tmpphi2]->Fill(hocro);
1315  }
1316  }
1317 
1318  if (tmpphi2 >=0 && tmpphi2 <nphimx) {
1319  crossg[tmpeta1][nphimx]->Fill(hocro);
1320  }
1321  if (m_combined) {
1322  com_crossg[iring2][tmprout2-1]->Fill(hocro);
1323  com_crossg[iring2][tmproutmx]->Fill(hocro);
1324  }
1325  }
1326 
1327  if (iselect==1) {
1328  for (int ij=0; ij<neffip; ij++) {
1329  if (ij==0) {
1330  sig_effi[ij]->Fill(ieta, iphi, 1.);
1331  } else {
1332  if (hosig[4] >ij*m_sigma) {
1333  sig_effi[ij]->Fill(ieta, iphi, 1.);
1334  }
1335  }
1336  }
1337 
1338  tmpphi1 = iphi - 1;
1339 
1340  if (sig_reg[tmpeta1][tmpphi1].size()<4000 ) {
1341  if (hosig[4]>-50&& hosig[4] <15) {
1342  sigrsg[tmpeta1][tmpphi1]->Fill(hosig[4]);
1343  if (!m_histfit && hosig[4]<=ahigh/2.) sig_reg[tmpeta1][tmpphi1].push_back(hosig[4]);
1344  invang[tmpeta1][tmpphi1] += 1./fabs(hoang);
1345  }
1346  }
1347 
1348  if (tmpphi1 >=0 && tmpphi1 <nphimx) { //GREN
1349  sigrsg[tmpeta1][nphimx]->Fill(hosig[4]);
1350  invang[tmpeta1][nphimx] += 1./fabs(hoang);
1351  }
1352 
1353  if (m_combined) { //#ifdef COMBINED
1354  com_sigrsg[iring2][tmprout-1]->Fill(hosig[4]);
1355  com_invang[iring2][tmprout-1] += 1./fabs(hoang);
1356 
1357  com_sigrsg[iring2][tmproutmx]->Fill(hosig[4]);
1358  com_invang[iring2][tmproutmx] += 1./fabs(hoang);
1359  } //m_combined #endif
1360 
1361  if (m_checkmap || m_correl) { //#ifdef CHECKMAP
1362  tmpeta = etamap[itag][npixel];
1363  tmpphi = phimap[itag][npixel];
1364  if (tmpeta>=0 && tmpphi >=0) {
1365  if (iflip !=0) tmpphi = abs(tmpphi-2);
1366  if (int((hocorsig[fact*tmpeta+tmpphi]-hosig[4])*10000)/10000.!=0) {
1367  iaxxx++;
1368  cout<<"iring2xxx "<<irun<<" "<<ievt<<" "<<isect<<" "<<iring<<" "<<tmpsect<<" "<<ieta<<" "<<iphi<<" "<<npixel<<" "<<tmpeta<<" "<<tmpphi<<" "<<tmpeta1<<" "<<tmpphi1<<" itag "<<itag<<" "<<iflip<<" "<<fact<<" "<<hocorsig[fact*tmpeta+tmpphi]<<" "<<fact*tmpeta+tmpphi<<" "<<hosig[4]<<" "<<hodx<<" "<<hody<<endl;
1369 
1370  for (int i=0; i<18; i++) {cout <<" "<<i<<" "<<hocorsig[i];}
1371  cout<<" ix "<<iaxxx<<" "<<ibxxx<<endl;
1372  } else { ibxxx++; }
1373 
1374  corrsgc[tmpeta1][tmpphi1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
1375  if (m_combined) { //#ifdef COMBINED
1376  com_corrsgc[iring2][tmpsect-1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
1377  } //m_combined #endif
1378  }
1379  } //m_checkmap #endif
1380 
1381  if (m_correl) { //#ifdef CORREL
1382  float allcorsig = 0.0;
1383 
1384  tmpeta = etamap[itag][npixleft[npixel]];
1385  tmpphi = phimap[itag][npixleft[npixel]];
1386 
1387  if (tmpeta>=0 && tmpphi >=0) {
1388  if (iflip !=0) tmpphi = abs(tmpphi-2);
1389  corrsgl[tmpeta1][tmpphi1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
1390  allcorsig +=hocorsig[fact*tmpeta+tmpphi];
1391  if (m_combined) { //#ifdef COMBINED
1392  com_corrsgl[iring2][tmpsect-1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
1393  } //m_combined #endif
1394  }
1395 
1396  tmpeta = etamap[itag][npixrigh[npixel]];
1397  tmpphi = phimap[itag][npixrigh[npixel]];
1398  if (tmpeta>=0 && tmpphi >=0) {
1399  if (iflip !=0) tmpphi = abs(tmpphi-2);
1400  corrsgr[tmpeta1][tmpphi1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
1401  allcorsig +=hocorsig[fact*tmpeta+tmpphi];
1402  if (m_combined) { // #ifdef COMBINED
1403  com_corrsgr[iring2][tmpsect-1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
1404  } //m_combined #endif
1405  }
1406 
1407  tmpeta = etamap[itag][npixlebt[npixel]];
1408  tmpphi = phimap[itag][npixlebt[npixel]];
1409  if (tmpeta>=0 && tmpphi >=0) {
1410  if (iflip !=0) tmpphi = abs(tmpphi-2);
1411  corrsglb[tmpeta1][tmpphi1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
1412  allcorsig +=hocorsig[fact*tmpeta+tmpphi];
1413  if (m_combined){ //#ifdef COMBINED
1414  com_corrsglb[iring2][tmpsect-1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
1415  } //m_combined #endif
1416  }
1417 
1418  tmpeta = etamap[itag][npixribt[npixel]];
1419  tmpphi = phimap[itag][npixribt[npixel]];
1420  if (tmpeta>=0 && tmpphi >=0) {
1421  if (iflip !=0) tmpphi = abs(tmpphi-2);
1422  corrsgrb[tmpeta1][tmpphi1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
1423  allcorsig +=hocorsig[fact*tmpeta+tmpphi];
1424  if (m_combined) { // #ifdef COMBINED
1425  com_corrsgrb[iring2][tmpsect-1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
1426  } //m_combined #endif
1427  }
1428 
1429  tmpeta = etamap[itag][npixleup[npixel]];
1430  tmpphi = phimap[itag][npixleup[npixel]];
1431  if (tmpeta>=0 && tmpphi >=0) {
1432  if (iflip !=0) tmpphi = abs(tmpphi-2);
1433  corrsglu[tmpeta1][tmpphi1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
1434  allcorsig +=hocorsig[fact*tmpeta+tmpphi];
1435  if (m_combined) {// #ifdef COMBINED
1436  com_corrsglu[iring2][tmpsect-1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
1437  } //m_combined #endif
1438  }
1439 
1440  tmpeta = etamap[itag][npixriup[npixel]];
1441  tmpphi = phimap[itag][npixriup[npixel]];
1442  if (tmpeta>=0 && tmpphi >=0) {
1443  if (iflip !=0) tmpphi = abs(tmpphi-2);
1444  corrsgru[tmpeta1][tmpphi1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
1445  allcorsig +=hocorsig[fact*tmpeta+tmpphi];
1446  if (m_combined) { // #ifdef COMBINED
1447  com_corrsgru[iring2][tmpsect-1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
1448  } //m_combined #endif
1449  }
1450  corrsgall[tmpeta1][tmpphi1]->Fill(allcorsig);
1451  if (m_combined) { // #ifdef COMBINED
1452  com_corrsgall[iring2][tmpsect-1]->Fill(allcorsig);
1453  } //m_combined #endif
1454 
1455 
1456  } //m_correl #endif
1457  for (int k=0; k<9; k++) {
1458  switch (iring)
1459  {
1460  case 2 : ho_sig2p[k]->Fill(hosig[k]); break;
1461  case 1 : ho_sig1p[k]->Fill(hosig[k]); break;
1462  case 0 : ho_sig00[k]->Fill(hosig[k]); break;
1463  case -1 : ho_sig1m[k]->Fill(hosig[k]); break;
1464  case -2 : ho_sig2m[k]->Fill(hosig[k]); break;
1465  }
1466  if (m_hbinfo) { // #ifdef HBINFO
1467  hbhe_sig[k]->Fill(hbhesig[k]);
1468  // cout <<"hbhe "<<k<<" "<<hbhesig[k]<<endl;
1469  } //m_hbinfo #endif
1470  }
1471  } //if (iselect==1)
1472  } //for (HOCalibVariableCollection::const_iterator hoC=(*HOCalib).begin(); hoC!=(*HOCalib).end(); hoC++){
1473 
1474  } //if (isCosMu)
1475 
1476 #ifdef THIS_IS_AN_EVENT_EXAMPLE
1477  Handle<ExampleData> pIn;
1478  iEvent.getByLabel("example",pIn);
1479 #endif
1480 
1481 #ifdef THIS_IS_AN_EVENTSETUP_EXAMPLE
1482  ESHandle<SetupData> pSetup;
1483  iSetup.get<SetupRecord>().get(pSetup);
1484 #endif
1485 }
1486 
1487 
1488 // ------------ method called once each job just before starting event loop ------------
1489 void
1490 HOCalibAnalyzer::beginJob()
1491 {
1492 }
1493 
1494 // ------------ method called once each job just after ending the event loop ------------
1495 void
1496 HOCalibAnalyzer::endJob() {
1497 
1498  theFile->cd();
1499 
1500  for (int i=0; i<nphimx; i++) {
1501  for (int j=0; j<netamx; j++) {
1502 
1503  nevsigrsg->Fill(netamx*i+j, sigrsg[j][i]->GetEntries());
1504  mnsigrsg->Fill(netamx*i+j, sigrsg[j][i]->GetMean());
1505  rmssigrsg->Fill(netamx*i+j, sigrsg[j][i]->GetRMS());
1506 
1507  nevcrossg->Fill(netamx*i+j, crossg[j][i]->GetEntries());
1508  mncrossg->Fill(netamx*i+j, crossg[j][i]->GetMean());
1509  rmscrossg->Fill(netamx*i+j, crossg[j][i]->GetRMS());
1510 
1511  if (m_correl) { //#ifdef CORREL
1512 
1513  nevcorrsglb->Fill(netamx*i+j, corrsglb[j][i]->GetEntries());
1514  mncorrsglb->Fill(netamx*i+j, corrsglb[j][i]->GetMean());
1515  rmscorrsglb->Fill(netamx*i+j, corrsglb[j][i]->GetRMS());
1516 
1517  nevcorrsgrb->Fill(netamx*i+j, corrsgrb[j][i]->GetEntries());
1518  mncorrsgrb->Fill(netamx*i+j, corrsgrb[j][i]->GetMean());
1519  rmscorrsgrb->Fill(netamx*i+j, corrsgrb[j][i]->GetRMS());
1520 
1521  nevcorrsglu->Fill(netamx*i+j, corrsglu[j][i]->GetEntries());
1522  mncorrsglu->Fill(netamx*i+j, corrsglu[j][i]->GetMean());
1523  rmscorrsglu->Fill(netamx*i+j, corrsglu[j][i]->GetRMS());
1524 
1525  nevcorrsgru->Fill(netamx*i+j, corrsgru[j][i]->GetEntries());
1526  mncorrsgru->Fill(netamx*i+j, corrsgru[j][i]->GetMean());
1527  rmscorrsgru->Fill(netamx*i+j, corrsgru[j][i]->GetRMS());
1528 
1529  nevcorrsgall->Fill(netamx*i+j, corrsgall[j][i]->GetEntries());
1530  mncorrsgall->Fill(netamx*i+j, corrsgall[j][i]->GetMean());
1531  rmscorrsgall->Fill(netamx*i+j, corrsgall[j][i]->GetRMS());
1532 
1533  nevcorrsgl->Fill(netamx*i+j, corrsgl[j][i]->GetEntries());
1534  mncorrsgl->Fill(netamx*i+j, corrsgl[j][i]->GetMean());
1535  rmscorrsgl->Fill(netamx*i+j, corrsgl[j][i]->GetRMS());
1536 
1537  nevcorrsgr->Fill(netamx*i+j, corrsgr[j][i]->GetEntries());
1538  mncorrsgr->Fill(netamx*i+j, corrsgr[j][i]->GetMean());
1539  rmscorrsgr->Fill(netamx*i+j, corrsgr[j][i]->GetRMS());
1540  } //m_correl #endif
1541  if (m_checkmap) { //#ifdef CHECKMAP
1542  nevcorrsgc->Fill(netamx*i+j, corrsgc[j][i]->GetEntries());
1543  mncorrsgc->Fill(netamx*i+j, corrsgc[j][i]->GetMean());
1544  rmscorrsgc->Fill(netamx*i+j, corrsgc[j][i]->GetRMS());
1545  } //m_checkmap #endif
1546  }
1547  }
1548 
1549  if (m_combined) { // #ifdef COMBINED
1550  for (int j=0; j<ringmx; j++) {
1551  for (int i=0; i<routmx; i++) {
1552  if (j!=2 && i>=rout12mx) continue;
1553  nevsigrsg->Fill(netamx*nphimx+ringmx*i+j, com_sigrsg[j][i]->GetEntries());
1554  mnsigrsg->Fill(netamx*nphimx+ringmx*i+j, com_sigrsg[j][i]->GetMean());
1555  rmssigrsg->Fill(netamx*nphimx+ringmx*i+j, com_sigrsg[j][i]->GetRMS());
1556 
1557  nevcrossg->Fill(netamx*nphimx+ringmx*i+j, com_crossg[j][i]->GetEntries());
1558  mncrossg->Fill(netamx*nphimx+ringmx*i+j, com_crossg[j][i]->GetMean());
1559  rmscrossg->Fill(netamx*nphimx+ringmx*i+j, com_crossg[j][i]->GetRMS());
1560  }
1561  }
1562 
1563  for (int i=0; i<sectmx; i++) {
1564  for (int j=0; j<ringmx; j++) {
1565  if (m_correl) { // #ifdef CORREL
1566  nevcorrsglb->Fill(netamx*nphimx+ringmx*i+j, com_corrsglb[j][i]->GetEntries());
1567  mncorrsglb->Fill(netamx*nphimx+ringmx*i+j, com_corrsglb[j][i]->GetMean());
1568  rmscorrsglb->Fill(netamx*nphimx+ringmx*i+j, com_corrsglb[j][i]->GetRMS());
1569 
1570  nevcorrsgrb->Fill(netamx*nphimx+ringmx*i+j, com_corrsgrb[j][i]->GetEntries());
1571  mncorrsgrb->Fill(netamx*nphimx+ringmx*i+j, com_corrsgrb[j][i]->GetMean());
1572  rmscorrsgrb->Fill(netamx*nphimx+ringmx*i+j, com_corrsgrb[j][i]->GetRMS());
1573 
1574  nevcorrsglu->Fill(netamx*nphimx+ringmx*i+j, com_corrsglu[j][i]->GetEntries());
1575  mncorrsglu->Fill(netamx*nphimx+ringmx*i+j, com_corrsglu[j][i]->GetMean());
1576  rmscorrsglu->Fill(netamx*nphimx+ringmx*i+j, com_corrsglu[j][i]->GetRMS());
1577 
1578  nevcorrsgru->Fill(netamx*nphimx+ringmx*i+j, com_corrsgru[j][i]->GetEntries());
1579  mncorrsgru->Fill(netamx*nphimx+ringmx*i+j, com_corrsgru[j][i]->GetMean());
1580  rmscorrsgru->Fill(netamx*nphimx+ringmx*i+j, com_corrsgru[j][i]->GetRMS());
1581 
1582  nevcorrsgall->Fill(netamx*nphimx+ringmx*i+j, com_corrsgall[j][i]->GetEntries());
1583  mncorrsgall->Fill(netamx*nphimx+ringmx*i+j, com_corrsgall[j][i]->GetMean());
1584  rmscorrsgall->Fill(netamx*nphimx+ringmx*i+j, com_corrsgall[j][i]->GetRMS());
1585 
1586  nevcorrsgl->Fill(netamx*nphimx+ringmx*i+j, com_corrsgl[j][i]->GetEntries());
1587  mncorrsgl->Fill(netamx*nphimx+ringmx*i+j, com_corrsgl[j][i]->GetMean());
1588  rmscorrsgl->Fill(netamx*nphimx+ringmx*i+j, com_corrsgl[j][i]->GetRMS());
1589 
1590  nevcorrsgr->Fill(netamx*nphimx+ringmx*i+j, com_corrsgr[j][i]->GetEntries());
1591  mncorrsgr->Fill(netamx*nphimx+ringmx*i+j, com_corrsgr[j][i]->GetMean());
1592  rmscorrsgr->Fill(netamx*nphimx+ringmx*i+j, com_corrsgr[j][i]->GetRMS());
1593  } //m_correl #endif
1594  if (m_checkmap) { // #ifdef CHECKMAP
1595  nevcorrsgc->Fill(netamx*nphimx+ringmx*i+j, com_corrsgc[j][i]->GetEntries());
1596  mncorrsgc->Fill(netamx*nphimx+ringmx*i+j, com_corrsgc[j][i]->GetMean());
1597  rmscorrsgc->Fill(netamx*nphimx+ringmx*i+j, com_corrsgc[j][i]->GetRMS());
1598  } //m_checkmap #endif
1599  }
1600  }
1601  } //m_combined #endif
1602 
1603 
1604  for (int i=1; i<neffip; i++) {
1605  sig_effi[i]->Divide(sig_effi[0]);
1606  }
1607  for (int ij=0; ij<netamx; ij++) {
1608  for (int jk = 0; jk <nphimx; jk++) {
1609  int ieta = (ij<15) ? ij+1 : 14-ij;
1610  int iphi = jk+1;
1611  double signal = sigrsg[ij][jk]->GetMean();
1612  mean_energy->Fill(ieta, iphi, signal);
1613  }
1614  }
1615 
1616  int irunold = irun;
1617 
1618 
1619  gStyle->SetOptLogy(0);
1620  gStyle->SetTitleFillColor(10);
1621  gStyle->SetStatColor(10);
1622 
1623  gStyle->SetCanvasColor(10);
1624  gStyle->SetOptStat(0); //1110);
1625  gStyle->SetOptTitle(1);
1626 
1627  gStyle->SetTitleColor(10);
1628  gStyle->SetTitleFontSize(0.09);
1629  gStyle->SetTitleOffset(-0.05);
1630  gStyle->SetTitleBorderSize(1);
1631 
1632  gStyle->SetPadColor(10);
1633  gStyle->SetPadBorderMode(0);
1634  gStyle->SetStatColor(10);
1635  gStyle->SetPadBorderMode(0);
1636  gStyle->SetStatBorderSize(1);
1637  gStyle->SetStatFontSize(.07);
1638 
1639  gStyle->SetStatStyle(1001);
1640  gStyle->SetOptFit(101);
1641  gStyle->SetCanvasColor(10);
1642  gStyle->SetCanvasBorderMode(0);
1643 
1644  gStyle->SetStatX(.99);
1645  gStyle->SetStatY(.99);
1646  gStyle->SetStatW(.45);
1647  gStyle->SetStatH(.16);
1648  gStyle->SetLabelSize(0.075,"XY");
1649  gStyle->SetLabelOffset(0.21,"XYZ");
1650  gStyle->SetTitleSize(0.065,"XY");
1651  gStyle->SetTitleOffset(0.06,"XYZ");
1652  gStyle->SetPadTopMargin(.09);
1653  gStyle->SetPadBottomMargin(0.11);
1654  gStyle->SetPadLeftMargin(0.12);
1655  gStyle->SetPadRightMargin(0.15);
1656  gStyle->SetPadGridX(3);
1657  gStyle->SetPadGridY(3);
1658  gStyle->SetGridStyle(2);
1659  gStyle->SetNdivisions(303,"XY");
1660 
1661  gStyle->SetMarkerSize(0.60);
1662  gStyle->SetMarkerColor(2);
1663  gStyle->SetMarkerStyle(20);
1664  gStyle->SetTitleFontSize(0.07);
1665 
1666  char out_file[200];
1667  int xsiz = 700;
1668  int ysiz = 500;
1669 
1670  TCanvas *c2 = new TCanvas("c2", "Statistics and efficiency", xsiz, ysiz);
1671  c2->Divide(2,1); //(3,2);
1672  for (int ij=0; ij<neffip; ij=ij+3) {
1673  sig_effi[ij]->GetXaxis()->SetTitle("#eta");
1674  sig_effi[ij]->GetXaxis()->SetTitleSize(0.075);
1675  sig_effi[ij]->GetXaxis()->SetTitleOffset(0.65); //0.85
1676  sig_effi[ij]->GetXaxis()->CenterTitle();
1677  sig_effi[ij]->GetXaxis()->SetLabelSize(0.055);
1678  sig_effi[ij]->GetXaxis()->SetLabelOffset(0.001);
1679 
1680  sig_effi[ij]->GetYaxis()->SetTitle("#phi");
1681  sig_effi[ij]->GetYaxis()->SetTitleSize(0.075);
1682  sig_effi[ij]->GetYaxis()->SetTitleOffset(0.9);
1683  sig_effi[ij]->GetYaxis()->CenterTitle();
1684  sig_effi[ij]->GetYaxis()->SetLabelSize(0.055);
1685  sig_effi[ij]->GetYaxis()->SetLabelOffset(0.01);
1686 
1687  c2->cd(int(ij/3.)+1); sig_effi[ij]->Draw("colz");
1688  }
1689  sprintf(out_file, "comb_hosig_evt_%i.jpg",irunold);
1690  c2->SaveAs(out_file);
1691 
1692  gStyle->SetTitleFontSize(0.045);
1693  gStyle->SetPadRightMargin(0.1);
1694  gStyle->SetPadLeftMargin(0.1);
1695  gStyle->SetPadBottomMargin(0.12);
1696 
1697  TCanvas *c1 = new TCanvas("c1", "Mean signal in each tower", xsiz, ysiz);
1698 
1699  mean_energy->GetXaxis()->SetTitle("#eta");
1700  mean_energy->GetXaxis()->SetTitleSize(0.075);
1701  mean_energy->GetXaxis()->SetTitleOffset(0.65); //0.6
1702  mean_energy->GetXaxis()->CenterTitle();
1703  mean_energy->GetXaxis()->SetLabelSize(0.045);
1704  mean_energy->GetXaxis()->SetLabelOffset(0.001);
1705 
1706  mean_energy->GetYaxis()->SetTitle("#phi");
1707  mean_energy->GetYaxis()->SetTitleSize(0.075);
1708  mean_energy->GetYaxis()->SetTitleOffset(0.5);
1709  mean_energy->GetYaxis()->CenterTitle();
1710  mean_energy->GetYaxis()->SetLabelSize(0.045);
1711  mean_energy->GetYaxis()->SetLabelOffset(0.01);
1712 
1713  mean_energy->Draw("colz");
1714  sprintf(out_file, "homean_energy_%i.jpg",irunold);
1715  c1->SaveAs(out_file);
1716 
1717  delete c1;
1718  delete c2;
1719 
1720  gStyle->SetPadBottomMargin(0.14);
1721  gStyle->SetPadLeftMargin(0.17);
1722  gStyle->SetPadRightMargin(0.03);
1723 
1724  gStyle->SetOptStat(1110);
1725 
1726  const int nsample =8;
1727  TF1* gx0[nsample]={0};
1728  TF1* ped0fun[nsample]={0};
1729  TF1* signal[nsample]={0};
1730  TF1* pedfun[nsample]={0};
1731  TF1* sigfun[nsample]={0};
1732  TF1* signalx[nsample]={0};
1733 
1734  TH1F* signall[nsample]={0};
1735  TH1F* pedstll[nsample]={0};
1736 
1737  if (m_constant) {
1738 
1739  gStyle->SetOptFit(101);
1740  gStyle->SetCanvasBorderMode(0);
1741  gStyle->SetPadBorderMode(0);
1742  gStyle->SetStatBorderSize(1);
1743  gStyle->SetStatStyle(1001);
1744  gStyle->SetTitleColor(10);
1745  gStyle->SetTitleFontSize(0.09);
1746  gStyle->SetTitleOffset(-0.05);
1747  gStyle->SetTitleBorderSize(1);
1748 
1749  gStyle->SetCanvasColor(10);
1750  gStyle->SetPadColor(10);
1751  gStyle->SetStatColor(10);
1752  gStyle->SetStatFontSize(.07);
1753  gStyle->SetStatX(0.99);
1754  gStyle->SetStatY(0.99);
1755  gStyle->SetStatW(0.30);
1756  gStyle->SetStatH(0.10);
1757  gStyle->SetTitleSize(0.065,"XYZ");
1758  gStyle->SetLabelSize(0.075,"XYZ");
1759  gStyle->SetLabelOffset(0.012,"XYZ");
1760  gStyle->SetPadGridX(1);
1761  gStyle->SetPadGridY(1);
1762  gStyle->SetGridStyle(3);
1763  gStyle->SetNdivisions(101,"XY");
1764  gStyle->SetOptLogy(0);
1765  int iiter = 0;
1766 
1767 
1768  std::ofstream file_out(theoutputtxtFile.c_str());
1769  // TPostScript* ps=0;
1770  int ips=111;
1771  TPostScript ps(theoutputpsFile.c_str(),ips);
1772  ps.Range(20,28);
1773 
1774  xsiz = 900; //900;
1775  ysiz = 1200; //600;
1776  TCanvas *c0 = new TCanvas("c0", " Pedestal vs signal", xsiz, ysiz);
1777 
1778 // Fix is done for eta-phi
1779 
1780  float mean_eta[netamx];
1781  float mean_phi[nphimx];
1782  float rms_eta[netamx];
1783  float rms_phi[nphimx];
1784 
1785  for (int ij=0; ij<nphimx; ++ij) {mean_phi[ij] = 0; rms_phi[ij] = 0;}
1786  for (int ij=0; ij<netamx; ++ij) {mean_eta[ij] = 0; rms_eta[ij] = 0;}
1787 
1788  int mxeta = 0;
1789  int mxphi = 0;
1790  int mneta = 0;
1791  int mnphi = 0;
1792 
1793  //iijj = 0 : Merging all ring
1794  // = 1 : Individual HPD
1795  //iijj = 2 : merging all phi
1796  // = 3 : Individual tower
1797 
1798  for (int iijj = 0; iijj <4; iijj++) {
1799  // if ((!mx_combined) && iijj==1) continue; //Use this only for combined data
1800  if (iijj==0){
1801  mxeta = ringmx; mxphi = 1; mneta = 0; mnphi = 0;
1802  } else if (iijj==1) {
1803  mxeta = ringmx; mxphi = routmx;
1804  mneta = 0; mnphi = 0;
1805  } else if (iijj==2) {
1806  mxeta = netamx; mxphi = 1; mneta = 0; mnphi = 0;
1807  } else if (iijj==3) {
1808  mxeta = netamx; mxphi = nphimx;
1809  mneta = 0; mnphi = 0;
1810  }
1811 
1812  for (int jk=mneta; jk<mxeta; jk++) {
1813  for (int ij=mnphi; ij<mxphi; ij++) {
1814  if (iijj==1) continue;
1815  if ((iijj==0 || iijj==1) && jk !=2 && ij >=rout12mx) continue;
1816  int izone = iiter%nsample;
1817 
1818  if (iijj==0) {
1819  int iread = (jk==2) ? routmx : rout12mx;
1820  signall[izone] = (TH1F*)com_sigrsg[jk][iread]->Clone("hnew");
1821  pedstll[izone] = (TH1F*)com_crossg[jk][iread]->Clone("hnew");
1822  } else if (iijj==1) {
1823  signall[izone] = (TH1F*)com_sigrsg[jk][ij]->Clone("hnew");
1824  pedstll[izone] = (TH1F*)com_crossg[jk][ij]->Clone("hnew");
1825  } else if (iijj==2) {
1826  signall[izone] = (TH1F*)sigrsg[jk][nphimx]->Clone("hnew");
1827  pedstll[izone] = (TH1F*)crossg[jk][nphimx]->Clone("hnew");
1828  } else if (iijj==3) {
1829  signall[izone] = (TH1F*)sigrsg[jk][ij]->Clone("hnew");
1830  pedstll[izone] = (TH1F*)crossg[jk][ij]->Clone("hnew");
1831  }
1832 
1833  pedstll[izone]->SetLineWidth(2);
1834  signall[izone]->SetLineWidth(2);
1835  pedstll[izone]->SetLineColor(2);
1836  signall[izone]->SetLineColor(4);
1837  pedstll[izone]->SetNdivisions(506,"XY");
1838  signall[izone]->SetNdivisions(506,"XY");
1839 
1840  signall[izone]->GetXaxis()->SetLabelSize(.065);
1841  signall[izone]->GetYaxis()->SetLabelSize(.06);
1842  if (m_digiInput) {
1843  signall[izone]->GetXaxis()->SetTitle("Signal (fC)");
1844  } else {
1845  signall[izone]->GetXaxis()->SetTitle("Signal (GeV)");
1846  }
1847  signall[izone]->GetXaxis()->SetTitleSize(.065);
1848  signall[izone]->GetXaxis()->CenterTitle();
1849 
1850  if (izone==0) { //iiter%8 ==0) {
1851  ps.NewPage();
1852  c0->Divide(4,4); //c0->Divide(2,4); // c0->Divide(1,2);
1853  }
1854  c0->cd(2*izone+1); // (iiter%8)+1); //c0->cd(iiter%8+1);
1855 
1856  /*
1857  if (iijj==0 && izone==0) {
1858  gStyle->SetOptLogy(1);
1859  gStyle->SetOptStat(0);
1860  gStyle->SetOptFit(0);
1861  c0->Divide(3,2);
1862  }
1863 
1864  if (iijj>0) {
1865  gStyle->SetOptLogy(0);
1866  gStyle->SetOptStat(1110);
1867  gStyle->SetOptFit(101);
1868 
1869  if (iiter==0) {
1870  int ips=111;
1871  ps = new TPostScript(theoutputpsFile.c_str(),ips);
1872  ps.Range(20,28);
1873  xsiz = 900; //900;
1874  ysiz = 1200; //600;
1875  c0 = new TCanvas("c0", " Pedestal vs signal", xsiz, ysiz);
1876  }
1877  if (izone==0) {
1878  ps.NewPage();
1879  c0->Divide(4,4);
1880  }
1881  }
1882  if (iijj==0) {c0->cd(izone+1); } else { c0->cd(2*izone+1);}
1883  */
1884 
1885  float mean = pedstll[izone]->GetMean();
1886  float rms = pedstll[izone]->GetRMS();
1887  if (m_digiInput) {
1888  if (rms <0.6) rms = 0.6;
1889  if (rms >1.2) rms = 1.2;
1890  if (mean >1.2) mean = 1.2;
1891  if (mean <-1.2) mean = -1.2;
1892  } else {
1893  if (rms <0.10) rms = 0.10;
1894  if (rms >0.15) rms=0.15;
1895  if (mean >0.20) mean = 0.20;
1896  if (mean <-0.20) mean = -0.20;
1897  }
1898  float xmn = mean-6.*rms;
1899  float xmx = mean+6.*rms;
1900 
1901  binwid = pedstll[izone]->GetBinWidth(1);
1902  if (xmx > pedstll[izone]->GetXaxis()->GetXmax()) xmx = pedstll[izone]->GetXaxis()->GetXmax()-0.5*binwid;
1903  if (xmn < pedstll[izone]->GetXaxis()->GetXmin()) xmn = pedstll[izone]->GetXaxis()->GetXmin()+0.5*binwid;
1904 
1905  float height = pedstll[izone]->GetEntries();
1906 
1907  double par[nbgpr] ={height, mean, 0.75*rms};
1908 
1909  double gaupr[nbgpr];
1910  double parer[nbgpr];
1911 
1912  ietafit = jk;
1913  iphifit = ij;
1914  pedstll[izone]->GetXaxis()->SetLabelSize(.065);
1915  pedstll[izone]->GetYaxis()->SetLabelSize(.06);
1916 
1917  // if (iijj==0) {
1918  // pedstll[izone]->GetXaxis()->SetRangeUser(alow, ahigh);
1919  // } else {
1920  pedstll[izone]->GetXaxis()->SetRangeUser(xmn, xmx);
1921  // }
1922 
1923  if (m_digiInput) {
1924  if (iijj==0) {
1925  pedstll[izone]->GetXaxis()->SetTitle("Pedestal/Signal (fC)");
1926  } else {
1927  pedstll[izone]->GetXaxis()->SetTitle("Pedestal (fC)");
1928  }
1929  } else {
1930  if (iijj==0) {
1931  pedstll[izone]->GetXaxis()->SetTitle("Pedestal/Signal (GeV)");
1932  } else {
1933  pedstll[izone]->GetXaxis()->SetTitle("Pedestal (GeV)");
1934  }
1935  }
1936  pedstll[izone]->GetXaxis()->SetTitleSize(.065);
1937  pedstll[izone]->GetXaxis()->CenterTitle();
1938  // pedstll[izone]->SetLineWidth(2);
1939 
1940  pedstll[izone]->Draw();
1941  if (m_pedsuppr && !m_digiInput) {
1942  gaupr[0] = 0;
1943  gaupr[1] = 0.0; // pedmean[ietafit][iphifit];
1944  if (m_digiInput) {
1945  gaupr[2] = 0.90; //GMA need from database pedwidth[ietafit][iphifit];
1946  } else {
1947  gaupr[2] = 0.15; //GMA need from database
1948  }
1949  parer[0] = parer[1] = parer[2] = 0;
1950  } else {
1951 
1952  if (pedstll[izone]->GetEntries() >5) {
1953 
1954  if ((iijj!=3) || m_histfit) {
1955  char temp[20];
1956  sprintf(temp, "gx0_%i",izone);
1957  gx0[izone] = new TF1(temp, gausX, xmn, xmx, nbgpr);
1958  gx0[izone]->SetParameters(par);
1959  gx0[izone]->SetLineWidth(1);
1960  pedstll[izone]->Fit(gx0[izone], "R+");
1961 
1962  for (int k=0; k<nbgpr; k++) {
1963  parer[k] = gx0[izone]->GetParError(k);
1964  gaupr[k] = gx0[izone]->GetParameter(k);
1965  }
1966  } else {
1967  double strt[nbgpr] = {height, mean, 0.75*rms};
1968  double step[nbgpr] = {1.0, 0.001, 0.001};
1969  double alowmn[nbgpr] = {0.5*height, mean-rms, 0.3*rms};
1970  double ahighmn[nbgpr] ={1.5*height, mean+rms, 1.5*rms};
1971 
1972  TMinuit *gMinuit = new TMinuit(nbgpr);
1973  gMinuit->SetFCN(fcnbg);
1974 
1975  double arglist[10];
1976  int ierflg = 0;
1977  arglist[0] =0.5;
1978  gMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
1979  char name[100];
1980  for (int k=0; k<nbgpr; k++) {
1981  sprintf(name, "pedpar%i",k);
1982  gMinuit->mnparm(k, name, strt[k], step[k], alowmn[k], ahighmn[k],ierflg);
1983  }
1984 
1985  arglist[0] = 0;
1986  gMinuit->mnexcm("SIMPLEX", arglist, 0, ierflg);
1987 
1988  arglist[0] = 0;
1989  gMinuit->mnexcm("IMPROVE", arglist, 0, ierflg);
1990 
1991  TString chnam;
1992  double parv,err,xlo,xup, plerr, mierr, eparab, gcc;
1993  int iuit;
1994 
1995  for (int k=0; k<nbgpr; k++) {
1996  if (step[k] >-10) {
1997  gMinuit->mnpout(k, chnam, parv, err, xlo, xup, iuit);
1998  gMinuit->mnerrs(k, plerr, mierr, eparab, gcc);
1999  // cout <<"k "<< k<<" "<<chnam<<" "<<parv<<" "<<err<<" "<<xlo<<" "<<xup<<" "<<plerr<<" "<<mierr<<" "<<eparab<<endl;
2000  if (k==0) {
2001  gaupr[k] = parv*binwid;
2002  parer[k] = err*binwid;
2003  } else {
2004  gaupr[k] = parv;
2005  parer[k] = err;
2006  }
2007  }
2008  }
2009 
2010  // gx0[izone]->SetParameters(gaupr);
2011 
2012  char temp[20];
2013  sprintf(temp, "ped0fun_%i",izone);
2014  ped0fun[izone] = new TF1(temp, gausX, xmn, xmx, nbgpr);
2015  ped0fun[izone]->SetParameters(gaupr);
2016  ped0fun[izone]->SetLineColor(3);
2017  ped0fun[izone]->SetLineWidth(1);
2018  ped0fun[izone]->Draw("same");
2019 
2020  delete gMinuit;
2021  }
2022  } else {
2023  for (int k=0; k<nbgpr; k++) {gaupr[k] = par[k]; }
2024  if (m_digiInput) { gaupr[2] = 0.90; } else { gaupr[2] = 0.15;}
2025  }
2026  }
2027  // if (iijj!=0)
2028  c0->cd(2*izone+2);
2029  if (signall[izone]->GetEntries() >5) {
2030  Double_t parall[nsgpr];
2031  double parserr[nsgpr];
2032  double fitres[nsgpr];
2033  double pedht = 0;
2034 
2035  char temp[20];
2036  sprintf(temp, "signal_%i",izone);
2037  xmn = signall[izone]->GetXaxis()->GetXmin();
2038  xmx = 0.5*signall[izone]->GetXaxis()->GetXmax();
2039  signal[izone] = new TF1(temp, totalfunc, xmn, xmx, nsgpr);
2040  xmx *=2.0;
2041  if ((iijj!=3) || m_histfit) {
2042  pedht = (signall[izone]->GetBinContent(nbn-1)+
2043  signall[izone]->GetBinContent(nbn)+
2044  signall[izone]->GetBinContent(nbn+1))/3.;
2045 
2046  if (m_pedsuppr && !m_digiInput) {
2047  parall[1] = 0.0; // pedmean[ietafit][iphifit];
2048  if (m_digiInput) { parall[2] = 0.90; } else { parall[2] = 0.15;}
2049  } else {
2050  for (int i=0; i<nbgpr; i++) {parall[i] = gaupr[i];}
2051  }
2052 
2053  set_mean(parall[1], m_digiInput);
2054  set_sigma(parall[2], m_digiInput);
2055 
2056  parall[0] = 0.9*pedht; //GM for Z-mumu, there is almost no pedestal
2057  parall[3] = 0.14;
2058  double area = binwid*signall[izone]->GetEntries();
2059  parall[5]= area;
2060 
2061  if (iijj==3) {
2062  parall[4] = fitprm[4][jk];
2063  parall[6] = fitprm[6][jk];
2064  } else {
2065  parall[4] = signall[izone]->GetMean();
2066  parall[6]=parall[2];
2067  }
2068 
2069  signal[izone]->SetParameters(parall);
2070  signal[izone]->FixParameter(1, parall[1]);
2071  signal[izone]->FixParameter(2, parall[2]);
2072  signal[izone]->SetParLimits(0, 0.00, 2.0*pedht+0.1);
2073  signal[izone]->FixParameter(3, 0.14);
2074 
2075  signal[izone]->SetParLimits(5, 0.40*area, 1.15*area);
2076  // if (m_histfit) { //GMA
2077  if (iijj==3) {
2078  signal[izone]->SetParLimits(4, 0.2*fitprm[4][jk], 2.0*fitprm[4][jk]);
2079  signal[izone]->SetParLimits(6, 0.2*fitprm[6][jk], 2.0*fitprm[6][jk]);
2080  } else {
2081  if (m_digiInput) {
2082  signal[izone]->SetParLimits(4, 0.6, 6.0);
2083  signal[izone]->SetParLimits(6, 0.60, 3.0);
2084  } else {
2085  signal[izone]->SetParLimits(4, 0.1, 1.0);
2086  signal[izone]->SetParLimits(6, 0.035, 0.3);
2087  }
2088  }
2089  signal[izone]->SetParNames("const", "mean", "sigma","Width","MP","Area","GSigma");
2090  signall[izone]->Fit(signal[izone], "0R+");
2091 
2092  signall[izone]->GetXaxis()->SetRangeUser(xmn,xmx);
2093  for (int k=0; k<nsgpr; k++) {
2094  fitres[k] = fitprm[k][jk] = signal[izone]->GetParameter(k);
2095  parserr[k] = signal[izone]->GetParError(k);
2096  }
2097 
2098  } else {
2099  double pedhtx = 0;
2100  for (unsigned i =0; i<sig_reg[ietafit][iphifit].size(); i++) {
2101  if (sig_reg[ietafit][iphifit][i] >gaupr[1]-3*gaupr[2] && sig_reg[ietafit][iphifit][i]<gaupr[1]+gaupr[2]) pedhtx++;
2102  }
2103 
2104  set_mean(gaupr[1], m_digiInput);
2105  set_sigma(gaupr[2], m_digiInput);
2106 
2107  TString name[nsgpr] = {"const", "mean", "sigma","Width","MP","Area","GSigma"};
2108  double strt[nsgpr] = {0.9*pedhtx, gaupr[1], gaupr[2], fitprm[3][jk], fitprm[4][jk], signall[izone]->GetEntries(), fitprm[6][jk]};
2109  double alowmn[nsgpr] = {0.1*pedhtx-0.1, gaupr[1]-0.1, gaupr[2]-0.1,0.07, 0.2*strt[4], 0.1*strt[5], 0.2*strt[6]};
2110  double ahighmn[nsgpr] ={1.2*pedhtx+0.1, gaupr[1]+0.1, gaupr[2]+0.1,0.20, 2.5*strt[4], 1.5*strt[5], 2.2*strt[6]};
2111  double step[nsgpr] = {1.0, 0.0, 0.0, 0.0, 0.001, 1.0, 0.002};
2112 
2113  TMinuit *gMinuit = new TMinuit(nsgpr);
2114  gMinuit->SetFCN(fcnsg);
2115 
2116  double arglist[10];
2117  int ierflg = 0;
2118  arglist[0] =0.5;
2119  gMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
2120 
2121  for (int k=0; k<nsgpr; k++) {
2122  gMinuit->mnparm(k, name[k], strt[k], step[k], alowmn[k], ahighmn[k],ierflg);
2123  }
2124 
2125  arglist[0] = 0;
2126  gMinuit->mnexcm("SIMPLEX", arglist, 0, ierflg);
2127 
2128  arglist[0] = 0;
2129  gMinuit->mnexcm("IMPROVE", arglist, 0, ierflg);
2130 
2131  TString chnam;
2132  double parv,err,xlo,xup, plerr, mierr, eparab, gcc;
2133  int iuit;
2134 
2135  for (int k=0; k<nsgpr; k++) {
2136  if (step[k] >-10) {
2137  gMinuit->mnpout(k, chnam, parv, err, xlo, xup, iuit);
2138  gMinuit->mnerrs(k, plerr, mierr, eparab, gcc);
2139  if (k==0 || k==5) {
2140  fitres[k] = parv*binwid;
2141  parserr[k]= err*binwid;
2142  } else {
2143  fitres[k] = parv;
2144  parserr[k]= err;
2145  }
2146 
2147  }
2148  }
2149 
2150  delete gMinuit;
2151  }
2152 
2153  // if (iijj==0) {
2154  // signall[izone]->Draw("same");
2155  // } else {
2156  signall[izone]->Draw();
2157  // }
2158 
2159  sprintf(temp, "pedfun_%i",izone);
2160  pedfun[izone] = new TF1(temp, gausX, xmn, xmx, nbgpr);
2161  pedfun[izone]->SetParameters(fitres);
2162  pedfun[izone]->SetLineColor(3);
2163  pedfun[izone]->SetLineWidth(1);
2164  pedfun[izone]->Draw("same");
2165 
2166  sprintf(temp, "signalfun_%i",izone);
2167  sigfun[izone] = new TF1(temp, langaufun, xmn, xmx, nsgpr-nbgpr);
2168  sigfun[izone]->SetParameters(&fitres[3]);
2169  sigfun[izone]->SetLineWidth(1);
2170  sigfun[izone]->SetLineColor(4);
2171  sigfun[izone]->Draw("same");
2172 
2173  sprintf(temp, "total_%i",izone);
2174  signalx[izone] = new TF1(temp, totalfunc, xmn, xmx, nsgpr);
2175  signalx[izone]->SetParameters(fitres);
2176  signalx[izone]->SetLineWidth(1);
2177  signalx[izone]->Draw("same");
2178 
2179  int kl = (jk<15) ? jk+1 : 14-jk;
2180 
2181  cout<<"histinfo"<<iijj<<" fit "
2182  <<std::setw(3)<< kl<<" "
2183  <<std::setw(3)<< ij+1<<" "
2184  <<std::setw(5)<<pedstll[izone]->GetEntries()<<" "
2185  <<std::setw(6)<<pedstll[izone]->GetMean()<<" "
2186  <<std::setw(6)<<pedstll[izone]->GetRMS()<<" "
2187  <<std::setw(5)<<signall[izone]->GetEntries()<<" "
2188  <<std::setw(6)<<signall[izone]->GetMean()<<" "
2189  <<std::setw(6)<<signall[izone]->GetRMS()<<" "
2190  <<std::setw(6)<< signal[izone]->GetChisquare()<<" "
2191  <<std::setw(3)<< signal[izone]->GetNDF()<<endl;
2192 
2193  file_out<<"histinfo"<<iijj<<" fit "
2194  <<std::setw(3)<< kl<<" "
2195  <<std::setw(3)<< ij+1<<" "
2196  <<std::setw(5)<<pedstll[izone]->GetEntries()<<" "
2197  <<std::setw(6)<<pedstll[izone]->GetMean()<<" "
2198  <<std::setw(6)<<pedstll[izone]->GetRMS()<<" "
2199  <<std::setw(5)<<signall[izone]->GetEntries()<<" "
2200  <<std::setw(6)<<signall[izone]->GetMean()<<" "
2201  <<std::setw(6)<<signall[izone]->GetRMS()<<" "
2202  <<std::setw(6)<< signal[izone]->GetChisquare()<<" "
2203  <<std::setw(3)<< signal[izone]->GetNDF()<<endl;
2204 
2205  file_out <<"fitres x"<<iijj<<" "<<kl<<" "<<ij+1<<" "<< fitres[0]<<" "<< fitres[1]<<" "<< fitres[2]<<" "<< fitres[3]<<" "<< fitres[4]<<" "<< fitres[5]<<" "<< fitres[6]<<endl;
2206  file_out <<"parserr"<<iijj<<" "<<kl<<" "<<ij+1<<" "<< parserr[0]<<" "<< parserr[1]<<" "<< parserr[2]<<" "<< parserr[3]<<" "<< parserr[4]<<" "<< parserr[5]<<" "<< parserr[6]<<endl;
2207 
2208  double diff=fitres[4]-fitres[1];
2209  if (diff <=0) diff = 0.000001;
2210  double error=parserr[4]*parserr[4]+parer[2]*parer[2];
2211  error = pow(error,0.5);
2212 
2213  int ieta = (jk<15) ? (15+jk) : (29-jk);
2214  int ifl = nphimx*ieta + ij;
2215 
2216  if (iijj==3) {
2217  ped_evt->Fill(ifl,pedstll[izone]->GetEntries());
2218  ped_mean->Fill(ifl,gaupr[1]);
2219  ped_width->Fill(ifl,gaupr[2]);
2220  fit_chi->Fill(ifl,signal[izone]->GetChisquare());
2221  sig_evt->Fill(ifl, signall[izone]->GetEntries());
2222  fit_sigevt->Fill(ifl, fitres[5]);
2223  fit_bkgevt->Fill(ifl, fitres[0]*sqrt(2*acos(-1.))*gaupr[2]);
2224  sig_mean->Fill(ifl, fitres[4]);
2225  sig_diff->Fill(ifl, fitres[4]-fitres[1]);
2226  sig_width->Fill(ifl, fitres[3]);
2227  sig_sigma->Fill(ifl, fitres[6]);
2228  sig_meanerr->Fill(ifl, parserr[4]);
2229  if (fitres[4]-fitres[1] !=0) sig_meanerrp->Fill(ifl, 100*parserr[4]/(fitres[4]-fitres[1]));
2230  if (gaupr[2]!=0) sig_signf->Fill(ifl,(fitres[4]-fitres[1])/gaupr[2]);
2231 
2232  ped_statmean->Fill(ifl,pedstll[izone]->GetMean());
2233  sig_statmean->Fill(ifl,signall[izone]->GetMean());
2234  ped_rms->Fill(ifl,pedstll[izone]->GetRMS());
2235  sig_rms->Fill(ifl,signall[izone]->GetRMS());
2236  }
2237 
2238  if ((iijj==2) || (iijj==3) || (iijj==1)) {
2239  if (signall[izone]->GetEntries() >5 && fitres[4]>0.1) {
2240  //GMA need to put this==1 in future
2241  float fact=0.812;
2242  if (abs(kl)<=4) fact=0.895;
2243  if (!m_digiInput) fact *=0.19; //conversion factor for GeV/fC
2244 
2245  float fact2 = 0;
2246  if (iijj==2) fact2 = invang[jk][nphimx];
2247  if (iijj==3) fact2 = invang[jk][ij];
2248  if (iijj==1) fact2 = com_invang[jk][ij];
2249 
2250  float calibc = fact*fact2/(fitres[4]*signall[izone]->GetEntries());
2251  float caliberr= TMath::Abs(calibc*parserr[4]/max(0.001,fitres[4]));
2252 
2253  if (iijj==2) {
2254  int ieta = (jk<15) ? jk+1 : 14-jk;
2255  mean_phi_hst->Fill(ieta, calibc);
2256  mean_phi_hst->SetBinError(mean_phi_hst->FindBin(ieta), caliberr);
2257  file_out<<"intieta "<<jk<<" "<<ij<<" "<<ieta<<" "<<mean_phi_hst->FindBin(double(ieta))<<" "<<calibc<<" "<<caliberr<<endl;
2258  } else if (iijj==3) {
2259  const_eta[jk]->Fill(ij+1,calibc);
2260  const_eta[jk]->SetBinError(const_eta[jk]->FindBin(ij+1), caliberr);
2261 
2262  peak_eta[jk]->Fill(ij+1,fitres[4]);
2263  peak_eta[jk]->SetBinError(peak_eta[jk]->FindBin(ij+1),parserr[4]);
2264 
2265  int ieta = (jk<15) ? jk+1 : 14-jk;
2266  const_eta_phi->Fill(ieta, ij+1,calibc);
2267  file_out<<"intietax "<<jk<<" "<<ij<<" "<<ieta<<" "<<const_eta_phi->FindBin(ieta, ij+1)<<endl;
2268  if (caliberr >0) {
2269  const_eta_phi->SetBinError(const_eta_phi->FindBin(ieta, ij+1),caliberr);
2270 
2271  mean_eta[ij] +=calibc/(caliberr*caliberr);
2272  mean_phi[jk] +=calibc/(caliberr*caliberr);
2273 
2274  rms_eta[ij] +=1./(caliberr*caliberr);
2275  rms_phi[jk] +=1./(caliberr*caliberr);
2276 
2277  } else {
2278  const_eta_phi->SetBinError(const_eta_phi->FindBin(ieta, ij+1), 0.0);
2279  }
2280  } else if (iijj==1) {
2281  const_hpdrm[jk]->Fill(ij+1,calibc);
2282  const_hpdrm[jk]->SetBinError(const_hpdrm[jk]->FindBin(ij+1), caliberr);
2283 
2284  peak_hpdrm[jk]->Fill(ij+1,fitres[4]);
2285  peak_hpdrm[jk]->SetBinError(peak_hpdrm[jk]->FindBin(ij+1),parserr[4]);
2286  }
2287 
2288  file_out<<"HO 4 "<<iijj<<" "<< std::setw(3)<<kl<<" "<<std::setw(3)<<ij+1<<" "
2289  <<std::setw(7)<<calibc<<" "<<std::setw(7)<<caliberr<<endl;
2290  }
2291  }
2292 
2293  } else { //if (signall[izone]->GetEntries() >10) {
2294  signall[izone]->Draw();
2295  float varx = 0.000;
2296  int kl = (jk<15) ? jk+1 : 14-jk;
2297  file_out<<"histinfo"<<iijj<<" nof "
2298  <<std::setw(3)<< kl<<" "
2299  <<std::setw(3)<< ij+1<<" "
2300  <<std::setw(5)<<pedstll[izone]->GetEntries()<<" "
2301  <<std::setw(6)<<pedstll[izone]->GetMean()<<" "
2302  <<std::setw(6)<<pedstll[izone]->GetRMS()<<" "
2303  <<std::setw(5)<<signall[izone]->GetEntries()<<" "
2304  <<std::setw(6)<<signall[izone]->GetMean()<<" "
2305  <<std::setw(6)<<signall[izone]->GetRMS()<<" "
2306  <<std::setw(6)<< varx<<" "
2307  <<std::setw(3)<< varx<<endl;
2308 
2309  file_out <<"fitres x"<<iijj<<" "<<kl<<" "<<ij+1<<" "<< varx<<" "<< varx<<" "<< varx<<" "<< varx<<" "<< varx<<" "<< varx<<" "<< varx<<endl;
2310  file_out <<"parserr"<<iijj<<" "<<kl<<" "<<ij+1<<" "<< varx<<" "<< varx<<" "<< varx<<" "<< varx<<" "<< varx<<" "<< varx<<" "<< varx<<endl;
2311 
2312  }
2313  iiter++;
2314  if (iiter%nsample==0) {
2315  c0->Update();
2316 
2317  for (int kl=0; kl<nsample; kl++) {
2318  if (gx0[kl]) {delete gx0[kl];gx0[kl] = 0;}
2319  if (ped0fun[kl]) {delete ped0fun[kl];ped0fun[kl] = 0;}
2320  if (signal[kl]) {delete signal[kl];signal[kl] = 0;}
2321  if (pedfun[kl]) {delete pedfun[kl];pedfun[kl] = 0;}
2322  if (sigfun[kl]) {delete sigfun[kl];sigfun[kl] = 0;}
2323  if (signalx[kl]) {delete signalx[kl];signalx[kl] = 0;}
2324  if (signall[kl]) {delete signall[kl];signall[kl] = 0;}
2325  if (pedstll[kl]) {delete pedstll[kl];pedstll[kl] = 0;}
2326  }
2327 
2328  }
2329  } //for (int jk=0; jk<netamx; jk++) {
2330  } //for (int ij=0; ij<nphimx; ij++) {
2331 
2332  // if (iijj==0) {
2333  // sprintf(out_file, "comb_hosig_allring_%i.jpg", irunold);
2334  // c0->SaveAs(out_file);
2335  // iiter = 0;
2336  // } else {
2337  // // c0->Update();
2338  // }
2339 
2340  // iiter = 0;
2341  } //end of iijj
2342  if (iiter%nsample!=0) {
2343  c0->Update();
2344  for (int kl=0; kl<nsample; kl++) {
2345  if (gx0[kl]) {delete gx0[kl];gx0[kl] = 0;}
2346  if (ped0fun[kl]) {delete ped0fun[kl];ped0fun[kl] = 0;}
2347  if (signal[kl]) {delete signal[kl];signal[kl] = 0;}
2348  if (pedfun[kl]) {delete pedfun[kl];pedfun[kl] = 0;}
2349  if (sigfun[kl]) {delete sigfun[kl];sigfun[kl] = 0;}
2350  if (signalx[kl]) {delete signalx[kl];signalx[kl] = 0;}
2351  if (signall[kl]) {delete signall[kl];signall[kl] = 0;}
2352  if (pedstll[kl]) {delete pedstll[kl];pedstll[kl] = 0;}
2353  }
2354  }
2355 
2356  delete c0;
2357 
2358  xsiz = 600; //int xsiz = 600;
2359  ysiz = 800; //int ysiz = 800;
2360 
2361  gStyle->SetTitleFontSize(0.05);
2362  gStyle->SetTitleSize(0.025,"XYZ");
2363  gStyle->SetLabelSize(0.025,"XYZ");
2364  gStyle->SetStatFontSize(.045);
2365 
2366  gStyle->SetOptStat(0);
2367  ps.NewPage(); TCanvas *c1 = new TCanvas("c1", " Pedestal vs signal", xsiz, ysiz);
2368  ped_evt->Draw(); c1->Update();
2369 
2370  ps.NewPage();
2371  ped_statmean->Draw(); c1->Update();
2372 
2373  ps.NewPage();
2374  ped_rms->Draw(); c1->Update();
2375 
2376  ps.NewPage();
2377  ped_mean->Draw(); c1->Update();
2378 
2379  ps.NewPage();
2380  ped_width->Draw(); c1->Update();
2381 
2382  ps.NewPage();
2383  sig_evt->Draw(); c1->Update();
2384 
2385  ps.NewPage();
2386  sig_statmean->Draw(); c1->Update();
2387 
2388  ps.NewPage();
2389  sig_rms->Draw(); c1->Update();
2390 
2391  ps.NewPage();
2392  fit_chi->Draw(); c1->Update();
2393 
2394  ps.NewPage();
2395  fit_sigevt->Draw(); c1->Update();
2396 
2397  ps.NewPage();
2398  fit_bkgevt->Draw(); c1->Update();
2399 
2400  ps.NewPage();
2401  sig_mean->Draw(); c1->Update();
2402 
2403  ps.NewPage();
2404  sig_width->Draw(); c1->Update();
2405 
2406  ps.NewPage();
2407  sig_sigma->Draw(); c1->Update();
2408 
2409  ps.NewPage();
2410  sig_meanerr->Draw(); c1->Update();
2411 
2412  ps.NewPage();
2413  sig_meanerrp->Draw(); c1->Update();
2414 
2415  ps.NewPage();
2416  sig_signf->Draw(); c1->Update();
2417 
2418  ps.Close();
2419  delete c1;
2420 
2421  file_out.close();
2422 
2423  if (m_figure) {
2424  xsiz = 700;
2425  ysiz = 450;
2426 
2427  gStyle->SetTitleFontSize(0.09);
2428  gStyle->SetPadBottomMargin(0.17);
2429  gStyle->SetPadLeftMargin(0.18);
2430  gStyle->SetPadRightMargin(0.01);
2431  gStyle->SetOptLogy(0);
2432  gStyle->SetOptStat(0);
2433 
2434  TCanvas *c2 = new TCanvas("c2", "runfile", xsiz, ysiz);
2435  c2->Divide(5,3);
2436 
2437  for (int side=0; side <2; side++) {
2438  gStyle->SetNdivisions(303,"XY");
2439  gStyle->SetPadRightMargin(0.01);
2440  int nmn = 0;
2441  int nmx = netamx/2;
2442  if (side==1) {
2443  nmn = netamx/2;
2444  nmx = netamx;
2445  }
2446 
2447  int nzone = 0;
2448 
2449  for (int ij=nmn; ij<nmx; ij++) {
2450 
2451  c2->cd(nzone+1);
2452  const_eta[ij]->GetXaxis()->SetTitle("#phi index");
2453  const_eta[ij]->GetXaxis()->SetTitleSize(.08);
2454  const_eta[ij]->GetXaxis()->CenterTitle();
2455  const_eta[ij]->GetXaxis()->SetTitleOffset(0.9);
2456  const_eta[ij]->GetXaxis()->SetLabelSize(.085);
2457  const_eta[ij]->GetXaxis()->SetLabelOffset(.01);
2458 
2459  const_eta[ij]->GetYaxis()->SetLabelSize(.08);
2460  const_eta[ij]->GetYaxis()->SetLabelOffset(.01);
2461  if (m_digiInput) {
2462  const_eta[ij]->GetYaxis()->SetTitle("GeV/fC");
2463  } else {
2464  const_eta[ij]->GetYaxis()->SetTitle("GeV/MIP-GeV!!");
2465  }
2466 
2467  const_eta[ij]->GetYaxis()->SetTitleSize(.085);
2468  const_eta[ij]->GetYaxis()->CenterTitle();
2469  const_eta[ij]->GetYaxis()->SetTitleOffset(1.3);
2470  const_eta[ij]->SetMarkerSize(0.60);
2471  const_eta[ij]->SetMarkerColor(2);
2472  const_eta[ij]->SetMarkerStyle(20);
2473 
2474 
2475  const_eta[ij]->Draw();
2476  nzone++;
2477  }
2478 
2479  sprintf(out_file, "calibho_%i_side%i.eps", irunold, side);
2480  c2->SaveAs(out_file);
2481 
2482  sprintf(out_file, "calibho_%i_side%i.jpg", irunold, side);
2483  c2->SaveAs(out_file);
2484 
2485  nzone = 0;
2486  for (int ij=nmn; ij<nmx; ij++) {
2487  c2->cd(nzone+1);
2488  peak_eta[ij]->GetXaxis()->SetTitle("#phi index");
2489  peak_eta[ij]->GetXaxis()->SetTitleSize(.08);
2490  peak_eta[ij]->GetXaxis()->CenterTitle();
2491  peak_eta[ij]->GetXaxis()->SetTitleOffset(0.90);
2492  peak_eta[ij]->GetXaxis()->SetLabelSize(.08);
2493  peak_eta[ij]->GetXaxis()->SetLabelOffset(.01);
2494 
2495  peak_eta[ij]->GetYaxis()->SetLabelSize(.08);
2496  peak_eta[ij]->GetYaxis()->SetLabelOffset(.01);
2497  if (m_digiInput) {
2498  peak_eta[ij]->GetYaxis()->SetTitle("fC");
2499  } else {
2500  peak_eta[ij]->GetYaxis()->SetTitle("GeV");
2501  }
2502 
2503  peak_eta[ij]->GetYaxis()->SetTitleSize(.085);
2504  peak_eta[ij]->GetYaxis()->CenterTitle();
2505  peak_eta[ij]->GetYaxis()->SetTitleOffset(1.3);
2506 
2507  peak_eta[ij]->SetMarkerSize(0.60);
2508  peak_eta[ij]->SetMarkerColor(2);
2509  peak_eta[ij]->SetMarkerStyle(20);
2510 
2511  peak_eta[ij]->Draw();
2512  nzone++;
2513  }
2514 
2515  sprintf(out_file, "peakho_%i_side%i.eps", irunold, side);
2516  c2->SaveAs(out_file);
2517 
2518  sprintf(out_file, "peakho_%i_side%i.jpg", irunold, side);
2519  c2->SaveAs(out_file);
2520  }
2521  delete c2;
2522 
2523  // if (m_combined) {
2524  gStyle->SetTitleFontSize(0.045);
2525  gStyle->SetPadRightMargin(0.13);
2526  gStyle->SetPadBottomMargin(0.15);
2527  gStyle->SetPadLeftMargin(0.1);
2528  gStyle->SetOptStat(0);
2529  xsiz = 700;
2530  ysiz = 600;
2531  TCanvas *c1 = new TCanvas("c1", "Fitted const in each tower", xsiz, ysiz);
2532  const_eta_phi->GetXaxis()->SetTitle("#eta");
2533  const_eta_phi->GetXaxis()->SetTitleSize(0.065);
2534  const_eta_phi->GetXaxis()->SetTitleOffset(0.85); //6);
2535  const_eta_phi->GetXaxis()->CenterTitle();
2536  const_eta_phi->GetXaxis()->SetLabelSize(0.045);
2537  const_eta_phi->GetXaxis()->SetLabelOffset(0.01);
2538 
2539  const_eta_phi->GetYaxis()->SetTitle("#phi");
2540  const_eta_phi->GetYaxis()->SetTitleSize(0.075);
2541  const_eta_phi->GetYaxis()->SetTitleOffset(0.5);
2542  const_eta_phi->GetYaxis()->CenterTitle();
2543  const_eta_phi->GetYaxis()->SetLabelSize(0.045);
2544  const_eta_phi->GetYaxis()->SetLabelOffset(0.01);
2545 
2546  const_eta_phi->Draw("colz");
2547  sprintf(out_file, "high_hoconst_eta_phi_%i.jpg",irunold);
2548  c1->SaveAs(out_file);
2549 
2550  delete c1;
2551 
2552  for (int jk=0; jk<netamx; jk++) {
2553  int ieta = (jk<15) ? jk+1 : 14-jk;
2554  if (rms_phi[jk]>0) {
2555  mean_phi_ave->Fill(ieta, mean_phi[jk]/rms_phi[jk]);
2556  mean_phi_ave->SetBinError(mean_phi_ave->FindBin(ieta), pow(double(rms_phi[jk]), -0.5));
2557  }
2558  }
2559 
2560  for (int ij=0; ij<nphimx; ij++) {
2561  if (rms_eta[ij] >0) {
2562  mean_eta_ave->Fill(ij+1, mean_eta[ij]/rms_eta[ij]);
2563  mean_eta_ave->SetBinError(mean_eta_ave->FindBin(ij+1), pow(double(rms_eta[ij]), -0.5));
2564  }
2565  }
2566 
2567  ysiz =450;
2568  gStyle->SetPadLeftMargin(0.13);
2569  gStyle->SetPadRightMargin(0.03);
2570 
2571 
2572  TCanvas *c2y = new TCanvas("c2", "Avearge signal in eta and phi", xsiz, ysiz);
2573  c2y->Divide(2,1);
2574  mean_eta_ave->GetXaxis()->SetTitle("#phi");
2575  mean_eta_ave->GetXaxis()->SetTitleSize(0.085);
2576  mean_eta_ave->GetXaxis()->SetTitleOffset(0.65);
2577  mean_eta_ave->GetXaxis()->CenterTitle();
2578  mean_eta_ave->GetXaxis()->SetLabelSize(0.05);
2579  mean_eta_ave->GetXaxis()->SetLabelOffset(0.001);
2580 
2581  mean_eta_ave->GetYaxis()->SetTitle("Signal (GeV)/MIP");
2582  mean_eta_ave->GetYaxis()->SetTitleSize(0.055);
2583  mean_eta_ave->GetYaxis()->SetTitleOffset(1.3);
2584  mean_eta_ave->GetYaxis()->CenterTitle();
2585  mean_eta_ave->GetYaxis()->SetLabelSize(0.045);
2586  mean_eta_ave->GetYaxis()->SetLabelOffset(0.01);
2587  mean_eta_ave->SetMarkerSize(0.60);
2588  mean_eta_ave->SetMarkerColor(2);
2589  mean_eta_ave->SetMarkerStyle(20);
2590 
2591  c2y->cd(1); mean_eta_ave->Draw();
2592 
2593  mean_phi_ave->GetXaxis()->SetTitle("#eta");
2594  mean_phi_ave->GetXaxis()->SetTitleSize(0.085);
2595  mean_phi_ave->GetXaxis()->SetTitleOffset(0.65); //55);
2596  mean_phi_ave->GetXaxis()->CenterTitle();
2597  mean_phi_ave->GetXaxis()->SetLabelSize(0.05);
2598  mean_phi_ave->GetXaxis()->SetLabelOffset(0.001);
2599 
2600  mean_phi_ave->GetYaxis()->SetTitle("Signal (GeV)/MIP");
2601  mean_phi_ave->GetYaxis()->SetTitleSize(0.055);
2602  mean_phi_ave->GetYaxis()->SetTitleOffset(1.3);
2603  mean_phi_ave->GetYaxis()->CenterTitle();
2604  mean_phi_ave->GetYaxis()->SetLabelSize(0.045);
2605  mean_phi_ave->GetYaxis()->SetLabelOffset(0.01);
2606  mean_phi_ave->SetMarkerSize(0.60);
2607  mean_phi_ave->SetMarkerColor(2);
2608  mean_phi_ave->SetMarkerStyle(20);
2609 
2610  c2y->cd(2); mean_phi_ave->Draw();
2611 
2612  sprintf(out_file, "high_hoaverage_eta_phi_%i.jpg",irunold);
2613  c2y->SaveAs(out_file);
2614 
2615  delete c2y;
2616  // } else { //m_combined
2617 
2618  xsiz = 800;
2619  ysiz = 450;
2620  TCanvas *c3 = new TCanvas("c3", "Avearge signal in eta and phi", xsiz, ysiz);
2621  c3->Divide(2,1);
2622  mean_phi_hst->GetXaxis()->SetTitle("#eta");
2623  mean_phi_hst->GetXaxis()->SetTitleSize(0.065);
2624  mean_phi_hst->GetXaxis()->SetTitleOffset(0.9);
2625  mean_phi_hst->GetXaxis()->CenterTitle();
2626  mean_phi_hst->GetXaxis()->SetLabelSize(0.065);
2627  mean_phi_hst->GetXaxis()->SetLabelOffset(0.001);
2628 
2629  mean_phi_hst->GetYaxis()->SetTitle("GeV/MIP");
2630  mean_phi_hst->GetYaxis()->SetTitleSize(0.055);
2631  mean_phi_hst->GetYaxis()->SetTitleOffset(0.9);
2632  mean_phi_hst->GetYaxis()->CenterTitle();
2633  mean_phi_hst->GetYaxis()->SetLabelSize(0.065);
2634  mean_phi_hst->GetYaxis()->SetLabelOffset(0.01);
2635 
2636  mean_phi_hst->SetMarkerColor(4);
2637  mean_phi_hst->SetMarkerSize(0.8);
2638  mean_phi_hst->SetMarkerStyle(20);
2639  mean_phi_hst->Draw();
2640 
2641  sprintf(out_file, "low_mean_phi_hst_%i.jpg",irunold);
2642  c3->SaveAs(out_file);
2643 
2644  delete c3;
2645 
2646  // } //m_combined
2647 
2648 
2649  gStyle->SetOptLogy(1);
2650  gStyle->SetPadTopMargin(.1);
2651  gStyle->SetPadLeftMargin(.15);
2652  xsiz = 800;
2653  ysiz = 500;
2654  TCanvas *c0x = new TCanvas("c0x", "Signal in each ring", xsiz, ysiz);
2655 
2656  c0x->Divide(3,2);
2657  for (int ij=0; ij<ringmx; ij++) {
2658  int iread = (ij==2) ? routmx : rout12mx;
2659  if (m_digiInput) {
2660  com_sigrsg[ij][iread]->GetXaxis()->SetTitle("Signal/ped (fC)");
2661  } else {
2662  com_sigrsg[ij][iread]->GetXaxis()->SetTitle("Signal/ped (GeV)");
2663  }
2664  com_sigrsg[ij][iread]->GetXaxis()->SetTitleSize(0.060);
2665  com_sigrsg[ij][iread]->GetXaxis()->SetTitleOffset(1.05);
2666  com_sigrsg[ij][iread]->GetXaxis()->CenterTitle();
2667  com_sigrsg[ij][iread]->GetXaxis()->SetLabelSize(0.065);
2668  com_sigrsg[ij][iread]->GetXaxis()->SetLabelOffset(0.01);
2669 
2670  com_sigrsg[ij][iread]->GetYaxis()->SetLabelSize(0.065);
2671  com_sigrsg[ij][iread]->GetYaxis()->SetLabelOffset(0.01);
2672 
2673 
2674  com_sigrsg[ij][iread]->SetLineWidth(3);
2675  com_sigrsg[ij][iread]->SetLineColor(4);
2676 
2677  c0x->cd(ij+1); com_sigrsg[ij][iread]->Draw();
2678 
2679  com_crossg[ij][iread]->SetLineWidth(2);
2680  com_crossg[ij][iread]->SetLineColor(2);
2681  com_crossg[ij][iread]->Draw("same");
2682  }
2683  sprintf(out_file, "hosig_ring_%i.jpg",irunold);
2684  c0x->SaveAs(out_file);
2685  delete c0x;
2686 
2687  gStyle->SetTitleFontSize(0.06);
2688  gStyle->SetOptStat(0);
2689  gStyle->SetOptLogy(0);
2690 
2691  TCanvas *c0 = new TCanvas("c0", "Signal in each ring", xsiz, ysiz);
2692 
2693  c0->Divide(3,2);
2694  for (int jk=0; jk<ringmx; jk++) {
2695  peak_hpdrm[jk]->GetXaxis()->SetTitle("RM #");
2696  peak_hpdrm[jk]->GetXaxis()->SetTitleSize(0.070);
2697  peak_hpdrm[jk]->GetXaxis()->SetTitleOffset(1.0);
2698  peak_hpdrm[jk]->GetXaxis()->CenterTitle();
2699  peak_hpdrm[jk]->GetXaxis()->SetLabelSize(0.065);
2700  peak_hpdrm[jk]->GetXaxis()->SetLabelOffset(0.01);
2701 
2702  peak_hpdrm[jk]->GetYaxis()->SetTitle("Peak(GeV)/MIP");
2703 
2704  peak_hpdrm[jk]->GetYaxis()->SetTitleSize(0.07);
2705  peak_hpdrm[jk]->GetYaxis()->SetTitleOffset(1.3);
2706  peak_hpdrm[jk]->GetYaxis()->CenterTitle();
2707  peak_hpdrm[jk]->GetYaxis()->SetLabelSize(0.065);
2708  peak_hpdrm[jk]->GetYaxis()->SetLabelOffset(0.01);
2709  // peak_hpdrm[jk]->SetLineWidth(3);
2710  // peak_hpdrm[jk]->SetLineColor(4);
2711  peak_hpdrm[jk]->SetMarkerSize(0.60);
2712  peak_hpdrm[jk]->SetMarkerColor(2);
2713  peak_hpdrm[jk]->SetMarkerStyle(20);
2714 
2715 
2716  c0->cd(jk+1); peak_hpdrm[jk]->Draw();
2717  }
2718  sprintf(out_file, "comb_peak_hpdrm_%i.jpg",irunold);
2719  c0->SaveAs(out_file);
2720 
2721  delete c0;
2722 
2723  TCanvas *c1y = new TCanvas("c1y", "Signal in each ring", xsiz, ysiz);
2724 
2725  c1y->Divide(3,2);
2726  for (int jk=0; jk<ringmx; jk++) {
2727  const_hpdrm[jk]->GetXaxis()->SetTitle("RM #");
2728  const_hpdrm[jk]->GetXaxis()->SetTitleSize(0.070);
2729  const_hpdrm[jk]->GetXaxis()->SetTitleOffset(1.3);
2730  const_hpdrm[jk]->GetXaxis()->CenterTitle();
2731  const_hpdrm[jk]->GetXaxis()->SetLabelSize(0.065);
2732  const_hpdrm[jk]->GetXaxis()->SetLabelOffset(0.01);
2733 
2734  if (m_digiInput) {
2735  const_hpdrm[jk]->GetYaxis()->SetTitle("Peak(fC)");
2736  } else {
2737  const_hpdrm[jk]->GetYaxis()->SetTitle("Peak(GeV)");
2738  }
2739  const_hpdrm[jk]->GetYaxis()->SetTitleSize(0.065);
2740  const_hpdrm[jk]->GetYaxis()->SetTitleOffset(1.0);
2741  const_hpdrm[jk]->GetYaxis()->CenterTitle();
2742  const_hpdrm[jk]->GetYaxis()->SetLabelSize(0.065);
2743  const_hpdrm[jk]->GetYaxis()->SetLabelOffset(0.01);
2744  // const_hpdrm[jk]->SetLineWidth(3);
2745  // const_hpdrm[jk]->SetLineColor(4);
2746  const_hpdrm[jk]->SetMarkerSize(0.60);
2747  const_hpdrm[jk]->SetMarkerColor(2);
2748  const_hpdrm[jk]->SetMarkerStyle(20);
2749 
2750  c1y->cd(jk+1); const_hpdrm[jk]->Draw();
2751  }
2752 
2753  sprintf(out_file, "comb_const_hpdrm_%i.jpg",irunold);
2754  c1y->SaveAs(out_file);
2755 
2756  delete c1y;
2757 
2758  } //if (m_figure) {
2759 
2760  // ps.Close();
2761  // file_out.close();
2762 
2763  }// if (m_constant){
2764 
2765 
2766  if (m_figure) {
2767  for (int ij=0; ij<nphimx; ij++) {
2768  for (int jk=0; jk<netamx; jk++) {
2769  stat_eta[jk]->Fill(ij+1,sigrsg[jk][ij]->GetEntries());
2770  statmn_eta[jk]->Fill(ij+1,sigrsg[jk][ij]->GetMean());
2771  }
2772  }
2773 
2774  xsiz = 700;
2775  ysiz = 450;
2776  gStyle->SetTitleFontSize(0.09);
2777  gStyle->SetPadBottomMargin(0.14);
2778  gStyle->SetPadLeftMargin(0.17);
2779  gStyle->SetPadRightMargin(0.01);
2780  gStyle->SetNdivisions(303,"XY");
2781  gStyle->SetOptLogy(1);
2782 
2783  TCanvas *c2x = new TCanvas("c2x", "runfile", xsiz, ysiz);
2784  c2x->Divide(5,3);
2785  for (int side=0; side <2; side++) {
2786  int nmn = 0;
2787  int nmx = netamx/2;
2788  if (side==1) {
2789  nmn = netamx/2;
2790  nmx = netamx;
2791  }
2792  int nzone = 0;
2793  char name[200];
2794 
2795  for (int ij=nmn; ij<nmx; ij++) {
2796  int ieta = (ij<15) ? ij+1 : 14-ij;
2797  c2x->cd(nzone+1);
2798  if (m_digiInput) {
2799  sprintf(name,"fC(#eta=%i)",ieta);
2800  } else {
2801  sprintf(name,"GeV(#eta=%i)",ieta);
2802  }
2803  sigrsg[ij][nphimx]->GetXaxis()->SetTitle(name);
2804  sigrsg[ij][nphimx]->GetXaxis()->SetTitleSize(.08);
2805  sigrsg[ij][nphimx]->GetXaxis()->CenterTitle();
2806  sigrsg[ij][nphimx]->GetXaxis()->SetTitleOffset(0.90);
2807  sigrsg[ij][nphimx]->GetXaxis()->SetLabelSize(.08);
2808  sigrsg[ij][nphimx]->GetXaxis()->SetLabelOffset(.01);
2809 
2810  sigrsg[ij][nphimx]->GetYaxis()->SetLabelSize(.08);
2811  sigrsg[ij][nphimx]->GetYaxis()->SetLabelOffset(.01);
2812  sigrsg[ij][nphimx]->SetLineWidth(2);
2813  sigrsg[ij][nphimx]->SetLineColor(4);
2814  sigrsg[ij][nphimx]->Draw();
2815  crossg[ij][nphimx]->SetLineWidth(2);
2816  crossg[ij][nphimx]->SetLineColor(2);
2817  crossg[ij][nphimx]->Draw("same");
2818  nzone++;
2819  }
2820 
2821  sprintf(out_file, "sig_ho_%i_side%i.eps", irunold, side);
2822  c2x->SaveAs(out_file);
2823 
2824  sprintf(out_file, "sig_ho_%i_side%i.jpg", irunold, side);
2825  c2x->SaveAs(out_file);
2826  }
2827 
2828  gStyle->SetOptLogy(0);
2829  c2x = new TCanvas("c2x", "runfile", xsiz, ysiz);
2830  c2x->Divide(5,3);
2831  for (int side=0; side <2; side++) {
2832  int nmn = 0;
2833  int nmx = netamx/2;
2834  if (side==1) {
2835  nmn = netamx/2;
2836  nmx = netamx;
2837  }
2838  int nzone = 0;
2839 
2840  nzone = 0;
2841  for (int ij=nmn; ij<nmx; ij++) {
2842  c2x->cd(nzone+1);
2843  statmn_eta[ij]->SetLineWidth(2);
2844  statmn_eta[ij]->SetLineColor(4);
2845  statmn_eta[ij]->GetXaxis()->SetTitle("#phi index");
2846  statmn_eta[ij]->GetXaxis()->SetTitleSize(.08);
2847  statmn_eta[ij]->GetXaxis()->CenterTitle();
2848  statmn_eta[ij]->GetXaxis()->SetTitleOffset(0.9);
2849  statmn_eta[ij]->GetYaxis()->SetLabelSize(.08);
2850  statmn_eta[ij]->GetYaxis()->SetLabelOffset(.01);
2851  statmn_eta[ij]->GetXaxis()->SetLabelSize(.08);
2852  statmn_eta[ij]->GetXaxis()->SetLabelOffset(.01);
2853  if (m_digiInput) {
2854  statmn_eta[ij]->GetYaxis()->SetTitle("fC");
2855  } else {
2856  statmn_eta[ij]->GetYaxis()->SetTitle("GeV");
2857  }
2858  statmn_eta[ij]->GetYaxis()->SetTitleSize(.075);
2859  statmn_eta[ij]->GetYaxis()->CenterTitle();
2860  statmn_eta[ij]->GetYaxis()->SetTitleOffset(1.30);
2861 
2862  statmn_eta[ij]->Draw();
2863  nzone++;
2864  }
2865 
2866  sprintf(out_file, "statmnho_%i_side%i.eps", irunold, side);
2867  c2x->SaveAs(out_file);
2868 
2869  sprintf(out_file, "statmnho_%i_side%i.jpg", irunold, side);
2870  c2x->SaveAs(out_file);
2871 
2872  gStyle->SetOptLogy(1);
2873  gStyle->SetNdivisions(203,"XY");
2874 
2875  nzone = 0;
2876  for (int ij=nmn; ij<nmx; ij++) {
2877  c2x->cd(nzone+1);
2878  stat_eta[ij]->SetLineWidth(2);
2879  stat_eta[ij]->SetLineColor(4);
2880  stat_eta[ij]->GetXaxis()->SetTitle("#phi index");
2881  stat_eta[ij]->GetXaxis()->SetTitleSize(.08);
2882  stat_eta[ij]->GetXaxis()->CenterTitle();
2883  stat_eta[ij]->GetXaxis()->SetTitleOffset(0.80);
2884  stat_eta[ij]->GetXaxis()->SetLabelSize(.08);
2885  stat_eta[ij]->GetXaxis()->SetLabelOffset(.01);
2886  stat_eta[ij]->GetYaxis()->SetLabelSize(.08);
2887  stat_eta[ij]->GetYaxis()->SetLabelOffset(.01);
2888 
2889  stat_eta[ij]->Draw();
2890  nzone++;
2891  }
2892 
2893  sprintf(out_file, "statho_%i_side%i.eps", irunold, side);
2894  c2x->SaveAs(out_file);
2895 
2896  sprintf(out_file, "statho_%i_side%i.jpg", irunold, side);
2897  c2x->SaveAs(out_file);
2898  }
2899  delete c2x;
2900 
2901  } //if (m_figure) {
2902 
2903  if (!m_constant) { //m_constant
2904  for (int j=0; j<netamx; j++) {
2905  for (int i=0; i<nphimx; i++) {
2906  if (crossg[j][i]) { delete crossg[j][i];}
2907  if (sigrsg[j][i]) { delete sigrsg[j][i];}
2908  }
2909  }
2910  }
2911 }
2912 
2913 //define this as a plug-in
2914 DEFINE_FWK_MODULE(HOCalibAnalyzer);
2915 
2916 
2917 
edm::EventID::run
RunNumber_t run() const
Definition: EventID.h:39
fitWZ.arglist
tuple arglist
Definition: fitWZ.py:38
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
edm::EventID::event
EventNumber_t event() const
Definition: EventID.h:41
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
i
int i
Definition: DBlmapReader.cc:9
HOCalibAnalyzer::ringmx
static const int ringmx
Definition: HOCalibAnalyzer.cc:348
HOCalibAnalyzer::theRootFileName
std::string theRootFileName
Definition: HOCalibAnalyzer.cc:230
npixrigh
static const int npixrigh[21]
Definition: HOCalibAnalyzer.cc:90
mapx2
static const int mapx2[6][3]
Definition: HOCalibAnalyzer.cc:73
npixleft
static const int npixleft[21]
Definition: HOCalibAnalyzer.cc:89
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:449
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
HOCalibAnalyzer::HOCalibAnalyzer
HOCalibAnalyzer(const edm::ParameterSet &)
Definition: HOCalibAnalyzer.cc:463
HOCalibAnalyzer::rout12mx
static const int rout12mx
Definition: HOCalibAnalyzer.cc:351
HOCalibAnalyzer::beginJob
virtual void beginJob() override
Definition: HOCalibAnalyzer.cc:1490
phimap
static const int phimap[4][21]
Definition: HOCalibAnalyzer.cc:83
TFileService::make
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
HOCalibAnalyzer::endJob
virtual void endJob() override
Definition: HOCalibAnalyzer.cc:1496
mapx1
static const int mapx1[6][3]
Definition: HOCalibAnalyzer.cc:71
HOCalibAnalyzer::analyze
virtual void analyze(const edm::Event &, const edm::EventSetup &) override
Definition: HOCalibAnalyzer.cc:917
npixriup
static const int npixriup[21]
Definition: HOCalibAnalyzer.cc:94
HOCalibAnalyzer::neffip
static const int neffip
Definition: HOCalibAnalyzer.cc:352
HOCalibAnalyzer::hoCalibVariableCollectionTag
edm::InputTag hoCalibVariableCollectionTag
Definition: HOCalibAnalyzer.cc:439
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
set_sigma
void set_sigma(double &x, bool mdigi)
Definition: HOCalibAnalyzer.cc:205
bk::beginJob
void beginJob()
Definition: Breakpoints.cc:15
HOCalibAnalyzer::theoutputtxtFile
std::string theoutputtxtFile
Definition: HOCalibAnalyzer.cc:231
iEvent
int iEvent
Definition: GenABIO.cc:230
np
int np
Definition: AMPTWrapper.h:33
totalfunc
Double_t totalfunc(Double_t *x, Double_t *par)
Definition: HOCalibAnalyzer.cc:169
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
sig_reg
vector< float > sig_reg[netamx][nphimx+1]
Definition: HOCalibAnalyzer.cc:103
HOCalibAnalyzer::theoutputpsFile
std::string theoutputpsFile
Definition: HOCalibAnalyzer.cc:232
npixribt
static const int npixribt[21]
Definition: HOCalibAnalyzer.cc:92
Abs
T Abs(T a)
Definition: MathUtil.h:49
iphifit
int iphifit
Definition: HOCalibAnalyzer.cc:102
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
j
int j
Definition: DBlmapReader.cc:9
fcnsg
void fcnsg(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t flag)
Definition: HOCalibAnalyzer.cc:184
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
nbgpr
const int nbgpr
Definition: HOCalibAnalyzer.cc:98
min
T min(T a, T b)
Definition: MathUtil.h:58
npixleup
static const int npixleup[21]
Definition: HOCalibAnalyzer.cc:93
gausX
Double_t gausX(Double_t *x, Double_t *par)
Definition: HOCalibAnalyzer.cc:112
nsgpr
const int nsgpr
Definition: HOCalibAnalyzer.cc:99
HOCalibAnalyzer::routmx
static const int routmx
Definition: HOCalibAnalyzer.cc:350
fact
const double fact
Definition: NuclearInteractionFTFSimulator.cc:61
mapx0m
static const int mapx0m[9][2]
Definition: HOCalibAnalyzer.cc:76
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:55
set_mean
void set_mean(double &x, bool mdigi)
Definition: HOCalibAnalyzer.cc:195
HOCalibAnalyzer::sectmx
static const int sectmx
Definition: HOCalibAnalyzer.cc:349
relval_parameters_module.step
string step
Definition: relval_parameters_module.py:36
cro_ssg
vector< float > cro_ssg[netamx][nphimx+1]
Definition: HOCalibAnalyzer.cc:104
edm::EventBase::id
edm::EventID id() const
Definition: EventBase.h:60
HOCalibAnalyzer::tok_ho_
edm::EDGetTokenT< HOCalibVariableCollection > tok_ho_
Definition: HOCalibAnalyzer.cc:440
nphimx
const int nphimx
Definition: AlCaHOCalibProducer.cc:152
netamx
const int netamx
Definition: AlCaHOCalibProducer.cc:153
mapx0p
static const int mapx0p[9][2]
Definition: HOCalibAnalyzer.cc:75
fitWZ.gMinuit
tuple gMinuit
Definition: fitWZ.py:35
ietafit
int ietafit
Definition: HOCalibAnalyzer.cc:101
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
fcnbg
void fcnbg(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t flag)
Definition: HOCalibAnalyzer.cc:173
x
Definition: DDAxes.h:10
langaufun
Double_t langaufun(Double_t *x, Double_t *par)
Definition: HOCalibAnalyzer.cc:116
etamap
static const int etamap[4][21]
Definition: HOCalibAnalyzer.cc:78
relval_steps.k
list k
Definition: relval_steps.py:1518
npixlebt
static const int npixlebt[21]
Definition: HOCalibAnalyzer.cc:91
citk::for
for(const auto &isodef:isoDefs)
Definition: CITKPFIsolationSumProducer.cc:77
findQualityFiles.size
tuple size
Write out results.
Definition: findQualityFiles.py:442
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
cmsBatch.log
tuple log
Definition: cmsBatch.py:347
Generated for CMSSW Reference Manual by doxygen 1.8.5