CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/Calibration/HcalCalibAlgos/src/HOCalibAnalyzer.cc

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    HOCalibAnalyzer
00004 // Class:      HOCalibAnalyzer
00005 // 
00013 //
00014 // Original Author:  Gobinda Majumder
00015 //         Created:  Sat Jul  7 09:51:31 CEST 2007
00016 // $Id: HOCalibAnalyzer.cc,v 1.12 2012/09/27 20:32:04 wdd Exp $
00017 // $Id: HOCalibAnalyzer.cc,v 1.12 2012/09/27 20:32:04 wdd Exp $
00018 //
00019 //
00020 
00021 
00022 // system include files
00023 #include <memory>
00024 
00025 // user include files
00026 #include "FWCore/Framework/interface/Frameworkfwd.h"
00027 #include "FWCore/Framework/interface/EDAnalyzer.h"
00028 
00029 #include "FWCore/Framework/interface/Event.h"
00030 #include "FWCore/Framework/interface/MakerMacros.h"
00031 
00032 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00033 //#include "DataFormats/HOCalibHit/interface/HOCalibVariables.h"
00034 #include "DataFormats/HcalCalibObjects/interface/HOCalibVariables.h"
00035 
00036 #include "FWCore/ServiceRegistry/interface/Service.h"
00037 #include "FWCore/Utilities/interface/InputTag.h"
00038 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00039 
00040 #include "TMath.h"
00041 #include "TFile.h"
00042 #include "TH1F.h"
00043 #include "TH2F.h"
00044 #include "TTree.h"
00045 #include "TProfile.h"
00046 #include "TPostScript.h"
00047 #include "TCanvas.h"
00048 #include "TF1.h"
00049 #include "TStyle.h"
00050 #include "TMinuit.h"
00051 #include "TMath.h"
00052 
00053 #include <string>
00054 
00055 #include <iostream>
00056 #include <fstream>
00057 #include <iomanip>
00058 //#include <sstream>
00059 
00060 
00061 using namespace std;
00062 using namespace edm;
00063 
00064 
00065   //
00066   //  Look for nearby pixel through eta, phi informations for pixel cross-talk
00067   // 1. Look PIXEL code from (eta,phi)
00068   // 2. Go to nearby pixel code 
00069   // 3. Come back to (eta,phi) from pixel code
00070   // Though it works, it is a very ugly/crude way to get cross talk, need better algorithms
00071   //
00072 
00073 static const int mapx1[6][3]={{1,4,8}, {12,7,3}, {5,9,13}, {11,6,2}, {16,15,14}, {19,18,17}}; 
00074 
00075 static const int mapx2[6][3]={{1,4,8}, {12,7,3}, {5,9,13}, {11,6,2}, {16,15,14}, {-1,-1,-1}};
00076 
00077 static const int mapx0p[9][2]={{3,1}, {7,4}, {6,5},  {12,8}, {0,0}, {11,9}, {16,13}, {15,14}, {19,17}};
00078 static const int mapx0m[9][2]={{17,19}, {14,15}, {13,16}, {9,11}, {0,0}, {8,12}, {5,6}, {4,7}, {1,3}};
00079 
00080 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
00081                    {-1, 0,3,1, 0,2,3, 1,0,2, -1, 3,1,2, 4,4,4, 5,5,5, -1},    //etamap1 
00082                    {-1, 0,-1,0, 1,2,2, 1,3,5, -1, 5,3,6, 7,7,6, 8,-1,8, -1},  //etamap0p
00083                    {-1, 8,-1,8, 7,6,6, 7,5,3, -1, 3,5,2, 1,1,2, 0,-1,0, -1}}; //etamap0m
00084 
00085 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
00086                     {-1, 0,2,2, 1,0,1, 1,2,1, -1, 0,0,2, 2,1,0, 2,1,0, -1},    //phimap1
00087                     {-1, 1,-1,0, 1,1,0, 0,1,1, -1, 0,0,1, 1,0,0, 1,-1,0, -1},  //phimap0p
00088                     {-1, 0,-1,1, 0,0,1, 1,0,0, -1, 1,1,0, 0,1,1, 0,-1,1, -1}};  //phimap0m
00089 //swapped phi map for R0+/R0- (15/03/07)  
00090 
00091 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};
00092 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};
00093 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};
00094 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};
00095 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};
00096 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};
00097 
00098 static const int netamx=30;
00099 static const int nphimx =72;
00100   const int nbgpr = 3;
00101   const int nsgpr = 7;
00102 
00103 int ietafit;
00104 int iphifit;
00105 vector<float>sig_reg[netamx][nphimx+1];
00106 vector<float>cro_ssg[netamx][nphimx+1];
00107 
00108 
00109 //#define CORREL
00110 //
00111 // class decleration
00112 //
00113 
00114 Double_t gausX(Double_t* x, Double_t* par){
00115   return par[0]*(TMath::Gaus(x[0], par[1], par[2], kTRUE));
00116 }
00117 
00118 Double_t langaufun(Double_t *x, Double_t *par) {
00119 
00120   //Fit parameters:
00121   //par[0]*par[1]=Width (scale) parameter of Landau density
00122   //par[1]=Most Probable (MP, location) parameter of Landau density
00123   //par[2]=Total area (integral -inf to inf, normalization constant)
00124   //par[3]=Width (sigma) of convoluted Gaussian function
00125   //
00126   //In the Landau distribution (represented by the CERNLIB approximation), 
00127   //the maximum is located at x=-0.22278298 with the location parameter=0.
00128   //This shift is corrected within this function, so that the actual
00129   //maximum is identical to the MP parameter.
00130   // /*
00131   // Numeric constants
00132   Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
00133   Double_t mpshift  = -0.22278298;       // Landau maximum location
00134   
00135   // Control constants
00136   Double_t np = 100.0;      // number of convolution steps
00137   Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas
00138   
00139   // Variables
00140   Double_t xx;
00141   Double_t mpc;
00142   Double_t fland;
00143   Double_t sum = 0.0;
00144   Double_t xlow,xupp;
00145   Double_t step;
00146   Double_t i;
00147   
00148   
00149   // MP shift correction
00150   mpc = par[1] - mpshift * par[0]*par[1]; 
00151   
00152   // Range of convolution integral
00153   xlow = x[0] - sc * par[3];
00154   xupp = x[0] + sc * par[3];
00155   
00156   step = (xupp-xlow) / np;
00157   
00158   // Convolution integral of Landau and Gaussian by sum
00159   for(i=1.0; i<=np/2; i++) {
00160     xx = xlow + (i-.5) * step;
00161     fland = TMath::Landau(xx,mpc,par[0]*par[1], kTRUE); // / par[0];
00162     sum += fland * TMath::Gaus(x[0],xx,par[3]);
00163     xx = xupp - (i-.5) * step;
00164     fland = TMath::Landau(xx,mpc,par[0]*par[1], kTRUE); // / par[0];
00165     sum += fland * TMath::Gaus(x[0],xx,par[3]);
00166   }
00167   
00168   return (par[2] * step * sum * invsq2pi / par[3]);
00169 }
00170 
00171 Double_t totalfunc(Double_t* x, Double_t* par){
00172   return gausX(x, par) + langaufun(x, &par[3]);
00173 }
00174 
00175 void fcnbg(Int_t &npar, Double_t* gin, Double_t &f, Double_t* par, Int_t flag) {
00176   
00177   double fval = -par[0];
00178   for (unsigned i=0; i<cro_ssg[ietafit][iphifit].size(); i++) {
00179     double xval = (double)cro_ssg[ietafit][iphifit][i];
00180     fval +=log(max(1.e-30,par[0]*TMath::Gaus(xval, par[1], par[2], 1)));
00181     //    fval +=log(par[0]*TMath::Gaus(xval, par[1], par[2], 1));
00182   }
00183   f = -fval;
00184 }
00185 
00186 void fcnsg(Int_t &npar, Double_t* gin, Double_t &f, Double_t* par, Int_t flag) {
00187 
00188   double xval[2];
00189   double fval = -(par[0]+par[5]);
00190   for (unsigned i=0; i<sig_reg[ietafit][iphifit].size(); i++) {
00191     xval[0] = (double)sig_reg[ietafit][iphifit][i];
00192     fval +=log(totalfunc(xval, par));
00193   }
00194   f = -fval;
00195 }
00196 
00197 void set_mean(double& x, bool mdigi) {
00198   if(mdigi) {
00199     x = min(x, 0.5);
00200     x = max(x, -0.5);
00201   } else {
00202     x = min(x, 0.1);
00203     x = max(x, -0.1);
00204   }
00205 }
00206 
00207 void set_sigma(double& x, bool mdigi) {
00208   if(mdigi) {
00209     x = min(x, 1.2);
00210     x = max(x, -1.2);
00211   } else {
00212     x = min(x, 0.24);
00213     x = max(x, -0.24);
00214   }
00215 }
00216 
00217 
00218 
00219 class HOCalibAnalyzer : public edm::EDAnalyzer {
00220    public:
00221       explicit HOCalibAnalyzer(const edm::ParameterSet&);
00222       ~HOCalibAnalyzer();
00223 
00224 
00225    private:
00226 
00227       virtual void beginJob() ;
00228       virtual void analyze(const edm::Event&, const edm::EventSetup&);
00229       virtual void endJob() ;
00230 
00231       TFile* theFile;
00232       std::string theRootFileName;
00233       std::string theoutputtxtFile;
00234       std::string theoutputpsFile;
00235 
00236       bool m_hotime;
00237       bool m_hbtime;
00238       bool m_correl;
00239       bool m_checkmap;
00240       bool m_hbinfo;
00241       bool m_combined;
00242       bool m_constant;
00243       bool m_figure;
00244       bool m_digiInput;  
00245       bool m_cosmic; 
00246       bool m_histfit;
00247       bool m_pedsuppr;
00248       double m_sigma;
00249 
00250   static const int ncut = 13;
00251   static const int mypow_2_ncut = 8192; // 2^13, should be changed to match ncut
00252 
00253   int ipass;
00254  
00255   TTree* T1;
00256 
00257   TH1F* muonnm;
00258   TH1F* muonmm;
00259   TH1F* muonth;
00260   TH1F* muonph;    
00261   TH1F* muonch;    
00262 
00263   TH1F* sel_muonnm;
00264   TH1F* sel_muonmm;
00265   TH1F* sel_muonth;
00266   TH1F* sel_muonph;    
00267   TH1F* sel_muonch; 
00268 
00269   //  if (m_hotime) { // #ifdef HOTIME
00270   TProfile* hotime[netamx][nphimx];
00271   TProfile* hopedtime[netamx][nphimx];
00272   //  } //  if (m_hotime) #endif
00273   //  if (m_hbtime) { // #ifdef HBTIME
00274   TProfile* hbtime[netamx][nphimx];
00275   //  } // m_hbtime #endif
00276 
00277   //  double pedval[netamx][nphimx]; // ={0};
00278 
00279   //  if (m_correl) { // #ifdef CORREL
00280   TH1F* corrsglb[netamx][nphimx];
00281   TH1F* corrsgrb[netamx][nphimx];
00282   TH1F* corrsglu[netamx][nphimx];
00283   TH1F* corrsgru[netamx][nphimx];
00284   TH1F* corrsgall[netamx][nphimx];
00285 
00286   TH1F* corrsgl[netamx][nphimx];
00287   TH1F* corrsgr[netamx][nphimx];
00288 
00289   TH1F* mncorrsglb;
00290   TH1F* mncorrsgrb;
00291   TH1F* mncorrsglu;
00292   TH1F* mncorrsgru;
00293   TH1F* mncorrsgall;
00294 
00295   TH1F* mncorrsgl;
00296   TH1F* mncorrsgr;
00297 
00298   TH1F* rmscorrsglb;
00299   TH1F* rmscorrsgrb;
00300   TH1F* rmscorrsglu;
00301   TH1F* rmscorrsgru;
00302   TH1F* rmscorrsgall;
00303 
00304   TH1F* rmscorrsgl;
00305   TH1F* rmscorrsgr;
00306 
00307   TH1F* nevcorrsglb;
00308   TH1F* nevcorrsgrb;
00309   TH1F* nevcorrsglu;
00310   TH1F* nevcorrsgru;
00311   TH1F* nevcorrsgall;
00312 
00313   TH1F* nevcorrsgl;
00314   TH1F* nevcorrsgr;
00315 
00316   //  } //m_correl #endif
00317   //  if (m_checkmap) { // #ifdef CHECKMAP
00318   TH1F* corrsgc[netamx][nphimx];
00319   TH1F* mncorrsgc;
00320   TH1F* rmscorrsgc;
00321   TH1F* nevcorrsgc;
00322 
00323   //  } //m_checkmap #endif
00324 
00325   TH1F* sigrsg[netamx][nphimx+1];
00326   TH1F* crossg[netamx][nphimx+1];
00327   float invang[netamx][nphimx+1];
00328 
00329   TH1F* mnsigrsg;
00330   TH1F* mncrossg;
00331 
00332   TH1F* rmssigrsg;
00333   TH1F* rmscrossg;
00334 
00335   TH1F* nevsigrsg;
00336   TH1F* nevcrossg;
00337 
00338   TH1F* ho_sig2p[9];
00339   TH1F* ho_sig1p[9];
00340   TH1F* ho_sig00[9];
00341   TH1F* ho_sig1m[9];
00342   TH1F* ho_sig2m[9];
00343 
00344   //  if (m_hbinfo) { // #ifdef HBINFO
00345   TH1F* hbhe_sig[9];
00346   //  } // m_hbinfo #endif
00347 
00348 
00349   //  if (m_combined) { //#ifdef COMBINED
00350   static const int ringmx=5;
00351   static const int sectmx=12;
00352   static const int routmx=36;
00353   static const int rout12mx=24;
00354   static const int neffip=6;
00355 
00356   //  if (m_hotime) { // #ifdef HOTIME
00357   TProfile* com_hotime[ringmx][sectmx];
00358   TProfile* com_hopedtime[ringmx][sectmx];
00359   //  } //m_hotime #endif
00360   //  if (m_hbtime) { #ifdef HBTIME
00361   TProfile* com_hbtime[ringmx][sectmx];
00362     //  } // m_hbtime #endif
00363   //  if (m_correl) { //#ifdef CORREL
00364   TH1F* com_corrsglb[ringmx][sectmx];
00365   TH1F* com_corrsgrb[ringmx][sectmx];
00366   TH1F* com_corrsglu[ringmx][sectmx];
00367   TH1F* com_corrsgru[ringmx][sectmx];
00368   TH1F* com_corrsgall[ringmx][sectmx];
00369 
00370   TH1F* com_corrsgl[ringmx][sectmx];
00371   TH1F* com_corrsgr[ringmx][sectmx];
00372   //  } //m_correl #endif
00373   //  if (m_checkmap) { #ifdef CHECKMAP
00374   TH1F* com_corrsgc[ringmx][sectmx];
00375   //  } // m_checkmap #endif
00376 
00377   TH1F* com_sigrsg[ringmx][routmx+1];
00378   TH1F* com_crossg[ringmx][routmx+1];
00379   float com_invang[ringmx][routmx+1];
00380   
00381 
00382   //  } //m_combined #endif
00383 
00384   TH1F* ped_evt;
00385   TH1F* ped_mean;
00386   TH1F* ped_width;
00387   TH1F* fit_chi;
00388   TH1F* sig_evt;
00389   TH1F* fit_sigevt;
00390   TH1F* fit_bkgevt;
00391   TH1F* sig_mean;
00392   TH1F* sig_diff;
00393   TH1F* sig_width;
00394   TH1F* sig_sigma;
00395   TH1F* sig_meanerr;
00396   TH1F* sig_meanerrp;
00397   TH1F* sig_signf;
00398 
00399   TH1F* ped_statmean;
00400   TH1F* sig_statmean;
00401   TH1F* ped_rms;
00402   TH1F* sig_rms;
00403 
00404   TH2F* const_eta_phi;
00405 
00406   TH1F* const_eta[netamx];
00407   TH1F* stat_eta[netamx];
00408   TH1F* statmn_eta[netamx];  
00409   TH1F* peak_eta[netamx]; 
00410   
00411   TH1F* const_hpdrm[ringmx];
00412   //  TH1F* stat_hpdrm[ringmx];
00413   //  TH1F* statmn_hpdrm[ringmx];  
00414   TH1F* peak_hpdrm[ringmx]; 
00415 
00416   TH1F* mean_eta_ave;
00417   TH1F* mean_phi_ave;
00418   TH1F* mean_phi_hst;
00419 
00420   TH2F* sig_effi[neffip];
00421   TH2F* mean_energy;
00422 
00423   double fitprm[nsgpr][netamx];
00424 
00425 
00426   TProfile* sigvsevt[15][ncut];
00427 
00428 
00429   //  int   irun, ievt, itrg1, itrg2, isect, nrecht, nfound, nlost, ndof, nmuon;
00430   int   irun, ievt, itrg1, itrg2, isect, ndof, nmuon;
00431   float trkdr, trkdz, trkvx, trkvy, trkvz, trkmm, trkth, trkph, chisq, therr, pherr, hodx, hody, 
00432     hoang, htime, hosig[9], hocorsig[18], hocro, hbhesig[9], caloen[3];
00433 
00434   int Nevents; 
00435   int nbn;
00436   float alow;
00437   float ahigh;
00438   float binwid;
00439   int irunold;
00440 
00441   edm::InputTag hoCalibVariableCollectionTag;
00442 
00443       // ----------member data ---------------------------
00444 
00445 };
00446 
00447 const int HOCalibAnalyzer::ringmx;
00448 const int HOCalibAnalyzer::sectmx;
00449 const int HOCalibAnalyzer::routmx;
00450 const int HOCalibAnalyzer::rout12mx;
00451 const int HOCalibAnalyzer::neffip;
00452 
00453 //
00454 // constants, enums and typedefs
00455 //
00456 
00457 //
00458 // static data member definitions
00459 //
00460 
00461 //
00462 // constructors and destructor
00463 //
00464 HOCalibAnalyzer::HOCalibAnalyzer(const edm::ParameterSet& iConfig) :
00465   hoCalibVariableCollectionTag(iConfig.getParameter<edm::InputTag>("hoCalibVariableCollectionTag"))
00466   // It is very likely you want the following in your configuration
00467   // hoCalibVariableCollectionTag = cms.InputTag('hoCalibProducer', 'HOCalibVariableCollection')
00468 {
00469    //now do what ever initialization is needed
00470   ipass = 0;
00471   Nevents = 0;
00472   
00473   theRootFileName = iConfig.getUntrackedParameter<string>("RootFileName", "test.root");
00474   theoutputtxtFile = iConfig.getUntrackedParameter<string>("txtFileName", "test.txt");
00475   theoutputpsFile = iConfig.getUntrackedParameter<string>("psFileName", "test.ps");
00476 
00477   m_hbinfo = iConfig.getUntrackedParameter<bool>("hbinfo", false);
00478   m_hbtime = iConfig.getUntrackedParameter<bool>("hbtime", false);
00479   m_hotime = iConfig.getUntrackedParameter<bool>("hotime", false);
00480   m_correl = iConfig.getUntrackedParameter<bool>("correl", false);
00481   m_checkmap = iConfig.getUntrackedParameter<bool>("checkmap", false);
00482   m_combined = iConfig.getUntrackedParameter<bool>("combined", false);
00483   m_constant = iConfig.getUntrackedParameter<bool>("get_constant", false);
00484   m_figure = iConfig.getUntrackedParameter<bool>("get_figure", true);
00485   m_digiInput = iConfig.getUntrackedParameter<bool>("digiInput", true);
00486   m_histfit = iConfig.getUntrackedParameter<bool>("histFit", true);
00487   m_pedsuppr = iConfig.getUntrackedParameter<bool>("pedSuppr", true);
00488   m_cosmic = iConfig.getUntrackedParameter<bool>("cosmic", true);
00489   m_sigma = iConfig.getUntrackedParameter<double>("sigma", 1.0);
00490   
00491   edm::Service<TFileService> fs;
00492 
00493   theFile = new TFile(theRootFileName.c_str(), "RECREATE");
00494   theFile->cd();
00495   
00496   T1 = new TTree("T1", "DT+CSC+HO");
00497   
00498   T1->Branch("irun",&irun,"irun/I");  
00499   T1->Branch("ievt",&ievt,"ievt/I");  
00500   T1->Branch("itrg1",&itrg1,"itrg1/I");  
00501   T1->Branch("itrg2",&itrg2,"itrg2/I");  
00502 
00503   T1->Branch("isect",&isect,"isect/I"); 
00504   //  T1->Branch("nrecht",&nrecht,"nrecht/I"); 
00505   T1->Branch("ndof",&ndof,"ndof/I"); 
00506   T1->Branch("nmuon",&nmuon,"nmuon/I"); 
00507 
00508   T1->Branch("trkdr",&trkdr,"trkdr/F");
00509   T1->Branch("trkdz",&trkdz,"trkdz/F");
00510  
00511   T1->Branch("trkvx",&trkvx,"trkvx/F");
00512   T1->Branch("trkvy",&trkvy,"trkvy/F");
00513   T1->Branch("trkvz",&trkvz,"trkvz/F");
00514   T1->Branch("trkmm",&trkmm,"trkmm/F");
00515   T1->Branch("trkth",&trkth,"trkth/F");
00516   T1->Branch("trkph",&trkph,"trkph/F");
00517 
00518   T1->Branch("chisq",&chisq,"chisq/F");
00519   T1->Branch("therr",&therr,"therr/F");
00520   T1->Branch("pherr",&pherr,"pherr/F");
00521   T1->Branch("hodx",&hodx,"hodx/F");
00522   T1->Branch("hody",&hody,"hody/F");
00523   T1->Branch("hoang",&hoang,"hoang/F");
00524 
00525   T1->Branch("htime",&htime,"htime/F");  
00526   T1->Branch("hosig",hosig,"hosig[9]/F");
00527   T1->Branch("hocro",&hocro,"hocro/F");
00528   T1->Branch("hocorsig",hocorsig,"hocorsig[18]/F");
00529   T1->Branch("caloen",caloen,"caloen[3]/F");
00530 
00531   if (m_hbinfo) { // #ifdef HBINFO
00532     T1->Branch("hbhesig",hbhesig,"hbhesig[9]/F");
00533   } //m_hbinfo #endif
00534 
00535   muonnm = fs->make<TH1F>("muonnm", "No of muon", 10, -0.5, 9.5);
00536   muonmm = fs->make<TH1F>("muonmm", "P_{mu}", 200, -100., 100.);
00537   muonth = fs->make<TH1F>("muonth", "{Theta}_{mu}", 180, 0., 180.);
00538   muonph = fs->make<TH1F>("muonph", "{Phi}_{mu}", 180, -180., 180.);
00539   muonch = fs->make<TH1F>("muonch", "{chi^2}/ndf", 100, 0., 1000.);
00540 
00541   sel_muonnm = fs->make<TH1F>("sel_muonnm", "No of muon(sel)", 10, -0.5, 9.5);
00542   sel_muonmm = fs->make<TH1F>("sel_muonmm", "P_{mu}(sel)", 200, -100., 100.);
00543   sel_muonth = fs->make<TH1F>("sel_muonth", "{Theta}_{mu}(sel)", 180, 0., 180.);
00544   sel_muonph = fs->make<TH1F>("sel_muonph", "{Phi}_{mu}(sel)", 180, -180., 180.);
00545   sel_muonch = fs->make<TH1F>("sel_muonch", "{chi^2}/ndf(sel)", 100, 0., 1000.);
00546 
00547 
00548   int nbin = 50; //40;// 45; //50; //55; //60; //55; //45; //40; //50;
00549   alow = -2.0;// -1.85; //-1.90; // -1.95; // -2.0;
00550   ahigh = 8.0;// 8.15; // 8.10; //  8.05; //  8.0;
00551 
00552   if (m_digiInput) {alow = -10.0; ahigh = 40.;}
00553 
00554   float tmpwid = (ahigh-alow)/nbin;
00555   nbn = int(-alow/tmpwid)+1;
00556   if (nbn <0) nbn = 0;
00557   if (nbn>nbin) nbn = nbin;
00558 
00559   char name[200];
00560   char title[200];
00561 
00562   cout <<"nbin "<< nbin<<" "<<alow<<" "<<ahigh<<" "<<tmpwid<<" "<<nbn<<endl;
00563 
00564   for (int i=0; i<15; i++) {
00565     
00566     sprintf(title, "sigvsndof_ring%i", i+1); 
00567     sigvsevt[i][0] = fs->make<TProfile>(title, title, 50, 0., 50.,-9., 20.);
00568     
00569     sprintf(title, "sigvschisq_ring%i", i+1); 
00570     sigvsevt[i][1] = fs->make<TProfile>(title, title, 50, 0., 30.,-9., 20.);
00571     
00572     sprintf(title, "sigvsth_ring%i", i+1); 
00573     sigvsevt[i][2] = fs->make<TProfile>(title, title, 50, .7, 2.4,-9., 20.);
00574     
00575     sprintf(title, "sigvsph_ring%i", i+1); 
00576     sigvsevt[i][3] = fs->make<TProfile>(title, title, 50, -2.4, -0.7,-9., 20.);
00577     
00578     sprintf(title, "sigvstherr_ring%i", i+1); 
00579     sigvsevt[i][4] = fs->make<TProfile>(title, title, 50, 0., 0.2,-9., 20.);
00580       
00581     sprintf(title, "sigvspherr_ring%i", i+1); 
00582     sigvsevt[i][5] = fs->make<TProfile>(title, title, 50, 0., 0.2,-9., 20.);
00583     
00584     sprintf(title, "sigvsdircos_ring%i", i+1); 
00585     sigvsevt[i][6] = fs->make<TProfile>(title, title, 50, 0.5, 1.,-9., 20.);
00586     
00587     sprintf(title, "sigvstrkmm_ring%i", i+1); 
00588     sigvsevt[i][7] = fs->make<TProfile>(title, title, 50, 0., 50.,-9., 20.);
00589     
00590     sprintf(title, "sigvsnmuon_ring%i", i+1); 
00591     sigvsevt[i][8] = fs->make<TProfile>(title, title, 5, 0.5, 5.5,-9., 20.);
00592     
00593     sprintf(title, "sigvserr_ring%i", i+1); 
00594     sigvsevt[i][9] = fs->make<TProfile>(title, title, 50, 0., .3, -9., 20.);
00595     
00596     sprintf(title, "sigvsaccx_ring%i", i+1); 
00597     sigvsevt[i][10] = fs->make<TProfile>(title, title, 100, -25., 25., -9., 20.);    
00598 
00599     sprintf(title, "sigvsaccy_ring%i", i+1); 
00600     sigvsevt[i][11] = fs->make<TProfile>(title, title, 100, -25., 25., -9., 20.);    
00601     
00602     sprintf(title, "sigvscalo_ring%i", i+1); 
00603     sigvsevt[i][12] = fs->make<TProfile>(title, title, 100, 0., 15., -9., 20.);    
00604   }
00605 
00606   for (int j=0; j<netamx; j++) {
00607     int ieta = (j<15) ? j+1 : 14-j;
00608     for (int i=0;i<nphimx+1;i++) {
00609       if (i==nphimx) {
00610         sprintf(title, "sig_eta%i_allphi", ieta);
00611       } else {
00612         sprintf(title, "sig_eta%i_phi%i", ieta,i+1);
00613       }
00614       sigrsg[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh); 
00615       if (i==nphimx) {
00616         sprintf(title, "ped_eta%i_allphi", ieta);
00617       } else {
00618         sprintf(title, "ped_eta%i_phi%i", ieta,i+1);
00619       }
00620       crossg[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh); 
00621     }
00622 
00623     for (int i=0;i<nphimx;i++) {
00624       if (m_hotime) { //#ifdef HOTIME
00625         sprintf(title, "hotime_eta%i_phi%i", (j<=14) ? j+1 : 14-j, i+1);
00626         hotime[j][i] = fs->make<TProfile>(title, title, 10, -0.5, 9.5, -1.0, 30.0);
00627         
00628         sprintf(title, "hopedtime_eta%i_phi%i", (j<=14) ? j+1 : 14-j, i+1);
00629         hopedtime[j][i] = fs->make<TProfile>(title, title, 10, -0.5, 9.5, -1.0, 30.0);
00630         
00631       } //m_hotime #endif
00632       if (m_hbtime) { //#ifdef HBTIME
00633         sprintf(title, "hbtime_eta%i_phi%i", (j<=15) ? j+1 : 15-j, i+1);
00634         hbtime[j][i] = fs->make<TProfile>(title, title, 10, -0.5, 9.5, -1.0, 30.0);
00635       } //m_hbtime #endif
00636 
00637       if (m_correl) { //#ifdef CORREL    
00638         sprintf(title, "corrsg_eta%i_phi%i_leftbottom", ieta,i+1);
00639         corrsglb[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);   
00640         
00641         sprintf(title, "corrsg_eta%i_phi%i_rightbottom", ieta,i+1);
00642         corrsgrb[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);   
00643         
00644         sprintf(title, "corrsg_eta%i_phi%i_leftup", ieta,i+1);
00645         corrsglu[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);   
00646         
00647         sprintf(title, "corrsg_eta%i_phi%i_rightup", ieta,i+1);
00648         corrsgru[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);   
00649         
00650         sprintf(title, "corrsg_eta%i_phi%i_all", ieta,i+1);
00651         corrsgall[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);   
00652         
00653         sprintf(title, "corrsg_eta%i_phi%i_left", ieta,i+1);
00654         corrsgl[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);   
00655         
00656         sprintf(title, "corrsg_eta%i_phi%i_right", ieta,i+1);
00657         corrsgr[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
00658       } //m_correl #endif
00659       if (m_checkmap) {// #ifdef CHECKMAP    
00660         sprintf(title, "corrsg_eta%i_phi%i_centrl", ieta,i+1);
00661         corrsgc[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);   
00662       } //m_checkmap #endif
00663     }
00664   }
00665 
00666   mnsigrsg = fs->make<TH1F>("mnsigrsg","mnsigrsg", netamx*nphimx+ringmx*routmx, -0.5, netamx*nphimx+ringmx*routmx -0.5);
00667   rmssigrsg = fs->make<TH1F>("rmssigrsg","rmssigrsg", netamx*nphimx+ringmx*routmx, -0.5, netamx*nphimx+ringmx*routmx -0.5);
00668   nevsigrsg = fs->make<TH1F>("nevsigrsg","nevsigrsg", netamx*nphimx+ringmx*routmx, -0.5, netamx*nphimx+ringmx*routmx -0.5);  
00669   
00670   mncrossg = fs->make<TH1F>("mncrossg","mncrossg", netamx*nphimx+ringmx*routmx, -0.5, netamx*nphimx+ringmx*routmx -0.5);
00671   rmscrossg = fs->make<TH1F>("rmscrossg","rmscrossg", netamx*nphimx+ringmx*routmx, -0.5, netamx*nphimx+ringmx*routmx -0.5);
00672   nevcrossg = fs->make<TH1F>("nevcrossg","nevcrossg", netamx*nphimx+ringmx*routmx, -0.5, netamx*nphimx+ringmx*routmx -0.5);  
00673 
00674   
00675   for (int i=0; i<neffip; i++) {
00676     if (i==0) {
00677       sprintf(title, "Total projected muon in tower"); 
00678       sprintf(name, "total_evt"); 
00679     } else {
00680       sprintf(title, "Efficiency with sig >%i #sigma", i); 
00681       sprintf(name, "Effi_with_gt%i_sig", i); 
00682     }
00683     sig_effi[i] = fs->make<TH2F>(name, title, netamx+1, -netamx/2-0.5, netamx/2+0.5, nphimx, 0.5, nphimx+0.5);
00684   }
00685 
00686   sprintf(title, "Mean Energy of all towers"); 
00687   sprintf(name, "mean_energy"); 
00688   mean_energy = fs->make<TH2F>(name, title, netamx+1, -netamx/2-0.5, netamx/2+0.5, nphimx, 0.5, nphimx+0.5);
00689 
00690   if (m_correl) { //#ifdef CORREL    
00691     mncorrsglb = fs->make<TH1F>("mncorrsglb","mncorrsglb", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
00692     rmscorrsglb = fs->make<TH1F>("rmscorrsglb","rmscorrsglb", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
00693     nevcorrsglb = fs->make<TH1F>("nevcorrsglb","nevcorrsglb", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
00694     
00695     mncorrsgrb = fs->make<TH1F>("mncorrsgrb","mncorrsgrb", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
00696     rmscorrsgrb = fs->make<TH1F>("rmscorrsgrb","rmscorrsgrb", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
00697     nevcorrsgrb = fs->make<TH1F>("nevcorrsgrb","nevcorrsgrb", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
00698 
00699     mncorrsglu = fs->make<TH1F>("mncorrsglu","mncorrsglu", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
00700     rmscorrsglu = fs->make<TH1F>("rmscorrsglu","rmscorrsglu", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
00701     nevcorrsglu = fs->make<TH1F>("nevcorrsglu","nevcorrsglu", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
00702     
00703     mncorrsgru = fs->make<TH1F>("mncorrsgru","mncorrsgru", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
00704     rmscorrsgru = fs->make<TH1F>("rmscorrsgru","rmscorrsgru", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
00705     nevcorrsgru = fs->make<TH1F>("nevcorrsgru","nevcorrsgru", netamx*nphimx+60, -0.5, netamx*nphimx+59.5); 
00706     
00707     mncorrsgall = fs->make<TH1F>("mncorrsgall","mncorrsgall", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
00708     rmscorrsgall = fs->make<TH1F>("rmscorrsgall","rmscorrsgall", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
00709     nevcorrsgall = fs->make<TH1F>("nevcorrsgall","nevcorrsgall", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
00710     
00711     mncorrsgl = fs->make<TH1F>("mncorrsgl","mncorrsgl", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
00712     rmscorrsgl = fs->make<TH1F>("rmscorrsgl","rmscorrsgl", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
00713     nevcorrsgl = fs->make<TH1F>("nevcorrsgl","nevcorrsgl", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);  
00714     
00715     mncorrsgr = fs->make<TH1F>("mncorrsgr","mncorrsgr", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
00716     rmscorrsgr = fs->make<TH1F>("rmscorrsgr","rmscorrsgr", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
00717     nevcorrsgr = fs->make<TH1F>("nevcorrsgr","nevcorrsgr", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);  
00718   } //m_correl #endif
00719   
00720   if (m_checkmap) { //#ifdef CHECKMAP    
00721     mncorrsgc = fs->make<TH1F>("mncorrsgc","mncorrsgc", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
00722     rmscorrsgc = fs->make<TH1F>("rmscorrsgc","rmscorrsgc", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);
00723     nevcorrsgc = fs->make<TH1F>("nevcorrsgc","nevcorrsgc", netamx*nphimx+60, -0.5, netamx*nphimx+59.5);  
00724   } //m_checkmap #endif
00725   
00726   if (m_combined) { //#ifdef COMBINED
00727     for (int j=0; j<ringmx; j++) {
00728 
00729       for (int i=0;i<routmx+1;i++) {
00730         if (j!=2 && i>rout12mx) continue;
00731         int phmn = 3*i-1;
00732         int phmx = 3*i+1;
00733         if (j==2) {phmn = 2*i-1; phmx=2*i;}
00734         if (phmn <=0) phmn = nphimx+phmn;
00735         if (phmx <=0) phmx = nphimx+phmx;
00736         
00737         if ((j==2 && i==routmx) || (j!=2 && i==rout12mx)) {
00738           sprintf(title, "sig_ring%i_allrm", j-2);
00739           sprintf(name, "sig_ring%i_allrm", j-2);
00740         } else {
00741           sprintf(title, "sig_ring%i_phi%i-%i", j-2,phmn,phmx);
00742           sprintf(name, "sig_ring%i_rout%i", j-2,i+1);
00743         }
00744         com_sigrsg[j][i] = fs->make<TH1F>(name, title, nbin, alow, ahigh);
00745         if ((j==2 && i==routmx) || (j!=2 && i==rout12mx)) {
00746           sprintf(title, "ped_ring%i_allrm", j-2);
00747           sprintf(name, "ped_ring%i_allrm", j-2);
00748         } else {
00749           sprintf(title, "ped_ring%i_phi%i-%i", j-2,phmn, phmx);
00750           sprintf(name, "ped_ring%i_rout%i", j-2,i+1);
00751         }
00752         com_crossg[j][i] = fs->make<TH1F>(name, title, nbin, alow, ahigh);   
00753       }
00754 
00755       for (int i=0;i<sectmx;i++) {
00756         if (m_hotime) { //#ifdef HOTIME
00757           sprintf(title, "com_hotime_ring%i_sect%i", j-2, i+1);
00758           com_hotime[j][i] = fs->make<TProfile>(title, title, 10, -0.5, 9.5, -1.0, 30.0);
00759           
00760           sprintf(title, "com_hopedtime_ring%i_sect%i", j-2, i+1);
00761           com_hopedtime[j][i] = fs->make<TProfile>(title, title, 10, -0.5, 9.5, -1.0, 30.0);
00762         } //m_hotime #endif      
00763         if (m_hbtime){ //#ifdef HBTIME
00764           sprintf(title, "_com_hbtime_ring%i_serrct%i", j-2, i+1);
00765           com_hbtime[j][i] = fs->make<TProfile>(title, title, 10, -0.5, 9.5, -1.0, 30.0);
00766         } //m_hbtime #endif
00767         
00768         if (m_correl) { //#ifdef CORREL    
00769           sprintf(title, "com_corrsg_ring%i_sect%i_leftbottom", j-2,i+1);
00770           com_corrsglb[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);   
00771           
00772           sprintf(title, "com_corrsg_ring%i_sect%i_rightbottom", j-2,i+1);
00773           com_corrsgrb[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);   
00774           
00775           sprintf(title, "com_corrsg_ring%i_sect%i_leftup", j-2,i+1);
00776           com_corrsglu[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);   
00777           
00778           sprintf(title, "com_corrsg_ring%i_sect%i_rightup", j-2,i+1);
00779           com_corrsgru[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);   
00780           
00781           sprintf(title, "com_corrsg_ring%i_sect%i_all", j-2,i+1);
00782           com_corrsgall[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);   
00783           
00784           sprintf(title, "com_corrsg_ring%i_sect%i_left", j-2,i+1);
00785           com_corrsgl[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);   
00786           
00787           sprintf(title, "com_corrsg_ring%i_sect%i_right", j-2,i+1);
00788           com_corrsgr[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);   
00789         } //m_correl #endif
00790         
00791         if (m_checkmap) { // #ifdef CHECKMAP    
00792           sprintf(title, "com_corrsg_ring%i_sect%i_centrl", j-2,i+1);
00793           com_corrsgc[j][i] = fs->make<TH1F>(title, title, nbin, alow, ahigh);   
00794         } //m_checkmap #endif
00795       }
00796     }
00797   } //m_combined #endif  
00798 
00799   for (int i=-1; i<=1; i++) {
00800     for (int j=-1; j<=1; j++) {    
00801       int k = 3*(i+1)+j+1;
00802       
00803       sprintf(title, "hosct2p_eta%i_phi%i", i, j);
00804       ho_sig2p[k] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
00805       
00806       sprintf(title, "hosct1p_eta%i_phi%i", i, j);
00807       ho_sig1p[k] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
00808 
00809       sprintf(title, "hosct00_eta%i_phi%i", i, j);
00810       ho_sig00[k] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
00811 
00812       sprintf(title, "hosct1m_eta%i_phi%i", i, j);
00813       ho_sig1m[k] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
00814 
00815       sprintf(title, "hosct2m_eta%i_phi%i", i, j);
00816       ho_sig2m[k] = fs->make<TH1F>(title, title, nbin, alow, ahigh);
00817 
00818       if (m_hbinfo) { // #ifdef HBINFO
00819         sprintf(title, "hbhesig_eta%i_phi%i", i, j);
00820         hbhe_sig[k] = fs->make<TH1F>(title, title, 51, -10.5, 40.5);
00821       } //m_hbinfo #endif
00822     }
00823   }
00824 
00825   if (m_constant) {
00826     ped_evt = fs->make<TH1F>("ped_evt", "ped_evt", netamx*nphimx, -0.5, netamx*nphimx-0.5);
00827     ped_mean = fs->make<TH1F>("ped_mean", "ped_mean", netamx*nphimx, -0.5, netamx*nphimx-0.5); 
00828     ped_width = fs->make<TH1F>("ped_width", "ped_width", netamx*nphimx, -0.5, netamx*nphimx-0.5); 
00829     
00830     fit_chi = fs->make<TH1F>("fit_chi", "fit_chi", netamx*nphimx, -0.5, netamx*nphimx-0.5);
00831     sig_evt = fs->make<TH1F>("sig_evt", "sig_evt", netamx*nphimx, -0.5, netamx*nphimx-0.5);
00832     fit_sigevt = fs->make<TH1F>("fit_sigevt", "fit_sigevt", netamx*nphimx, -0.5, netamx*nphimx-0.5);
00833     fit_bkgevt = fs->make<TH1F>("fit_bkgevt", "fit_bkgevt", netamx*nphimx, -0.5, netamx*nphimx-0.5);
00834     sig_mean = fs->make<TH1F>("sig_mean", "sig_mean", netamx*nphimx, -0.5, netamx*nphimx-0.5);       
00835     sig_diff = fs->make<TH1F>("sig_diff", "sig_diff", netamx*nphimx, -0.5, netamx*nphimx-0.5);       
00836     sig_width = fs->make<TH1F>("sig_width", "sig_width", netamx*nphimx, -0.5, netamx*nphimx-0.5);
00837     sig_sigma = fs->make<TH1F>("sig_sigma", "sig_sigma", netamx*nphimx, -0.5, netamx*nphimx-0.5);
00838     sig_meanerr = fs->make<TH1F>("sig_meanerr", "sig_meanerr", netamx*nphimx, -0.5, netamx*nphimx-0.5); 
00839     sig_meanerrp = fs->make<TH1F>("sig_meanerrp", "sig_meanerrp", netamx*nphimx, -0.5, netamx*nphimx-0.5); 
00840     sig_signf = fs->make<TH1F>("sig_signf", "sig_signf", netamx*nphimx, -0.5, netamx*nphimx-0.5); 
00841 
00842     ped_statmean = fs->make<TH1F>("ped_statmean", "ped_statmean", netamx*nphimx, -0.5, netamx*nphimx-0.5);
00843     sig_statmean = fs->make<TH1F>("sig_statmean", "sig_statmean", netamx*nphimx, -0.5, netamx*nphimx-0.5);
00844     ped_rms = fs->make<TH1F>("ped_rms", "ped_rms", netamx*nphimx, -0.5, netamx*nphimx-0.5);
00845     sig_rms = fs->make<TH1F>("sig_rms", "sig_rms", netamx*nphimx, -0.5, netamx*nphimx-0.5);
00846 
00847     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);
00848     
00849     for (int ij=0; ij<netamx; ij++) {
00850       int ieta = (ij<15) ? ij+1 : 14-ij;
00851       sprintf(title, "Cont_Eta_%i", ieta);
00852       const_eta[ij] = fs->make<TH1F>(title, title, nphimx, 0.5, nphimx+0.5);
00853       
00854       sprintf(title, "Peak_Eta_%i", ieta);
00855       peak_eta[ij] =  fs->make<TH1F>(title, title, nphimx, 0.5, nphimx+0.5);
00856     }
00857 
00858     for (int ij=0; ij<ringmx; ij++) {
00859       int iring = ij-2;
00860       int iread = (ij==2) ? routmx : rout12mx;
00861       sprintf(title, "Cont_hpdrm_%i", iring);
00862       const_hpdrm[ij] = fs->make<TH1F>(title, title, iread, 0.5, iread+0.5);
00863       
00864       sprintf(title, "Peak_hpdrm_%i", iring);
00865       peak_hpdrm[ij] =  fs->make<TH1F>(title, title, iread, 0.5, iread+0.5);
00866     }
00867 
00868     mean_phi_hst = fs->make<TH1F>("mean_phi_hst", "mean_phi_hst", netamx+1, -(netamx+1)/2., (netamx+1)/2.);
00869     mean_phi_ave = fs->make<TH1F>("mean_phi_ave", "mean_phi_ave", netamx+1, -(netamx+1)/2., (netamx+1)/2.);
00870 
00871     mean_eta_ave = fs->make<TH1F>("mean_eta_ave", "mean_eta_ave", nphimx, 0.5, nphimx+0.5);
00872 
00873   } //m_constant
00874  
00875   for (int ij=0; ij<netamx; ij++) {
00876     int ieta = (ij<15) ? ij+1 : 14-ij;
00877     
00878     sprintf(title, "Stat_Eta_%i", ieta);
00879     stat_eta[ij] =  fs->make<TH1F>(title, title, nphimx, 0.5, nphimx+0.5);
00880     
00881     sprintf(title, "#mu(stat)_Eta_%i", ieta);
00882     statmn_eta[ij] =  fs->make<TH1F>(title, title, nphimx, 0.5, nphimx+0.5);
00883   }
00884 
00885   for (int j=0; j<netamx; j++) {
00886     for (int i=0; i<nphimx; i++) {             
00887       invang[j][i] = 0.0;
00888     }
00889   }
00890   for (int j=0; j<ringmx; j++) {
00891     for (int i=0;i<routmx+1;i++) {             
00892       com_invang[j][i] = 0.0;
00893     }
00894   }
00895 
00896 }
00897 
00898 
00899 HOCalibAnalyzer::~HOCalibAnalyzer()
00900 {
00901  
00902    // do anything here that needs to be done at desctruction time
00903    // (e.g. close files, deallocate resources etc.)
00904 
00905   theFile->cd();
00906   theFile->Write();
00907   theFile->Close();
00908   cout <<" Selected events # is "<<ipass<<endl;
00909 }
00910 
00911 
00912 //
00913 // member functions
00914 //
00915 
00916 // ------------ method called to for each event  ------------
00917 void
00918 HOCalibAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00919 {
00920 
00921   // calcualte these once (and avoid the pow(int,int) ambiguities for c++)
00922   int mypow_2_0 = 0; // 2^0
00923   int mypow_2_1 = 2; // 2^1
00924   int mypow_2_2 = 4; // 2^2
00925 
00926   int mypow_2_3 = 8; // 2^3
00927   int mypow_2_4 = 16; // 2^4
00928   int mypow_2_5 = 32; // 2^5
00929   int mypow_2_6 = 64; // 2^6
00930   int mypow_2_7 = 128; // 2^7
00931   int mypow_2_8 = 256; // 2^8
00932   int mypow_2_9 = 512; // 2^9
00933   int mypow_2_10 = 1024; // 2^10
00934   int mypow_2_11 = 2048; // 2^11
00935   int mypow_2_12 = 4096; // 2^12
00936 
00937   /*
00938   //FIXGM Put this is initialiser
00939   int mapx1[6][3]={{1,4,8}, {12,7,3}, {5,9,13}, {11,6,2}, {16,15,14}, {19,18,17}}; 
00940   //    int etamap1[21]={-1, 0,3,1, 0,2,3, 1,0,2, -1, 3,1,2, 4,4,4, 5,5,5, -1};
00941   //  int phimap1[21]={-1, 0,2,2, 1,0,1, 1,2,1, -1, 0,0,2, 2,1,0, 2,1,0,-1};
00942 
00943     int mapx2[6][3]={{1,4,8}, {12,7,3}, {5,9,13}, {11,6,2}, {16,15,14}, {-1,-1,-1}};
00944   //  int etamap2[21]={-1, 0,3,1, 0,2,3, 1,0,2, -1, 3,1,2, 4,4,4, -1,-1,-1, -1};
00945   //  int phimap2[21]={-1, 0,2,2, 1,0,1, 1,2,1, -1, 0,0,2, 2,1,0,  2, 1, 0, -1};
00946 
00947     int mapx0p[9][2]={{3,1}, {7,4}, {6,5},  {12,8}, {0,0}, {11,9}, {16,13}, {15,14}, {19,17}};
00948     int mapx0m[9][2]={{17,19}, {14,15}, {13,16}, {9,11}, {0,0}, {8,12}, {5,6}, {4,7}, {1,3}};
00949 
00950   //  int etamap0p[21]={-1, 0,-1,0, 1,2,2, 1,3,5, -1, 5,3,6, 7,7,6, 8,-1,8, -1};
00951   //  int phimap0p[21]={-1, 1,-1,0, 1,1,0, 0,1,1, -1, 0,0,1, 1,0,0, 1,-1,0, -1};
00952 
00953   //  int etamap0m[21]={-1, 8,-1,8, 7,6,6, 7,5,3, -1, 3,5,2, 1,1,2, 0,-1,0, -1};
00954   //  int phimap0m[21]={-1, 0,-1,1, 0,0,1, 1,0,0, -1, 1,1,0, 0,1,1, 0,-1,1, -1};
00955 
00956     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
00957                        {-1, 0,3,1, 0,2,3, 1,0,2, -1, 3,1,2, 4,4,4, 5,5,5, -1},    //etamap1 
00958                        {-1, 0,-1,0, 1,2,2, 1,3,5, -1, 5,3,6, 7,7,6, 8,-1,8, -1},  //etamap0p
00959                        {-1, 8,-1,8, 7,6,6, 7,5,3, -1, 3,5,2, 1,1,2, 0,-1,0, -1}}; //etamap0m
00960 
00961     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
00962                         {-1, 0,2,2, 1,0,1, 1,2,1, -1, 0,0,2, 2,1,0, 2,1,0, -1},    //phimap1
00963                         {-1, 1,-1,0, 1,1,0, 0,1,1, -1, 0,0,1, 1,0,0, 1,-1,0, -1},  //phimap0p
00964                         {-1, 0,-1,1, 0,0,1, 1,0,0, -1, 1,1,0, 0,1,1, 0,-1,1, -1}};  //phimap0m
00965   //swapped phi map for R0+/R0- (15/03/07)  
00966   for (int i=0; i<4; i++) {
00967     for (int j=0; j<21; j++) {
00968       cout <<"ieta "<<i<<" "<<j<<" "<<etamap[i][j]<<endl;
00969     }
00970   }
00971 
00972   // Character convention for R+/-1/2
00973   //      int npixleft[21] = {-1, F, Q,-1, M, D, J,-1, T,-1, C,-1, R, P, H,-1, N, G};
00974   //      int npixrigh[21] = { Q, S,-1, D, J, L,-1, K,-1, E,-1, P, H, B,-1, G, A,-1};
00975   
00976   //      int npixlb1[21]={-1,-1,-1,-1, F, Q, S,-1, M, J, L, T, K,-1, C, R, P, H};
00977   //      int npixrb1[21]={-1,-1,-1, F, Q, S,-1, M, D, L,-1, K,-1, C, E, P, H, B};
00978   //      int npixlu1[21]={ M, D, J, T, K,-1, C,-1, R, H, B,-1, N, G, A,-1,-1,-1};
00979   //      int npixru1[21]={ D, J, L, K,-1, C, E, R, P, B,-1, N, G, A,-1,-1,-1,-1};
00980 
00981   int npixleft[21]={0, 0, 1, 2, 0, 4, 5, 6, 0, 8, 0, 0,11, 0,13,14,15, 0,17,18,0};
00982   int npixrigh[21]={0, 2, 3, 0, 5, 6, 7, 0, 9, 0, 0,12, 0,14,15,16, 0,18,19, 0,0};
00983   int npixlebt[21]={0, 0, 0, 0, 0, 1, 2, 3, 0, 4, 0, 6, 7, 8, 9, 0,11,13,14,15,0};
00984   int npixribt[21]={0, 0, 0, 0, 1, 2, 3, 0, 4, 5, 0, 7, 0, 9, 0,11,12,14,15,16,0};
00985   int npixleup[21]={0, 4, 5, 6, 8, 9, 0,11, 0,13, 0,15,16, 0,17,18,19, 0, 0, 0,0};
00986   int npixriup[21]={0, 5, 6, 7, 9, 0,11,12,13,14, 0,16, 0,17,18,19, 0, 0, 0, 0,0};
00987   */
00988   
00989   int iaxxx = 0;
00990   int ibxxx = 0;
00991   
00992   Nevents++;
00993 
00994   using namespace edm;
00995 
00996   float pival = acos(-1.);
00997   irunold = irun = iEvent.id().run();
00998   ievt = iEvent.id().event();
00999 
01000   //  cout <<"Nevents "<<Nevents<<endl;
01001   
01002   edm::Handle<HOCalibVariableCollection>HOCalib;
01003   bool isCosMu = true;
01004   try {
01005     iEvent.getByLabel(hoCalibVariableCollectionTag, HOCalib); 
01006     //    iEvent.getByLabel("hoCalibProducer","HOCalibVariableCollection",HOCalib);
01007 
01008   } catch ( cms::Exception &iEvent ) { isCosMu = false; } 
01009 
01010   if (isCosMu && (*HOCalib).size() >0 ) { 
01011     nmuon = (*HOCalib).size();
01012     if (Nevents%5000==1)   cout <<"nmuon event # "<<Nevents<<" Run # "<<iEvent.id().run()<<" Evt # "<<iEvent.id().event()<<" "<<nmuon<<endl;
01013     
01014     for (HOCalibVariableCollection::const_iterator hoC=(*HOCalib).begin(); hoC!=(*HOCalib).end(); hoC++){
01015       itrg1 = (*hoC).trig1;
01016       itrg2 = (*hoC).trig2;
01017 
01018       trkdr = (*hoC).trkdr;
01019       trkdz = (*hoC).trkdz;
01020 
01021       trkvx = (*hoC).trkvx;
01022       trkvy = (*hoC).trkvy;
01023       trkvz = (*hoC).trkvz;
01024       
01025       trkmm = (*hoC).trkmm;
01026       trkth = (*hoC).trkth;
01027       trkph = (*hoC).trkph;
01028 
01029       ndof   = (int)(*hoC).ndof;
01030       //      nrecht = (int)(*hoC).nrecht;
01031       chisq = (*hoC).chisq;
01032 
01033       therr = (*hoC).therr;
01034       pherr = (*hoC).pherr;
01035       trkph = (*hoC).trkph;
01036 
01037       isect = (*hoC).isect;
01038       hodx  = (*hoC).hodx;
01039       hody  = (*hoC).hody;
01040       hoang = (*hoC).hoang;
01041       htime = (*hoC).htime;
01042       for (int i=0; i<9; i++) { hosig[i] = (*hoC).hosig[i];} //cout<<"hosig "<<i<<" "<<hosig[i]<<endl;}
01043       for (int i=0; i<18; i++) { hocorsig[i] = (*hoC).hocorsig[i];} // cout<<"hocorsig "<<i<<" "<<hocorsig[i]<<endl;}    
01044       hocro = (*hoC).hocro;
01045       for (int i=0; i<3; i++) { caloen[i] = (*hoC).caloen[i];}
01046 
01047       if (m_hbinfo) { for (int i=0; i<9; i++) { hbhesig[i] = (*hoC).hbhesig[i];}} // cout<<"hbhesig "<<i<<" "<<hbhesig[i]<<endl;}}
01048 
01049 
01050       int ipsall=0;
01051       int ips0=0; 
01052       int ips1=0; 
01053       int ips2=0; 
01054       int ips3=0; 
01055       int ips4=0; 
01056       int ips5=0; 
01057       int ips6=0; 
01058       int ips7=0; 
01059       int ips8=0; 
01060       int ips9=0;
01061       int ips10 =0;
01062       int ips11 =0;
01063       int ips12 = 0;
01064 
01065       int iselect3 = 0;
01066       if (ndof >=15 && chisq <30) iselect3 = 1;
01067       
01068       if (isect <0) continue; //FIXGM Is it proper place ?
01069       if (fabs(trkth-pival/2)<0.000001) continue;   //22OCT07
01070 
01071       int ieta = int((abs(isect)%10000)/100.)-30; //an offset to acodate -ve eta values
01072       if (abs(ieta)>=16) continue;
01073       int iphi = abs(isect)%100;
01074 
01075       int tmpsect = int((iphi + 1)/6.) + 1;
01076       if (tmpsect >12) tmpsect = 1;
01077       
01078       int iring = 0;
01079       int tmpeta = ieta + 4; //For pixel mapping
01080       if (ieta >=-15 && ieta <=-11) {iring = -2; tmpeta =-11-ieta; } //abs(ieta)-11;} 
01081       if (ieta >=-10 && ieta <=-5)  {iring = -1; tmpeta =-5-ieta; } // abs(ieta)-5;}
01082       if (ieta >=  5 && ieta <= 10) {iring = 1; tmpeta  =ieta-5; }    
01083       if (ieta >= 11 && ieta <= 15) {iring = 2; tmpeta  =ieta-11;}   
01084       
01085       int iring2 = iring + 2;
01086       
01087       int tmprout = int((iphi + 1)/3.) + 1;
01088       int tmproutmx =routmx; 
01089       if (iring==0) {
01090         tmprout = int((iphi + 1)/2.) + 1;
01091         if (tmprout >routmx) tmprout = 1;
01092       } else {
01093         if (tmprout >rout12mx) tmprout = 1;
01094         tmproutmx =rout12mx; 
01095       }
01096 
01097       // CRUZET1
01098       if (m_cosmic) {
01099         /*  GMA temoparily change to increase event size at 3 & 9 O'clock position */
01100         if (abs(ndof) >=20 && abs(ndof) <40) {ips0 = (int)mypow_2_0; ipsall += ips0;}
01101         if (chisq >0 && chisq<15) {ips1 = (int)mypow_2_1; ipsall +=ips1;} //18Jan2008
01102         if (fabs(trkth-pival/2) <21.5) {ips2 = (int)mypow_2_2; ipsall +=ips2;} //No nead for pp evt
01103         if (fabs(trkph+pival/2) <21.5) {ips3 = (int)mypow_2_3; ipsall +=ips3;} //No nead for pp evt
01104         
01105         if (therr <0.02)             {ips4 = (int)mypow_2_4; ipsall +=ips4;}
01106         if (pherr <0.0002)             {ips5 = (int)mypow_2_5; ipsall +=ips5;}
01107         if (fabs(hoang) >0.30)             {ips6 = (int)mypow_2_6;  ipsall +=ips6;}
01108         if (fabs(trkmm) >0.100)        {ips7 = (int)mypow_2_7; ipsall +=ips7;}
01109         //      if (nmuon ==1)               {ips8 = (int)mypow_2_8;  ipsall +=ips8;}
01110         if (nmuon >=1 && nmuon <=4)        {ips8 = (int)mypow_2_8;  ipsall +=ips8;}
01111 
01112         if (iring2==2) {
01113           if (fabs(hodx)<100 && fabs(hodx)>2 && fabs(hocorsig[8]) <40 && fabs(hocorsig[8]) >2 )
01114             {ips10=(int)mypow_2_10;ipsall +=ips10;}
01115           
01116           if (fabs(hody)<100 && fabs(hody)>2 && fabs(hocorsig[9]) <40 && fabs(hocorsig[9]) >2 )
01117             {ips11=(int)mypow_2_11;ipsall +=ips11;}
01118           
01119         } else {
01120           if (fabs(hodx)<100 && fabs(hodx)>2 )
01121             {ips10=(int)mypow_2_10;ipsall +=ips10;}
01122           
01123           if (fabs(hody)<100 && fabs(hody)>2)
01124             {ips11=(int)mypow_2_11;ipsall +=ips11;}
01125         }
01126         if (caloen[0] ==0) { ips12=(int)mypow_2_12;ipsall +=ips12;}
01127       } else {
01128         //csa08
01129         if (abs(ndof) >=20 && abs(ndof) <40) {ips0 = (int)mypow_2_0; ipsall += ips0;}
01130         if (chisq >0 && chisq<15) {ips1 = (int)mypow_2_1; ipsall +=ips1;} //18Jan2008
01131         if (fabs(trkth-pival/2) <21.5) {ips2 = (int)mypow_2_2; ipsall +=ips2;} //No nead for pp evt
01132         if (fabs(trkph+pival/2) <21.5) {ips3 = (int)mypow_2_3; ipsall +=ips3;} //No nead for pp evt
01133         
01134         if (therr <0.02)             {ips4 = (int)mypow_2_4; ipsall +=ips4;}
01135         if (pherr <0.0002)             {ips5 = (int)mypow_2_5; ipsall +=ips5;}
01136         if (fabs(hoang) >0.30)             {ips6 = (int)mypow_2_6;  ipsall +=ips6;}
01137         if (fabs(trkmm) >4.0)        {ips7 = (int)mypow_2_7; ipsall +=ips7;}
01138         if (nmuon >=1 && nmuon <=2)        {ips8 = (int)mypow_2_8;  ipsall +=ips8;}
01139 
01140         if (iring2==2) {
01141           if (fabs(hodx)<100 && fabs(hodx)>2 && fabs(hocorsig[8]) <40 && fabs(hocorsig[8]) >2 )
01142             {ips10=(int)mypow_2_10;ipsall +=ips10;}
01143           
01144           if (fabs(hody)<100 && fabs(hody)>2 && fabs(hocorsig[9]) <40 && fabs(hocorsig[9]) >2 )
01145             {ips11=(int)mypow_2_11;ipsall +=ips11;}
01146           
01147         } else {
01148           if (fabs(hodx)<100 && fabs(hodx)>2 )
01149             {ips10=(int)mypow_2_10;ipsall +=ips10;}
01150           
01151           if (fabs(hody)<100 && fabs(hody)>2)
01152             {ips11=(int)mypow_2_11;ipsall +=ips11;}      
01153         }
01154         //      if (m_cosmic || (caloen[0] >0.5 && caloen[0]<5.0)) {ips12=(int)pow_2_12;ipsall +=ips12;}
01155         if (ndof >0 && caloen[0]<5.0) {ips12=(int)mypow_2_12;ipsall +=ips12;}
01156         /*      */
01157       }
01158       
01159       if (m_digiInput) {
01160         if (htime >2.5 && htime <6.5)        {ips9=(int)mypow_2_9;ipsall +=ips9;}
01161       } else {
01162         if (htime >-40 && htime <60)         {ips9=(int)mypow_2_9;ipsall +=ips9;}
01163       }    
01164 
01165       if (ipsall-ips0==mypow_2_ncut-mypow_2_0-1) sigvsevt[iring2][0]->Fill(abs(ndof), hosig[4]);
01166       if (ipsall-ips1==mypow_2_ncut-mypow_2_1-1) sigvsevt[iring2][1]->Fill(chisq, hosig[4]);
01167       if (ipsall-ips2==mypow_2_ncut-mypow_2_2-1) sigvsevt[iring2][2]->Fill(trkth, hosig[4]);
01168       if (ipsall-ips3==mypow_2_ncut-mypow_2_3-1) sigvsevt[iring2][3]->Fill(trkph, hosig[4]);
01169       if (ipsall-ips4==mypow_2_ncut-mypow_2_4-1) sigvsevt[iring2][4]->Fill(therr, hosig[4]);
01170       if (ipsall-ips5==mypow_2_ncut-mypow_2_5-1) sigvsevt[iring2][5]->Fill(pherr, hosig[4]);
01171       if (ipsall-ips6==mypow_2_ncut-mypow_2_6-1) sigvsevt[iring2][6]->Fill(hoang, hosig[4]);
01172       if (ipsall-ips7==mypow_2_ncut-mypow_2_7-1) sigvsevt[iring2][7]->Fill(fabs(trkmm), hosig[4]);
01173       if (ipsall-ips8==mypow_2_ncut-mypow_2_8-1) sigvsevt[iring2][8]->Fill(nmuon, hosig[4]);
01174       if (ipsall-ips9==mypow_2_ncut-mypow_2_9-1) sigvsevt[iring2][9]->Fill(htime, hosig[4]);
01175       if (ipsall-ips10==mypow_2_ncut-mypow_2_10-1) sigvsevt[iring2][10]->Fill(hodx, hosig[4]);
01176       if (ipsall-ips11==mypow_2_ncut-mypow_2_11-1) sigvsevt[iring2][11]->Fill(hody, hosig[4]);
01177       if (!m_cosmic) {
01178         if (ipsall-ips12==mypow_2_ncut-mypow_2_12-1) sigvsevt[iring2][12]->Fill(caloen[0], hosig[4]);
01179       }
01180       
01181       sigvsevt[iring2+5][0]->Fill(abs(ndof), hosig[4]);               
01182       if (ips0 >0) {
01183         sigvsevt[iring2+5][1]->Fill(chisq, hosig[4]);   
01184         if (ips1 >0) {
01185           sigvsevt[iring2+5][2]->Fill(trkth, hosig[4]);        
01186           if (ips2 >0) {
01187             sigvsevt[iring2+5][3]->Fill(trkph, hosig[4]);        
01188             if (ips3 >0) {
01189               sigvsevt[iring2+5][4]->Fill(therr, hosig[4]);        
01190               if (ips4 >0) {
01191                 sigvsevt[iring2+5][5]->Fill(pherr, hosig[4]);        
01192                 if (ips5 >0) {
01193                   sigvsevt[iring2+5][6]->Fill(hoang, hosig[4]);        
01194                   if (ips6 >0) {
01195                     sigvsevt[iring2+5][7]->Fill(fabs(trkmm), hosig[4]);   
01196                     if (ips7 >0) {
01197                       sigvsevt[iring2+5][8]->Fill(nmuon, hosig[4]); 
01198                       if (ips8 >0) {
01199                         sigvsevt[iring2+5][9]->Fill(htime, hosig[4]);
01200                         if (ips9 >0) {
01201                           sigvsevt[iring2+5][10]->Fill(hodx, hosig[4]);
01202                           if (ips10>0) {
01203                             sigvsevt[iring2+5][11]->Fill(hody, hosig[4]); 
01204                             if (ips11>0) {
01205                               if (!m_cosmic) sigvsevt[iring2+5][12]->Fill(caloen[0], hosig[4]);  
01206                             }
01207                           }
01208                         }
01209                       }
01210                     }
01211                   }
01212                 }
01213               }
01214             }
01215           }
01216         }
01217       }
01218       
01219       sigvsevt[iring2+10][0]->Fill(abs(ndof), hosig[4]);               
01220       sigvsevt[iring2+10][1]->Fill(chisq, hosig[4]);   
01221       sigvsevt[iring2+10][2]->Fill(trkth, hosig[4]);        
01222       sigvsevt[iring2+10][3]->Fill(trkph, hosig[4]);        
01223       sigvsevt[iring2+10][4]->Fill(therr, hosig[4]);        
01224       sigvsevt[iring2+10][5]->Fill(pherr, hosig[4]);        
01225       sigvsevt[iring2+10][6]->Fill(hoang, hosig[4]);        
01226       sigvsevt[iring2+10][7]->Fill(fabs(trkmm), hosig[4]);   
01227       sigvsevt[iring2+10][8]->Fill(nmuon, hosig[4]); 
01228       sigvsevt[iring2+10][9]->Fill(htime, hosig[4]);
01229       sigvsevt[iring2+10][10]->Fill(hodx, hosig[4]);
01230       sigvsevt[iring2+10][11]->Fill(hody, hosig[4]);     
01231       if (!m_cosmic) sigvsevt[iring2+10][12]->Fill(caloen[0], hosig[4]);  
01232 
01233       int iselect = (ipsall == mypow_2_ncut - 1) ? 1 : 0;
01234 
01235       if (hocro !=-100.0 && hocro <-50.0) hocro +=100.;
01236       
01237       muonnm->Fill(nmuon);
01238       muonmm->Fill(trkmm);
01239       muonth->Fill(trkth*180/pival);
01240       muonph->Fill(trkph*180/pival);
01241       muonch->Fill(chisq);
01242       
01243       if (iselect==1) { 
01244         ipass++;
01245         sel_muonnm->Fill(nmuon);
01246         sel_muonmm->Fill(trkmm);
01247         sel_muonth->Fill(trkth*180/pival);
01248         sel_muonph->Fill(trkph*180/pival);
01249         sel_muonch->Fill(chisq);
01250       }
01251 
01252       if (iselect3) T1->Fill();
01253 
01254       int tmpphi = (iphi + 1)%3; //pixel mapping
01255       int npixel = 0;
01256       int itag = -1;
01257       int iflip = 0;
01258       int fact = 2;
01259       
01260       if (iring ==0) { 
01261         tmpphi = (iphi+1)%2;
01262         if (tmpsect==2 || tmpsect==3 || tmpsect==6 || tmpsect==7 || tmpsect==10 || tmpsect==11) {
01263           npixel = mapx0p[tmpeta][tmpphi];
01264           itag = 2;
01265         } else {
01266           npixel = mapx0m[tmpeta][tmpphi];
01267           itag = 3;
01268         }
01269       } else { 
01270         fact = 3;
01271         if (tmpsect%2==1) iflip =1;
01272         if (abs(iring)==1) {
01273           npixel = mapx1[tmpeta][(iflip==0) ? tmpphi : abs(tmpphi-2)];
01274           itag = 1;
01275         } else {
01276           npixel = mapx2[tmpeta][(iflip==0) ? tmpphi : abs(tmpphi-2)];
01277           itag = 0;
01278         }
01279       }
01280 
01281       int tmpeta1 = (ieta>0) ? ieta -1 : -ieta +14; 
01282 
01283       //Histogram filling for noise study: phi shift according to DTChamberAnalysis
01284       int tmpphi1 = iphi -1 ; //(iphi+6 <=nphimx) ? iphi+5 : iphi+5-nphimx;
01285 
01286       int iselect2=0;
01287       if (hosig[4]!=-100) { 
01288         if (m_cosmic) {
01289           if (caloen[2]<=0.0) iselect2=1;
01290         } else {
01291           if (caloen[2]<=3.0) iselect2=1;
01292         }  
01293       }
01294 
01295       //      cout <<"cosmic "<<hosig[4]<<" "<<caloen[3]<<" "<<int(iselect2)<<" "<<int(m_cosmic)<<endl;
01296 
01297 
01298       if (iselect2==1) {
01299         int tmpphi2 = iphi - 1;
01300         if (!m_digiInput) { tmpphi2 =  (iphi+6 <=nphimx) ? iphi+5 : iphi+5-nphimx;}
01301         
01302         int tmpsect2 = int((tmpphi2 + 2)/6.) + 1;
01303         if (tmpsect2 >12) tmpsect2 = 1;
01304         
01305         int tmprout2 = int((tmpphi2 + 2)/3.) + 1;
01306         if (iring==0) {
01307           tmprout2 = int((tmpphi2 + 2)/2.) + 1;
01308           if (tmprout2 >routmx) tmprout2 = 1;
01309         } else {
01310           if (tmprout2 >rout12mx) tmprout2 = 1;
01311         }
01312 
01313         if (cro_ssg[tmpeta1][tmpphi2].size()<4000) {
01314           if (hocro>alow && hocro<ahigh) {
01315             if (!m_histfit) cro_ssg[tmpeta1][tmpphi2].push_back(hocro);
01316             crossg[tmpeta1][tmpphi2]->Fill(hocro);
01317           }
01318         }
01319 
01320         if (tmpphi2 >=0 && tmpphi2 <nphimx) {
01321           crossg[tmpeta1][nphimx]->Fill(hocro);
01322         }
01323         if (m_combined) {
01324           com_crossg[iring2][tmprout2-1]->Fill(hocro); 
01325           com_crossg[iring2][tmproutmx]->Fill(hocro); 
01326         }
01327       }
01328 
01329       if (iselect==1) { 
01330         for (int ij=0; ij<neffip; ij++) {
01331           if (ij==0) {
01332             sig_effi[ij]->Fill(ieta, iphi, 1.);
01333           } else {
01334             if (hosig[4] >ij*m_sigma) {
01335               sig_effi[ij]->Fill(ieta, iphi, 1.);
01336             } 
01337           }
01338         }
01339 
01340         tmpphi1 = iphi - 1;
01341 
01342         if (sig_reg[tmpeta1][tmpphi1].size()<4000 ) {
01343           if (hosig[4]>-50&& hosig[4] <15) {
01344             sigrsg[tmpeta1][tmpphi1]->Fill(hosig[4]);  
01345             if (!m_histfit && hosig[4]<=ahigh/2.) sig_reg[tmpeta1][tmpphi1].push_back(hosig[4]);
01346             invang[tmpeta1][tmpphi1] += 1./fabs(hoang);
01347           }
01348         }
01349         
01350         if (tmpphi1 >=0 && tmpphi1 <nphimx) { //GREN
01351           sigrsg[tmpeta1][nphimx]->Fill(hosig[4]);  
01352           invang[tmpeta1][nphimx] += 1./fabs(hoang);
01353         }
01354 
01355         if (m_combined) { //#ifdef COMBINED
01356           com_sigrsg[iring2][tmprout-1]->Fill(hosig[4]); 
01357           com_invang[iring2][tmprout-1] += 1./fabs(hoang);
01358           
01359           com_sigrsg[iring2][tmproutmx]->Fill(hosig[4]); 
01360           com_invang[iring2][tmproutmx] += 1./fabs(hoang);
01361         } //m_combined #endif
01362         
01363         if (m_checkmap || m_correl) { //#ifdef CHECKMAP
01364           tmpeta = etamap[itag][npixel];
01365           tmpphi = phimap[itag][npixel];
01366           if (tmpeta>=0 && tmpphi >=0) {
01367             if (iflip !=0) tmpphi = abs(tmpphi-2);
01368             if (int((hocorsig[fact*tmpeta+tmpphi]-hosig[4])*10000)/10000.!=0) {
01369               iaxxx++;
01370               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;
01371               
01372               for (int i=0; i<18; i++) {cout <<" "<<i<<" "<<hocorsig[i];}
01373               cout<<" ix "<<iaxxx<<" "<<ibxxx<<endl;
01374             } else { ibxxx++; }
01375             
01376             corrsgc[tmpeta1][tmpphi1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
01377             if (m_combined) { //#ifdef COMBINED
01378               com_corrsgc[iring2][tmpsect-1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
01379             } //m_combined #endif
01380           }
01381         } //m_checkmap #endif
01382         
01383         if (m_correl) { //#ifdef CORREL
01384           float allcorsig = 0.0;
01385           
01386           tmpeta = etamap[itag][npixleft[npixel]];
01387           tmpphi = phimap[itag][npixleft[npixel]];
01388           
01389           if (tmpeta>=0 && tmpphi >=0) {
01390             if (iflip !=0) tmpphi = abs(tmpphi-2);
01391             corrsgl[tmpeta1][tmpphi1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
01392             allcorsig +=hocorsig[fact*tmpeta+tmpphi];
01393             if (m_combined) { //#ifdef COMBINED
01394               com_corrsgl[iring2][tmpsect-1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
01395             } //m_combined #endif
01396           }         
01397           
01398           tmpeta = etamap[itag][npixrigh[npixel]];
01399           tmpphi = phimap[itag][npixrigh[npixel]];
01400           if (tmpeta>=0 && tmpphi >=0) {
01401             if (iflip !=0) tmpphi = abs(tmpphi-2);
01402             corrsgr[tmpeta1][tmpphi1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
01403             allcorsig +=hocorsig[fact*tmpeta+tmpphi];
01404             if (m_combined) { // #ifdef COMBINED
01405               com_corrsgr[iring2][tmpsect-1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
01406             } //m_combined #endif
01407           }
01408           
01409           tmpeta = etamap[itag][npixlebt[npixel]];
01410           tmpphi = phimap[itag][npixlebt[npixel]];
01411           if (tmpeta>=0 && tmpphi >=0) {
01412             if (iflip !=0) tmpphi = abs(tmpphi-2);
01413             corrsglb[tmpeta1][tmpphi1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
01414             allcorsig +=hocorsig[fact*tmpeta+tmpphi];
01415             if (m_combined){ //#ifdef COMBINED
01416               com_corrsglb[iring2][tmpsect-1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
01417             } //m_combined #endif
01418           }
01419           
01420           tmpeta = etamap[itag][npixribt[npixel]];
01421           tmpphi = phimap[itag][npixribt[npixel]];
01422           if (tmpeta>=0 && tmpphi >=0) {
01423             if (iflip !=0) tmpphi = abs(tmpphi-2);
01424             corrsgrb[tmpeta1][tmpphi1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
01425             allcorsig +=hocorsig[fact*tmpeta+tmpphi];
01426             if (m_combined) { // #ifdef COMBINED
01427               com_corrsgrb[iring2][tmpsect-1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
01428             } //m_combined #endif
01429           }
01430           
01431           tmpeta = etamap[itag][npixleup[npixel]];
01432           tmpphi = phimap[itag][npixleup[npixel]];
01433           if (tmpeta>=0 && tmpphi >=0) {
01434             if (iflip !=0) tmpphi = abs(tmpphi-2);
01435             corrsglu[tmpeta1][tmpphi1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
01436             allcorsig +=hocorsig[fact*tmpeta+tmpphi];
01437             if (m_combined) {// #ifdef COMBINED
01438               com_corrsglu[iring2][tmpsect-1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
01439             } //m_combined #endif
01440           }
01441           
01442           tmpeta = etamap[itag][npixriup[npixel]];
01443           tmpphi = phimap[itag][npixriup[npixel]];
01444           if (tmpeta>=0 && tmpphi >=0) {
01445             if (iflip !=0) tmpphi = abs(tmpphi-2);
01446             corrsgru[tmpeta1][tmpphi1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
01447             allcorsig +=hocorsig[fact*tmpeta+tmpphi];
01448             if (m_combined) { // #ifdef COMBINED
01449               com_corrsgru[iring2][tmpsect-1]->Fill(hocorsig[fact*tmpeta+tmpphi]);
01450             } //m_combined #endif
01451           }
01452           corrsgall[tmpeta1][tmpphi1]->Fill(allcorsig); 
01453           if (m_combined) { // #ifdef COMBINED
01454             com_corrsgall[iring2][tmpsect-1]->Fill(allcorsig);
01455           } //m_combined #endif
01456           
01457           
01458         } //m_correl #endif
01459         for (int k=0; k<9; k++) {
01460           switch (iring) 
01461             {
01462             case 2 : ho_sig2p[k]->Fill(hosig[k]); break;
01463             case 1 : ho_sig1p[k]->Fill(hosig[k]); break;
01464             case 0 : ho_sig00[k]->Fill(hosig[k]); break;
01465             case -1 : ho_sig1m[k]->Fill(hosig[k]); break;
01466             case -2 : ho_sig2m[k]->Fill(hosig[k]); break;
01467             }     
01468           if (m_hbinfo) { // #ifdef HBINFO
01469             hbhe_sig[k]->Fill(hbhesig[k]);
01470             //      cout <<"hbhe "<<k<<" "<<hbhesig[k]<<endl;
01471           } //m_hbinfo #endif
01472         }
01473       } //if (iselect==1)
01474     } //for (HOCalibVariableCollection::const_iterator hoC=(*HOCalib).begin(); hoC!=(*HOCalib).end(); hoC++){
01475     
01476   } //if (isCosMu) 
01477   
01478 #ifdef THIS_IS_AN_EVENT_EXAMPLE
01479    Handle<ExampleData> pIn;
01480    iEvent.getByLabel("example",pIn);
01481 #endif
01482    
01483 #ifdef THIS_IS_AN_EVENTSETUP_EXAMPLE
01484    ESHandle<SetupData> pSetup;
01485    iSetup.get<SetupRecord>().get(pSetup);
01486 #endif
01487 }
01488 
01489 
01490 // ------------ method called once each job just before starting event loop  ------------
01491 void 
01492 HOCalibAnalyzer::beginJob()
01493 {
01494 }
01495 
01496 // ------------ method called once each job just after ending the event loop  ------------
01497 void 
01498 HOCalibAnalyzer::endJob() {
01499 
01500   theFile->cd();
01501   
01502   for (int i=0; i<nphimx; i++) {
01503     for (int j=0; j<netamx; j++) {
01504       
01505       nevsigrsg->Fill(netamx*i+j, sigrsg[j][i]->GetEntries());
01506       mnsigrsg->Fill(netamx*i+j, sigrsg[j][i]->GetMean());
01507       rmssigrsg->Fill(netamx*i+j, sigrsg[j][i]->GetRMS());
01508       
01509       nevcrossg->Fill(netamx*i+j, crossg[j][i]->GetEntries());
01510       mncrossg->Fill(netamx*i+j, crossg[j][i]->GetMean());
01511       rmscrossg->Fill(netamx*i+j, crossg[j][i]->GetRMS());
01512       
01513       if (m_correl) { //#ifdef CORREL    
01514         
01515         nevcorrsglb->Fill(netamx*i+j, corrsglb[j][i]->GetEntries());
01516         mncorrsglb->Fill(netamx*i+j, corrsglb[j][i]->GetMean());
01517         rmscorrsglb->Fill(netamx*i+j, corrsglb[j][i]->GetRMS());
01518         
01519         nevcorrsgrb->Fill(netamx*i+j, corrsgrb[j][i]->GetEntries());
01520         mncorrsgrb->Fill(netamx*i+j, corrsgrb[j][i]->GetMean());
01521         rmscorrsgrb->Fill(netamx*i+j, corrsgrb[j][i]->GetRMS());
01522         
01523         nevcorrsglu->Fill(netamx*i+j, corrsglu[j][i]->GetEntries());
01524         mncorrsglu->Fill(netamx*i+j, corrsglu[j][i]->GetMean());
01525         rmscorrsglu->Fill(netamx*i+j, corrsglu[j][i]->GetRMS());
01526 
01527         nevcorrsgru->Fill(netamx*i+j, corrsgru[j][i]->GetEntries());
01528         mncorrsgru->Fill(netamx*i+j, corrsgru[j][i]->GetMean());
01529         rmscorrsgru->Fill(netamx*i+j, corrsgru[j][i]->GetRMS());
01530         
01531         nevcorrsgall->Fill(netamx*i+j, corrsgall[j][i]->GetEntries());
01532         mncorrsgall->Fill(netamx*i+j, corrsgall[j][i]->GetMean());
01533         rmscorrsgall->Fill(netamx*i+j, corrsgall[j][i]->GetRMS());
01534         
01535         nevcorrsgl->Fill(netamx*i+j, corrsgl[j][i]->GetEntries());
01536         mncorrsgl->Fill(netamx*i+j, corrsgl[j][i]->GetMean());
01537         rmscorrsgl->Fill(netamx*i+j, corrsgl[j][i]->GetRMS());
01538         
01539         nevcorrsgr->Fill(netamx*i+j, corrsgr[j][i]->GetEntries());
01540         mncorrsgr->Fill(netamx*i+j, corrsgr[j][i]->GetMean());
01541         rmscorrsgr->Fill(netamx*i+j, corrsgr[j][i]->GetRMS());
01542       } //m_correl #endif
01543       if (m_checkmap) { //#ifdef CHECKMAP
01544         nevcorrsgc->Fill(netamx*i+j, corrsgc[j][i]->GetEntries());
01545         mncorrsgc->Fill(netamx*i+j, corrsgc[j][i]->GetMean());
01546         rmscorrsgc->Fill(netamx*i+j, corrsgc[j][i]->GetRMS());
01547       } //m_checkmap #endif
01548     }
01549   }
01550 
01551   if (m_combined) { // #ifdef COMBINED
01552     for (int j=0; j<ringmx; j++) {
01553       for (int i=0; i<routmx; i++) {
01554         if (j!=2 && i>=rout12mx) continue;
01555         nevsigrsg->Fill(netamx*nphimx+ringmx*i+j, com_sigrsg[j][i]->GetEntries());
01556         mnsigrsg->Fill(netamx*nphimx+ringmx*i+j, com_sigrsg[j][i]->GetMean());
01557         rmssigrsg->Fill(netamx*nphimx+ringmx*i+j, com_sigrsg[j][i]->GetRMS());
01558         
01559         nevcrossg->Fill(netamx*nphimx+ringmx*i+j, com_crossg[j][i]->GetEntries());
01560         mncrossg->Fill(netamx*nphimx+ringmx*i+j, com_crossg[j][i]->GetMean());
01561         rmscrossg->Fill(netamx*nphimx+ringmx*i+j, com_crossg[j][i]->GetRMS());
01562       }
01563     }
01564 
01565     for (int i=0; i<sectmx; i++) {
01566       for (int j=0; j<ringmx; j++) {
01567         if (m_correl) { // #ifdef CORREL      
01568           nevcorrsglb->Fill(netamx*nphimx+ringmx*i+j, com_corrsglb[j][i]->GetEntries());
01569           mncorrsglb->Fill(netamx*nphimx+ringmx*i+j, com_corrsglb[j][i]->GetMean());
01570           rmscorrsglb->Fill(netamx*nphimx+ringmx*i+j, com_corrsglb[j][i]->GetRMS());
01571           
01572           nevcorrsgrb->Fill(netamx*nphimx+ringmx*i+j, com_corrsgrb[j][i]->GetEntries());
01573           mncorrsgrb->Fill(netamx*nphimx+ringmx*i+j, com_corrsgrb[j][i]->GetMean());
01574           rmscorrsgrb->Fill(netamx*nphimx+ringmx*i+j, com_corrsgrb[j][i]->GetRMS());
01575           
01576           nevcorrsglu->Fill(netamx*nphimx+ringmx*i+j, com_corrsglu[j][i]->GetEntries());
01577           mncorrsglu->Fill(netamx*nphimx+ringmx*i+j, com_corrsglu[j][i]->GetMean());
01578           rmscorrsglu->Fill(netamx*nphimx+ringmx*i+j, com_corrsglu[j][i]->GetRMS());
01579           
01580           nevcorrsgru->Fill(netamx*nphimx+ringmx*i+j, com_corrsgru[j][i]->GetEntries());
01581           mncorrsgru->Fill(netamx*nphimx+ringmx*i+j, com_corrsgru[j][i]->GetMean());
01582           rmscorrsgru->Fill(netamx*nphimx+ringmx*i+j, com_corrsgru[j][i]->GetRMS());
01583           
01584           nevcorrsgall->Fill(netamx*nphimx+ringmx*i+j, com_corrsgall[j][i]->GetEntries());
01585           mncorrsgall->Fill(netamx*nphimx+ringmx*i+j, com_corrsgall[j][i]->GetMean());
01586           rmscorrsgall->Fill(netamx*nphimx+ringmx*i+j, com_corrsgall[j][i]->GetRMS());
01587           
01588           nevcorrsgl->Fill(netamx*nphimx+ringmx*i+j, com_corrsgl[j][i]->GetEntries());
01589           mncorrsgl->Fill(netamx*nphimx+ringmx*i+j, com_corrsgl[j][i]->GetMean());
01590           rmscorrsgl->Fill(netamx*nphimx+ringmx*i+j, com_corrsgl[j][i]->GetRMS());
01591           
01592           nevcorrsgr->Fill(netamx*nphimx+ringmx*i+j, com_corrsgr[j][i]->GetEntries());
01593           mncorrsgr->Fill(netamx*nphimx+ringmx*i+j, com_corrsgr[j][i]->GetMean());
01594           rmscorrsgr->Fill(netamx*nphimx+ringmx*i+j, com_corrsgr[j][i]->GetRMS());
01595         } //m_correl #endif
01596         if (m_checkmap) { // #ifdef CHECKMAP
01597           nevcorrsgc->Fill(netamx*nphimx+ringmx*i+j, com_corrsgc[j][i]->GetEntries());
01598           mncorrsgc->Fill(netamx*nphimx+ringmx*i+j, com_corrsgc[j][i]->GetMean());
01599           rmscorrsgc->Fill(netamx*nphimx+ringmx*i+j, com_corrsgc[j][i]->GetRMS());
01600         } //m_checkmap #endif
01601       }
01602     }
01603   } //m_combined #endif
01604 
01605 
01606   for (int i=1; i<neffip; i++) {
01607     sig_effi[i]->Divide(sig_effi[0]);
01608   }
01609   for (int ij=0; ij<netamx; ij++) {
01610     for (int jk = 0; jk <nphimx; jk++) {
01611       int ieta = (ij<15) ? ij+1 : 14-ij;
01612       int iphi = jk+1;
01613       double signal = sigrsg[ij][jk]->GetMean();
01614       mean_energy->Fill(ieta, iphi, signal);
01615     }
01616   }      
01617 
01618  int irunold = irun;
01619 
01620 
01621   gStyle->SetOptLogy(0);
01622   gStyle->SetTitleFillColor(10);
01623   gStyle->SetStatColor(10);
01624   
01625   gStyle->SetCanvasColor(10);
01626   gStyle->SetOptStat(0); //1110);
01627   gStyle->SetOptTitle(1);
01628 
01629   gStyle->SetTitleColor(10);
01630   gStyle->SetTitleFontSize(0.09);
01631   gStyle->SetTitleOffset(-0.05);
01632   gStyle->SetTitleBorderSize(1);
01633 
01634   gStyle->SetPadColor(10);
01635   gStyle->SetPadBorderMode(0);
01636   gStyle->SetStatColor(10);
01637   gStyle->SetPadBorderMode(0);
01638   gStyle->SetStatBorderSize(1);
01639   gStyle->SetStatFontSize(.07);
01640 
01641   gStyle->SetStatStyle(1001);
01642   gStyle->SetOptFit(101);
01643   gStyle->SetCanvasColor(10);
01644   gStyle->SetCanvasBorderMode(0);
01645 
01646   gStyle->SetStatX(.99);
01647   gStyle->SetStatY(.99);
01648   gStyle->SetStatW(.45);
01649   gStyle->SetStatH(.16);
01650   gStyle->SetLabelSize(0.075,"XY");  
01651   gStyle->SetLabelOffset(0.21,"XYZ");
01652   gStyle->SetTitleSize(0.065,"XY");  
01653   gStyle->SetTitleOffset(0.06,"XYZ");
01654   gStyle->SetPadTopMargin(.09);
01655   gStyle->SetPadBottomMargin(0.11);
01656   gStyle->SetPadLeftMargin(0.12);
01657   gStyle->SetPadRightMargin(0.15);
01658   gStyle->SetPadGridX(3);
01659   gStyle->SetPadGridY(3);
01660   gStyle->SetGridStyle(2);
01661   gStyle->SetNdivisions(303,"XY");
01662 
01663   gStyle->SetMarkerSize(0.60);
01664   gStyle->SetMarkerColor(2);
01665   gStyle->SetMarkerStyle(20);
01666   gStyle->SetTitleFontSize(0.07);
01667 
01668   char out_file[200];
01669   int xsiz = 700;
01670   int ysiz = 500;
01671 
01672   TCanvas *c2 = new TCanvas("c2", "Statistics and efficiency", xsiz, ysiz);
01673   c2->Divide(2,1); //(3,2);
01674   for (int ij=0; ij<neffip; ij=ij+3) {
01675     sig_effi[ij]->GetXaxis()->SetTitle("#eta");
01676     sig_effi[ij]->GetXaxis()->SetTitleSize(0.075);
01677     sig_effi[ij]->GetXaxis()->SetTitleOffset(0.65); //0.85 
01678     sig_effi[ij]->GetXaxis()->CenterTitle();
01679     sig_effi[ij]->GetXaxis()->SetLabelSize(0.055);
01680     sig_effi[ij]->GetXaxis()->SetLabelOffset(0.001); 
01681     
01682     sig_effi[ij]->GetYaxis()->SetTitle("#phi");
01683     sig_effi[ij]->GetYaxis()->SetTitleSize(0.075);
01684     sig_effi[ij]->GetYaxis()->SetTitleOffset(0.9); 
01685     sig_effi[ij]->GetYaxis()->CenterTitle();
01686     sig_effi[ij]->GetYaxis()->SetLabelSize(0.055);
01687     sig_effi[ij]->GetYaxis()->SetLabelOffset(0.01); 
01688     
01689     c2->cd(int(ij/3.)+1); sig_effi[ij]->Draw("colz");
01690   }
01691   sprintf(out_file, "comb_hosig_evt_%i.jpg",irunold); 
01692   c2->SaveAs(out_file); 
01693   
01694   gStyle->SetTitleFontSize(0.045);
01695   gStyle->SetPadRightMargin(0.1);
01696   gStyle->SetPadLeftMargin(0.1);
01697   gStyle->SetPadBottomMargin(0.12);
01698 
01699   TCanvas *c1 = new TCanvas("c1", "Mean signal in each tower", xsiz, ysiz);  
01700   
01701   mean_energy->GetXaxis()->SetTitle("#eta");
01702   mean_energy->GetXaxis()->SetTitleSize(0.075);
01703   mean_energy->GetXaxis()->SetTitleOffset(0.65); //0.6 
01704   mean_energy->GetXaxis()->CenterTitle();
01705   mean_energy->GetXaxis()->SetLabelSize(0.045);
01706   mean_energy->GetXaxis()->SetLabelOffset(0.001); 
01707   
01708   mean_energy->GetYaxis()->SetTitle("#phi");
01709   mean_energy->GetYaxis()->SetTitleSize(0.075);
01710   mean_energy->GetYaxis()->SetTitleOffset(0.5); 
01711   mean_energy->GetYaxis()->CenterTitle();
01712   mean_energy->GetYaxis()->SetLabelSize(0.045);
01713   mean_energy->GetYaxis()->SetLabelOffset(0.01); 
01714   
01715   mean_energy->Draw("colz");
01716   sprintf(out_file, "homean_energy_%i.jpg",irunold); 
01717   c1->SaveAs(out_file); 
01718   
01719   delete c1; 
01720   delete c2;
01721 
01722   gStyle->SetPadBottomMargin(0.14);
01723   gStyle->SetPadLeftMargin(0.17);
01724   gStyle->SetPadRightMargin(0.03);
01725 
01726   gStyle->SetOptStat(1110);
01727 
01728   const int nsample =8;  
01729   TF1*  gx0[nsample]={0};
01730   TF1* ped0fun[nsample]={0};
01731   TF1* signal[nsample]={0};
01732   TF1* pedfun[nsample]={0};
01733   TF1* sigfun[nsample]={0};
01734   TF1* signalx[nsample]={0};
01735   
01736   TH1F* signall[nsample]={0};
01737   TH1F* pedstll[nsample]={0};
01738 
01739   if (m_constant) { 
01740 
01741     gStyle->SetOptFit(101);
01742     gStyle->SetCanvasBorderMode(0);
01743     gStyle->SetPadBorderMode(0);
01744     gStyle->SetStatBorderSize(1);
01745     gStyle->SetStatStyle(1001);
01746     gStyle->SetTitleColor(10);
01747     gStyle->SetTitleFontSize(0.09);
01748     gStyle->SetTitleOffset(-0.05);
01749     gStyle->SetTitleBorderSize(1);
01750   
01751     gStyle->SetCanvasColor(10);
01752     gStyle->SetPadColor(10);
01753     gStyle->SetStatColor(10);
01754     gStyle->SetStatFontSize(.07);
01755     gStyle->SetStatX(0.99);
01756     gStyle->SetStatY(0.99);
01757     gStyle->SetStatW(0.30);
01758     gStyle->SetStatH(0.10);
01759     gStyle->SetTitleSize(0.065,"XYZ");
01760     gStyle->SetLabelSize(0.075,"XYZ");
01761     gStyle->SetLabelOffset(0.012,"XYZ");
01762     gStyle->SetPadGridX(1);
01763     gStyle->SetPadGridY(1);
01764     gStyle->SetGridStyle(3);
01765     gStyle->SetNdivisions(101,"XY");
01766     gStyle->SetOptLogy(0);
01767     int iiter = 0;
01768 
01769 
01770     ofstream file_out(theoutputtxtFile.c_str());
01771     //    TPostScript* ps=0;
01772     int ips=111;
01773     TPostScript ps(theoutputpsFile.c_str(),ips);
01774     ps.Range(20,28);
01775 
01776     xsiz = 900; //900;
01777     ysiz = 1200; //600;
01778     TCanvas *c0 = new TCanvas("c0", " Pedestal vs signal", xsiz, ysiz);
01779 
01780     float mean_eta[nphimx];
01781     float mean_phi[netamx];
01782     float rms_eta[nphimx];
01783     float rms_phi[netamx];
01784 
01785     for (int ij=0; ij<nphimx; ij++) {mean_phi[ij] = rms_phi[ij] =0;}
01786     for (int ij=0; ij<netamx; ij++) {mean_eta[ij] = rms_eta[ij] =0;}
01787 
01788     int mxeta = 0;
01789     int mxphi = 0;
01790     int mneta = 0;
01791     int mnphi = 0;
01792 
01793     //iijj = 0 : Merging all ring
01794     //     = 1 : Individual HPD
01795     //iijj = 2 : merging all phi
01796     //     = 3 : Individual tower
01797 
01798     for (int iijj = 0; iijj <4; iijj++) {
01799       //      if ((!mx_combined) && iijj==1) continue; //Use this only for combined data  
01800       if (iijj==0){
01801         mxeta = ringmx; mxphi = 1;      mneta = 0; mnphi = 0;
01802       } else if (iijj==1) {
01803         mxeta = ringmx; mxphi = routmx;
01804         mneta = 0; mnphi = 0;
01805       } else if (iijj==2) {
01806         mxeta = netamx; mxphi = 1; mneta = 0; mnphi = 0;  
01807       } else if (iijj==3) {
01808         mxeta = netamx; mxphi = nphimx;
01809         mneta = 0; mnphi = 0;
01810       }
01811       
01812       for (int jk=mneta; jk<mxeta; jk++) {
01813         for (int ij=mnphi; ij<mxphi; ij++) {
01814           if (iijj==1) continue;
01815           if ((iijj==0 || iijj==1) && jk !=2 && ij >=rout12mx) continue;
01816           int izone = iiter%nsample;
01817 
01818           if (iijj==0) {
01819             int iread = (jk==2) ? routmx : rout12mx;
01820             signall[izone] = (TH1F*)com_sigrsg[jk][iread]->Clone("hnew");
01821             pedstll[izone] = (TH1F*)com_crossg[jk][iread]->Clone("hnew");
01822           } else if (iijj==1) {
01823             signall[izone] = (TH1F*)com_sigrsg[jk][ij]->Clone("hnew");
01824             pedstll[izone] = (TH1F*)com_crossg[jk][ij]->Clone("hnew");
01825           } else if (iijj==2) {
01826             signall[izone] = (TH1F*)sigrsg[jk][nphimx]->Clone("hnew");
01827             pedstll[izone] = (TH1F*)crossg[jk][nphimx]->Clone("hnew");
01828           } else if (iijj==3) {
01829             signall[izone] = (TH1F*)sigrsg[jk][ij]->Clone("hnew");
01830             pedstll[izone] = (TH1F*)crossg[jk][ij]->Clone("hnew");
01831           }
01832 
01833           pedstll[izone]->SetLineWidth(2);
01834           signall[izone]->SetLineWidth(2);
01835           pedstll[izone]->SetLineColor(2);
01836           signall[izone]->SetLineColor(4);
01837           pedstll[izone]->SetNdivisions(506,"XY");
01838           signall[izone]->SetNdivisions(506,"XY");
01839 
01840           signall[izone]->GetXaxis()->SetLabelSize(.065);
01841           signall[izone]->GetYaxis()->SetLabelSize(.06);            
01842           if (m_digiInput) {
01843             signall[izone]->GetXaxis()->SetTitle("Signal (fC)");
01844           } else {
01845             signall[izone]->GetXaxis()->SetTitle("Signal (GeV)");
01846           }
01847           signall[izone]->GetXaxis()->SetTitleSize(.065);
01848           signall[izone]->GetXaxis()->CenterTitle(); 
01849           
01850           if (izone==0) { //iiter%8 ==0) { 
01851             ps.NewPage();
01852             c0->Divide(4,4); //c0->Divide(2,4); // c0->Divide(1,2);
01853           }
01854           c0->cd(2*izone+1); // (iiter%8)+1); //c0->cd(iiter%8+1);
01855 
01856           /*
01857           if (iijj==0 && izone==0) {
01858             gStyle->SetOptLogy(1);
01859             gStyle->SetOptStat(0);
01860             gStyle->SetOptFit(0);
01861             c0->Divide(3,2);
01862           }
01863 
01864           if (iijj>0) {
01865             gStyle->SetOptLogy(0);
01866             gStyle->SetOptStat(1110);
01867             gStyle->SetOptFit(101);
01868             
01869             if (iiter==0) {
01870               int ips=111;
01871               ps = new TPostScript(theoutputpsFile.c_str(),ips);
01872               ps.Range(20,28);
01873               xsiz = 900; //900;
01874               ysiz = 1200; //600;
01875               c0 = new TCanvas("c0", " Pedestal vs signal", xsiz, ysiz);
01876             }
01877             if (izone==0) {
01878               ps.NewPage();
01879               c0->Divide(4,4);
01880             }
01881           }
01882           if (iijj==0) {c0->cd(izone+1); } else { c0->cd(2*izone+1);}
01883           */
01884 
01885           float mean = pedstll[izone]->GetMean();
01886           float rms = pedstll[izone]->GetRMS();
01887           if (m_digiInput) {
01888             if (rms <0.6) rms = 0.6;
01889             if (rms >1.2) rms = 1.2;
01890             if (mean >1.2) mean = 1.2;
01891             if (mean <-1.2) mean = -1.2;
01892           } else {
01893             if (rms <0.10) rms = 0.10;
01894             if (rms >0.15) rms=0.15;
01895             if (mean >0.20) mean = 0.20;
01896             if (mean <-0.20) mean = -0.20;
01897           }
01898           float xmn = mean-6.*rms;
01899           float xmx = mean+6.*rms;
01900           
01901           binwid =      pedstll[izone]->GetBinWidth(1);
01902           if (xmx > pedstll[izone]->GetXaxis()->GetXmax()) xmx = pedstll[izone]->GetXaxis()->GetXmax()-0.5*binwid;
01903           if (xmn < pedstll[izone]->GetXaxis()->GetXmin()) xmn = pedstll[izone]->GetXaxis()->GetXmin()+0.5*binwid;
01904           
01905           float height = pedstll[izone]->GetEntries();
01906           
01907           double par[nbgpr] ={height, mean, 0.75*rms};
01908 
01909           double gaupr[nbgpr];
01910           double parer[nbgpr];
01911           
01912           ietafit = jk;
01913           iphifit = ij;
01914           pedstll[izone]->GetXaxis()->SetLabelSize(.065);
01915           pedstll[izone]->GetYaxis()->SetLabelSize(.06);
01916 
01917           //      if (iijj==0) {
01918           //        pedstll[izone]->GetXaxis()->SetRangeUser(alow, ahigh);
01919           //      } else {
01920             pedstll[izone]->GetXaxis()->SetRangeUser(xmn, xmx);
01921             //    }
01922 
01923           if (m_digiInput) {
01924             if (iijj==0) {
01925               pedstll[izone]->GetXaxis()->SetTitle("Pedestal/Signal (fC)");
01926             } else {
01927               pedstll[izone]->GetXaxis()->SetTitle("Pedestal (fC)");
01928             }
01929           } else {
01930             if (iijj==0) {
01931               pedstll[izone]->GetXaxis()->SetTitle("Pedestal/Signal (GeV)");
01932             } else {
01933               pedstll[izone]->GetXaxis()->SetTitle("Pedestal (GeV)");
01934             }
01935           }
01936           pedstll[izone]->GetXaxis()->SetTitleSize(.065);
01937           pedstll[izone]->GetXaxis()->CenterTitle(); 
01938           //      pedstll[izone]->SetLineWidth(2);
01939 
01940           pedstll[izone]->Draw();
01941           if (m_pedsuppr && !m_digiInput) {
01942             gaupr[0] = 0;
01943             gaupr[1] = 0.0; // pedmean[ietafit][iphifit];
01944             if (m_digiInput) {
01945               gaupr[2] = 0.90; //GMA need from database pedwidth[ietafit][iphifit];
01946             } else {
01947               gaupr[2] = 0.15; //GMA need from database
01948             }
01949             parer[0] = parer[1] = parer[2] = 0;
01950           } else {
01951             
01952             if (pedstll[izone]->GetEntries() >5) {
01953 
01954               if ((iijj!=3) || m_histfit) {
01955                 char temp[20];
01956                 sprintf(temp, "gx0_%i",izone);
01957                 gx0[izone] = new TF1(temp, gausX, xmn, xmx, nbgpr);   
01958                 gx0[izone]->SetParameters(par);
01959                 gx0[izone]->SetLineWidth(1);
01960                 pedstll[izone]->Fit(gx0[izone], "R+");
01961                 
01962                 for (int k=0; k<nbgpr; k++) {
01963                   parer[k] = gx0[izone]->GetParError(k);
01964                   gaupr[k] = gx0[izone]->GetParameter(k);
01965                 }
01966               } else {
01967                 double strt[nbgpr] = {height, mean, 0.75*rms};
01968                 double step[nbgpr] = {1.0, 0.001, 0.001};
01969                 double alowmn[nbgpr] = {0.5*height, mean-rms, 0.3*rms};
01970                 double ahighmn[nbgpr] ={1.5*height, mean+rms, 1.5*rms};
01971                 
01972                 TMinuit *gMinuit = new TMinuit(nbgpr);
01973                 gMinuit->SetFCN(fcnbg);
01974                 
01975                 double arglist[10];
01976                 int ierflg = 0;
01977                 arglist[0] =0.5;
01978                 gMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
01979                 char name[100];
01980                 for (int k=0; k<nbgpr; k++) {
01981                   sprintf(name, "pedpar%i",k);
01982                   gMinuit->mnparm(k, name, strt[k], step[k], alowmn[k], ahighmn[k],ierflg);
01983                 }
01984                 
01985                 arglist[0] = 0;
01986                 gMinuit->mnexcm("SIMPLEX", arglist, 0, ierflg);
01987                 
01988                 arglist[0] = 0;
01989                 gMinuit->mnexcm("IMPROVE", arglist, 0, ierflg);
01990                 
01991                 TString chnam;
01992                 double parv,err,xlo,xup, plerr, mierr, eparab, gcc;
01993                 int iuit;
01994                 
01995                 for (int k=0; k<nbgpr; k++) {
01996                   if (step[k] >-10) {
01997                     gMinuit->mnpout(k, chnam, parv, err, xlo, xup, iuit);
01998                     gMinuit->mnerrs(k, plerr, mierr, eparab, gcc);
01999                     //              cout <<"k "<< k<<" "<<chnam<<" "<<parv<<" "<<err<<" "<<xlo<<" "<<xup<<" "<<plerr<<" "<<mierr<<" "<<eparab<<endl;
02000                     if (k==0) {
02001                       gaupr[k] = parv*binwid;
02002                       parer[k] = err*binwid;
02003                     } else {
02004                       gaupr[k] = parv;
02005                       parer[k] = err;   
02006                     }
02007                   }
02008                 }
02009 
02010                 //              gx0[izone]->SetParameters(gaupr);
02011                 
02012                 char temp[20];
02013                 sprintf(temp, "ped0fun_%i",izone);
02014                 ped0fun[izone] = new TF1(temp, gausX, xmn, xmx, nbgpr);
02015                 ped0fun[izone]->SetParameters(gaupr);
02016                 ped0fun[izone]->SetLineColor(3);
02017                 ped0fun[izone]->SetLineWidth(1);
02018                 ped0fun[izone]->Draw("same");   
02019                 
02020                 delete  gMinuit;
02021               }
02022             } else {
02023               for (int k=0; k<nbgpr; k++) {gaupr[k] = par[k]; }
02024               if (m_digiInput) { gaupr[2] = 0.90; } else { gaupr[2] = 0.15;}
02025             }
02026           }
02027           //      if (iijj!=0) 
02028           c0->cd(2*izone+2);
02029           if (signall[izone]->GetEntries() >5) {
02030             Double_t parall[nsgpr];
02031             double parserr[nsgpr];
02032             double fitres[nsgpr];
02033             double pedht = 0;
02034             
02035             char temp[20];
02036             sprintf(temp, "signal_%i",izone);
02037             xmn = signall[izone]->GetXaxis()->GetXmin();
02038             xmx = 0.5*signall[izone]->GetXaxis()->GetXmax();
02039             signal[izone] = new TF1(temp, totalfunc, xmn, xmx, nsgpr);
02040             xmx *=2.0;
02041             if ((iijj!=3) || m_histfit) {
02042               pedht = (signall[izone]->GetBinContent(nbn-1)+
02043                        signall[izone]->GetBinContent(nbn)+
02044                        signall[izone]->GetBinContent(nbn+1))/3.;
02045               
02046               if (m_pedsuppr && !m_digiInput) {
02047                 parall[1] = 0.0; // pedmean[ietafit][iphifit];
02048                 if (m_digiInput) { parall[2] = 0.90; } else { parall[2] = 0.15;}
02049               } else {
02050                 for (int i=0; i<nbgpr; i++) {parall[i] = gaupr[i];}
02051               }
02052               
02053               set_mean(parall[1], m_digiInput);
02054               set_sigma(parall[2], m_digiInput);
02055 
02056               parall[0] = 0.9*pedht; //GM for Z-mumu, there is almost no pedestal
02057               parall[3] = 0.14;
02058               double area = binwid*signall[izone]->GetEntries();
02059               parall[5]= area;
02060               
02061               if (iijj==3) {
02062                 parall[4] = fitprm[4][jk];
02063                 parall[6] = fitprm[6][jk];
02064               } else {
02065                 parall[4] = signall[izone]->GetMean();
02066                 parall[6]=parall[2];
02067               }
02068               
02069               signal[izone]->SetParameters(parall);
02070               signal[izone]->FixParameter(1, parall[1]);
02071               signal[izone]->FixParameter(2, parall[2]); 
02072               signal[izone]->SetParLimits(0, 0.00, 2.0*pedht+0.1);
02073               signal[izone]->FixParameter(3, 0.14);
02074               
02075               signal[izone]->SetParLimits(5, 0.40*area, 1.15*area);
02076               //              if (m_histfit) { //GMA
02077               if (iijj==3) {
02078                 signal[izone]->SetParLimits(4, 0.2*fitprm[4][jk], 2.0*fitprm[4][jk]);
02079                 signal[izone]->SetParLimits(6, 0.2*fitprm[6][jk], 2.0*fitprm[6][jk]);
02080               } else {
02081                 if (m_digiInput) {
02082                   signal[izone]->SetParLimits(4, 0.6, 6.0);
02083                   signal[izone]->SetParLimits(6, 0.60, 3.0);
02084                 } else {
02085                   signal[izone]->SetParLimits(4, 0.1, 1.0); 
02086                   signal[izone]->SetParLimits(6, 0.035, 0.3);
02087                 }
02088               }
02089               signal[izone]->SetParNames("const", "mean", "sigma","Width","MP","Area","GSigma");   
02090               signall[izone]->Fit(signal[izone], "0R+");
02091               
02092               signall[izone]->GetXaxis()->SetRangeUser(xmn,xmx);
02093               for (int k=0; k<nsgpr; k++) {
02094                 fitres[k] = fitprm[k][jk] = signal[izone]->GetParameter(k);
02095                 parserr[k] = signal[izone]->GetParError(k);
02096               }
02097               
02098             } else {
02099               double pedhtx = 0;
02100               for (unsigned i =0; i<sig_reg[ietafit][iphifit].size(); i++) {
02101                 if (sig_reg[ietafit][iphifit][i] >gaupr[1]-3*gaupr[2] && sig_reg[ietafit][iphifit][i]<gaupr[1]+gaupr[2]) pedhtx++;
02102               }
02103               
02104               set_mean(gaupr[1], m_digiInput);
02105               set_sigma(gaupr[2], m_digiInput);
02106 
02107               TString name[nsgpr] = {"const", "mean", "sigma","Width","MP","Area","GSigma"};
02108               double strt[nsgpr] = {0.9*pedhtx, gaupr[1], gaupr[2], fitprm[3][jk], fitprm[4][jk], signall[izone]->GetEntries(), fitprm[6][jk]};
02109               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]};
02110               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]};
02111               double step[nsgpr] = {1.0, 0.0, 0.0, 0.0, 0.001, 1.0, 0.002};
02112               
02113               TMinuit *gMinuit = new TMinuit(nsgpr);
02114               gMinuit->SetFCN(fcnsg);
02115               
02116               double arglist[10];
02117               int ierflg = 0;
02118               arglist[0] =0.5;
02119               gMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
02120               
02121               for (int k=0; k<nsgpr; k++) {
02122                 gMinuit->mnparm(k, name[k], strt[k], step[k], alowmn[k], ahighmn[k],ierflg);
02123               }
02124               
02125               arglist[0] = 0;
02126               gMinuit->mnexcm("SIMPLEX", arglist, 0, ierflg);
02127               
02128               arglist[0] = 0;
02129               gMinuit->mnexcm("IMPROVE", arglist, 0, ierflg);
02130               
02131               TString chnam;
02132               double parv,err,xlo,xup, plerr, mierr, eparab, gcc;
02133               int iuit;
02134               
02135               for (int k=0; k<nsgpr; k++) {
02136                 if (step[k] >-10) {
02137                   gMinuit->mnpout(k, chnam, parv, err, xlo, xup, iuit);
02138                   gMinuit->mnerrs(k, plerr, mierr, eparab, gcc);
02139                   if (k==0 || k==5) { 
02140                     fitres[k] = parv*binwid;
02141                     parserr[k]= err*binwid;
02142                   } else {
02143                     fitres[k] = parv;
02144                     parserr[k]= err;
02145                   }
02146                   
02147                 }
02148               }
02149               
02150               delete gMinuit;
02151             }
02152 
02153             //      if (iijj==0) {
02154             //        signall[izone]->Draw("same");
02155             //      } else {
02156               signall[izone]->Draw();
02157               //            }
02158             
02159             sprintf(temp, "pedfun_%i",izone);
02160             pedfun[izone] = new TF1(temp, gausX, xmn, xmx, nbgpr);
02161             pedfun[izone]->SetParameters(fitres);
02162             pedfun[izone]->SetLineColor(3);
02163             pedfun[izone]->SetLineWidth(1);
02164             pedfun[izone]->Draw("same");        
02165             
02166             sprintf(temp, "signalfun_%i",izone);
02167             sigfun[izone] = new TF1(temp, langaufun, xmn, xmx, nsgpr-nbgpr);
02168             sigfun[izone]->SetParameters(&fitres[3]);
02169             sigfun[izone]->SetLineWidth(1);
02170             sigfun[izone]->SetLineColor(4);
02171             sigfun[izone]->Draw("same");
02172             
02173             sprintf(temp, "total_%i",izone);
02174             signalx[izone] = new TF1(temp, totalfunc, xmn, xmx, nsgpr);
02175             signalx[izone]->SetParameters(fitres);
02176             signalx[izone]->SetLineWidth(1);
02177             signalx[izone]->Draw("same");
02178             
02179             int kl = (jk<15) ? jk+1 : 14-jk;
02180 
02181             cout<<"histinfo"<<iijj<<" fit "
02182                 <<std::setw(3)<< kl<<" "
02183                 <<std::setw(3)<< ij+1<<" "
02184                 <<std::setw(5)<<pedstll[izone]->GetEntries()<<" "
02185                 <<std::setw(6)<<pedstll[izone]->GetMean()<<" "
02186                 <<std::setw(6)<<pedstll[izone]->GetRMS()<<" "
02187                 <<std::setw(5)<<signall[izone]->GetEntries()<<" "
02188                 <<std::setw(6)<<signall[izone]->GetMean()<<" "
02189                 <<std::setw(6)<<signall[izone]->GetRMS()<<" "
02190                 <<std::setw(6)<< signal[izone]->GetChisquare()<<" "
02191                 <<std::setw(3)<< signal[izone]->GetNDF()<<endl;
02192             
02193             file_out<<"histinfo"<<iijj<<" fit "
02194                     <<std::setw(3)<< kl<<" "
02195                     <<std::setw(3)<< ij+1<<" "
02196                     <<std::setw(5)<<pedstll[izone]->GetEntries()<<" "
02197                     <<std::setw(6)<<pedstll[izone]->GetMean()<<" "
02198                     <<std::setw(6)<<pedstll[izone]->GetRMS()<<" "
02199                     <<std::setw(5)<<signall[izone]->GetEntries()<<" "
02200                     <<std::setw(6)<<signall[izone]->GetMean()<<" "
02201                     <<std::setw(6)<<signall[izone]->GetRMS()<<" "
02202                     <<std::setw(6)<< signal[izone]->GetChisquare()<<" "
02203                     <<std::setw(3)<< signal[izone]->GetNDF()<<endl;
02204             
02205             file_out <<"fitres x"<<iijj<<" "<<kl<<" "<<ij+1<<" "<< fitres[0]<<" "<< fitres[1]<<" "<< fitres[2]<<" "<< fitres[3]<<" "<< fitres[4]<<" "<< fitres[5]<<" "<< fitres[6]<<endl;
02206             file_out <<"parserr"<<iijj<<" "<<kl<<" "<<ij+1<<" "<< parserr[0]<<" "<< parserr[1]<<" "<< parserr[2]<<" "<< parserr[3]<<" "<< parserr[4]<<" "<< parserr[5]<<" "<< parserr[6]<<endl;    
02207 
02208             double diff=fitres[4]-fitres[1];
02209             if (diff <=0) diff = 0.000001;
02210             double error=parserr[4]*parserr[4]+parer[2]*parer[2];
02211             error = pow(error,0.5);
02212             
02213             int ieta = (jk<15) ? (15+jk) : (29-jk);
02214             int ifl = nphimx*ieta + ij;
02215             
02216             if (iijj==3) {
02217               ped_evt->Fill(ifl,pedstll[izone]->GetEntries());
02218               ped_mean->Fill(ifl,gaupr[1]);
02219               ped_width->Fill(ifl,gaupr[2]);
02220               fit_chi->Fill(ifl,signal[izone]->GetChisquare());
02221               sig_evt->Fill(ifl, signall[izone]->GetEntries());
02222               fit_sigevt->Fill(ifl, fitres[5]);
02223               fit_bkgevt->Fill(ifl, fitres[0]*sqrt(2*acos(-1.))*gaupr[2]);
02224               sig_mean->Fill(ifl, fitres[4]);
02225               sig_diff->Fill(ifl, fitres[4]-fitres[1]);
02226               sig_width->Fill(ifl, fitres[3]);
02227               sig_sigma->Fill(ifl, fitres[6]);
02228               sig_meanerr->Fill(ifl, parserr[4]);
02229               if (fitres[4]-fitres[1] !=0) sig_meanerrp->Fill(ifl, 100*parserr[4]/(fitres[4]-fitres[1]));
02230               if (gaupr[2]!=0) sig_signf->Fill(ifl,(fitres[4]-fitres[1])/gaupr[2]); 
02231               
02232               ped_statmean->Fill(ifl,pedstll[izone]->GetMean());
02233               sig_statmean->Fill(ifl,signall[izone]->GetMean());
02234               ped_rms->Fill(ifl,pedstll[izone]->GetRMS());
02235               sig_rms->Fill(ifl,signall[izone]->GetRMS());
02236             }
02237             
02238             if ((iijj==2) || (iijj==3) || (iijj==1)) {
02239               if (signall[izone]->GetEntries() >5 && fitres[4]>0.1) {
02240                 //GMA need to put this==1 in future
02241                 float fact=0.812;
02242                 if (abs(kl)<=4) fact=0.895;
02243                 if (!m_digiInput) fact *=0.19; //conversion factor for GeV/fC
02244 
02245                 float fact2 = 0;
02246                 if (iijj==2) fact2 = invang[jk][nphimx];
02247                 if (iijj==3) fact2 = invang[jk][ij];
02248                 if (iijj==1) fact2 = com_invang[jk][ij];
02249 
02250                 float calibc = fact*fact2/(fitres[4]*signall[izone]->GetEntries());
02251                 float caliberr= TMath::Abs(calibc*parserr[4]/max(0.001,fitres[4]));
02252 
02253                 if (iijj==2) {
02254                   int ieta = (jk<15) ? jk+1 : 14-jk;
02255                   mean_phi_hst->Fill(ieta, calibc);
02256                   mean_phi_hst->SetBinError(mean_phi_hst->FindBin(ieta), caliberr);
02257                   file_out<<"intieta "<<jk<<" "<<ij<<" "<<ieta<<" "<<mean_phi_hst->FindBin(double(ieta))<<" "<<calibc<<" "<<caliberr<<endl;
02258                 } else if (iijj==3) {
02259                   const_eta[jk]->Fill(ij+1,calibc);
02260                   const_eta[jk]->SetBinError(const_eta[jk]->FindBin(ij+1), caliberr);
02261                   
02262                   peak_eta[jk]->Fill(ij+1,fitres[4]);
02263                   peak_eta[jk]->SetBinError(peak_eta[jk]->FindBin(ij+1),parserr[4]);
02264 
02265                   int ieta = (jk<15) ? jk+1 : 14-jk;
02266                   const_eta_phi->Fill(ieta, ij+1,calibc);
02267                   file_out<<"intietax "<<jk<<" "<<ij<<" "<<ieta<<" "<<const_eta_phi->FindBin(ieta, ij+1)<<endl;
02268                   if (caliberr >0) {
02269                     const_eta_phi->SetBinError(const_eta_phi->FindBin(ieta, ij+1),caliberr);
02270 
02271                     mean_eta[ij] +=calibc/(caliberr*caliberr);
02272                     mean_phi[jk] +=calibc/(caliberr*caliberr);
02273                     
02274                     rms_eta[ij] +=1./(caliberr*caliberr);
02275                     rms_phi[jk] +=1./(caliberr*caliberr);
02276 
02277                   } else {
02278                     const_eta_phi->SetBinError(const_eta_phi->FindBin(ieta, ij+1), 0.0);
02279                   }
02280                 } else if (iijj==1) {
02281                   const_hpdrm[jk]->Fill(ij+1,calibc);
02282                   const_hpdrm[jk]->SetBinError(const_hpdrm[jk]->FindBin(ij+1), caliberr);
02283                   
02284                   peak_hpdrm[jk]->Fill(ij+1,fitres[4]);
02285                   peak_hpdrm[jk]->SetBinError(peak_hpdrm[jk]->FindBin(ij+1),parserr[4]);
02286                 }
02287 
02288                 file_out<<"HO  4 "<<iijj<<" "<< std::setw(3)<<kl<<" "<<std::setw(3)<<ij+1<<" "
02289                         <<std::setw(7)<<calibc<<" "<<std::setw(7)<<caliberr<<endl;
02290               }
02291             }
02292             
02293           } else {   //if (signall[izone]->GetEntries() >10) {
02294             signall[izone]->Draw();
02295             float varx = 0.000;
02296             int kl = (jk<15) ? jk+1 : 14-jk;
02297             file_out<<"histinfo"<<iijj<<" nof "
02298                     <<std::setw(3)<< kl<<" "
02299                     <<std::setw(3)<< ij+1<<" "
02300                     <<std::setw(5)<<pedstll[izone]->GetEntries()<<" "
02301                     <<std::setw(6)<<pedstll[izone]->GetMean()<<" "
02302                     <<std::setw(6)<<pedstll[izone]->GetRMS()<<" "
02303                     <<std::setw(5)<<signall[izone]->GetEntries()<<" "
02304                     <<std::setw(6)<<signall[izone]->GetMean()<<" "
02305                     <<std::setw(6)<<signall[izone]->GetRMS()<<" "
02306                     <<std::setw(6)<< varx<<" "
02307                     <<std::setw(3)<< varx<<endl;
02308             
02309             file_out <<"fitres x"<<iijj<<" "<<kl<<" "<<ij+1<<" "<< varx<<" "<< varx<<" "<< varx<<" "<< varx<<" "<< varx<<" "<< varx<<" "<< varx<<endl;
02310             file_out <<"parserr"<<iijj<<" "<<kl<<" "<<ij+1<<" "<< varx<<" "<< varx<<" "<< varx<<" "<< varx<<" "<< varx<<" "<< varx<<" "<< varx<<endl;
02311             
02312           }
02313           iiter++;
02314           if (iiter%nsample==0) { 
02315             c0->Update();   
02316 
02317             for (int kl=0; kl<nsample; kl++) {
02318               if (gx0[kl]) {delete gx0[kl];gx0[kl] = 0;}
02319               if (ped0fun[kl]) {delete ped0fun[kl];ped0fun[kl] = 0;}
02320               if (signal[kl]) {delete signal[kl];signal[kl] = 0;}
02321               if (pedfun[kl]) {delete pedfun[kl];pedfun[kl] = 0;}
02322               if (sigfun[kl]) {delete sigfun[kl];sigfun[kl] = 0;}
02323               if (signalx[kl]) {delete signalx[kl];signalx[kl] = 0;}
02324               if (signall[kl]) {delete signall[kl];signall[kl] = 0;}
02325               if (pedstll[kl]) {delete pedstll[kl];pedstll[kl] = 0;}
02326             }
02327 
02328           }
02329         } //for (int jk=0; jk<netamx; jk++) { 
02330       } //for (int ij=0; ij<nphimx; ij++) {
02331       
02332       //      if (iijj==0) {
02333       //        sprintf(out_file, "comb_hosig_allring_%i.jpg", irunold);
02334       //        c0->SaveAs(out_file);
02335       //        iiter = 0;
02336       //      } else {
02337       //        //      c0->Update(); 
02338       //      }
02339 
02340       //      iiter = 0;
02341     } //end of iijj
02342     if (iiter%nsample!=0) { 
02343       c0->Update(); 
02344       for (int kl=0; kl<nsample; kl++) {
02345         if (gx0[kl]) {delete gx0[kl];gx0[kl] = 0;}
02346         if (ped0fun[kl]) {delete ped0fun[kl];ped0fun[kl] = 0;}
02347         if (signal[kl]) {delete signal[kl];signal[kl] = 0;}
02348         if (pedfun[kl]) {delete pedfun[kl];pedfun[kl] = 0;}
02349         if (sigfun[kl]) {delete sigfun[kl];sigfun[kl] = 0;}
02350         if (signalx[kl]) {delete signalx[kl];signalx[kl] = 0;}
02351         if (signall[kl]) {delete signall[kl];signall[kl] = 0;}
02352         if (pedstll[kl]) {delete pedstll[kl];pedstll[kl] = 0;}
02353       }
02354     }
02355 
02356     delete c0;
02357 
02358     xsiz = 600; //int xsiz = 600;
02359     ysiz = 800; //int ysiz = 800;
02360     
02361     gStyle->SetTitleFontSize(0.05);
02362     gStyle->SetTitleSize(0.025,"XYZ");
02363     gStyle->SetLabelSize(0.025,"XYZ");
02364     gStyle->SetStatFontSize(.045);
02365   
02366     gStyle->SetOptStat(0);
02367     ps.NewPage(); TCanvas *c1 = new TCanvas("c1", " Pedestal vs signal", xsiz, ysiz);
02368     ped_evt->Draw(); c1->Update();
02369     
02370     ps.NewPage();
02371     ped_statmean->Draw(); c1->Update();
02372     
02373     ps.NewPage();
02374     ped_rms->Draw(); c1->Update();
02375     
02376     ps.NewPage();
02377     ped_mean->Draw(); c1->Update();
02378     
02379     ps.NewPage();
02380     ped_width->Draw(); c1->Update();
02381     
02382     ps.NewPage();
02383     sig_evt->Draw(); c1->Update();
02384     
02385     ps.NewPage();
02386     sig_statmean->Draw(); c1->Update();
02387     
02388     ps.NewPage();
02389     sig_rms->Draw(); c1->Update();
02390     
02391     ps.NewPage();
02392     fit_chi->Draw(); c1->Update();
02393     
02394     ps.NewPage();
02395     fit_sigevt->Draw(); c1->Update();
02396     
02397     ps.NewPage();
02398     fit_bkgevt->Draw(); c1->Update();
02399     
02400     ps.NewPage();
02401     sig_mean->Draw(); c1->Update();
02402     
02403     ps.NewPage();
02404     sig_width->Draw(); c1->Update();
02405     
02406     ps.NewPage(); 
02407     sig_sigma->Draw(); c1->Update();
02408     
02409     ps.NewPage(); 
02410     sig_meanerr->Draw(); c1->Update();
02411     
02412     ps.NewPage(); 
02413     sig_meanerrp->Draw(); c1->Update();
02414     
02415     ps.NewPage(); 
02416     sig_signf->Draw(); c1->Update();
02417     
02418     ps.Close();
02419     delete c1;
02420     
02421     file_out.close();
02422 
02423     if (m_figure) {
02424       xsiz = 700;
02425       ysiz = 450;
02426 
02427       gStyle->SetTitleFontSize(0.09);
02428       gStyle->SetPadBottomMargin(0.17);
02429       gStyle->SetPadLeftMargin(0.18);
02430       gStyle->SetPadRightMargin(0.01);
02431       gStyle->SetOptLogy(0);
02432       gStyle->SetOptStat(0);
02433       
02434       TCanvas *c2 = new TCanvas("c2", "runfile", xsiz, ysiz); 
02435       c2->Divide(5,3);
02436       
02437       for (int side=0; side <2; side++) {
02438         gStyle->SetNdivisions(303,"XY");
02439         gStyle->SetPadRightMargin(0.01);
02440         int nmn = 0;
02441         int nmx = netamx/2;
02442         if (side==1) {
02443           nmn = netamx/2;
02444           nmx = netamx;
02445         }
02446         
02447         int nzone = 0;
02448         
02449         for (int ij=nmn; ij<nmx; ij++) {
02450           
02451           c2->cd(nzone+1); 
02452           const_eta[ij]->GetXaxis()->SetTitle("#phi index");
02453           const_eta[ij]->GetXaxis()->SetTitleSize(.08);
02454           const_eta[ij]->GetXaxis()->CenterTitle(); 
02455           const_eta[ij]->GetXaxis()->SetTitleOffset(0.9);
02456           const_eta[ij]->GetXaxis()->SetLabelSize(.085);
02457           const_eta[ij]->GetXaxis()->SetLabelOffset(.01);
02458           
02459           const_eta[ij]->GetYaxis()->SetLabelSize(.08);
02460           const_eta[ij]->GetYaxis()->SetLabelOffset(.01);       
02461           if (m_digiInput) {
02462             const_eta[ij]->GetYaxis()->SetTitle("GeV/fC");
02463           } else {
02464             const_eta[ij]->GetYaxis()->SetTitle("GeV/MIP-GeV!!");
02465           }
02466           
02467           const_eta[ij]->GetYaxis()->SetTitleSize(.085);
02468           const_eta[ij]->GetYaxis()->CenterTitle(); 
02469           const_eta[ij]->GetYaxis()->SetTitleOffset(1.3);
02470           const_eta[ij]->SetMarkerSize(0.60);
02471           const_eta[ij]->SetMarkerColor(2);
02472           const_eta[ij]->SetMarkerStyle(20);
02473 
02474 
02475           const_eta[ij]->Draw();
02476           nzone++;
02477         } 
02478         
02479         sprintf(out_file, "calibho_%i_side%i.eps", irunold, side);
02480         c2->SaveAs(out_file);
02481         
02482         sprintf(out_file, "calibho_%i_side%i.jpg", irunold, side);
02483         c2->SaveAs(out_file);
02484         
02485         nzone = 0;
02486         for (int ij=nmn; ij<nmx; ij++) {
02487           c2->cd(nzone+1); 
02488           peak_eta[ij]->GetXaxis()->SetTitle("#phi index");
02489           peak_eta[ij]->GetXaxis()->SetTitleSize(.08);
02490           peak_eta[ij]->GetXaxis()->CenterTitle(); 
02491           peak_eta[ij]->GetXaxis()->SetTitleOffset(0.90);
02492           peak_eta[ij]->GetXaxis()->SetLabelSize(.08);
02493           peak_eta[ij]->GetXaxis()->SetLabelOffset(.01);
02494           
02495           peak_eta[ij]->GetYaxis()->SetLabelSize(.08);
02496           peak_eta[ij]->GetYaxis()->SetLabelOffset(.01);        
02497           if (m_digiInput) {
02498             peak_eta[ij]->GetYaxis()->SetTitle("fC");
02499           } else {
02500             peak_eta[ij]->GetYaxis()->SetTitle("GeV");
02501           }
02502           
02503           peak_eta[ij]->GetYaxis()->SetTitleSize(.085);
02504           peak_eta[ij]->GetYaxis()->CenterTitle(); 
02505           peak_eta[ij]->GetYaxis()->SetTitleOffset(1.3);
02506           
02507           peak_eta[ij]->SetMarkerSize(0.60);
02508           peak_eta[ij]->SetMarkerColor(2);
02509           peak_eta[ij]->SetMarkerStyle(20);
02510 
02511           peak_eta[ij]->Draw();
02512           nzone++;
02513         } 
02514         
02515         sprintf(out_file, "peakho_%i_side%i.eps", irunold, side);
02516         c2->SaveAs(out_file);
02517         
02518         sprintf(out_file, "peakho_%i_side%i.jpg", irunold, side);
02519         c2->SaveAs(out_file);
02520       }
02521       delete c2;
02522 
02523       //      if (m_combined) {
02524       gStyle->SetTitleFontSize(0.045);
02525       gStyle->SetPadRightMargin(0.13);
02526       gStyle->SetPadBottomMargin(0.15);
02527       gStyle->SetPadLeftMargin(0.1);
02528       gStyle->SetOptStat(0);
02529       xsiz = 700;
02530       ysiz = 600;
02531       TCanvas *c1 = new TCanvas("c1", "Fitted const in each tower", xsiz, ysiz);  
02532       const_eta_phi->GetXaxis()->SetTitle("#eta");
02533       const_eta_phi->GetXaxis()->SetTitleSize(0.065);
02534       const_eta_phi->GetXaxis()->SetTitleOffset(0.85); //6); 
02535       const_eta_phi->GetXaxis()->CenterTitle();
02536       const_eta_phi->GetXaxis()->SetLabelSize(0.045);
02537       const_eta_phi->GetXaxis()->SetLabelOffset(0.01); 
02538       
02539       const_eta_phi->GetYaxis()->SetTitle("#phi");
02540       const_eta_phi->GetYaxis()->SetTitleSize(0.075);
02541       const_eta_phi->GetYaxis()->SetTitleOffset(0.5); 
02542       const_eta_phi->GetYaxis()->CenterTitle();
02543       const_eta_phi->GetYaxis()->SetLabelSize(0.045);
02544       const_eta_phi->GetYaxis()->SetLabelOffset(0.01); 
02545       
02546       const_eta_phi->Draw("colz");
02547       sprintf(out_file, "high_hoconst_eta_phi_%i.jpg",irunold); 
02548       c1->SaveAs(out_file); 
02549       
02550       delete c1; 
02551       
02552       for (int jk=0; jk<netamx; jk++) {
02553         int ieta = (jk<15) ? jk+1 : 14-jk;
02554         if (rms_phi[jk]>0) {
02555           mean_phi_ave->Fill(ieta, mean_phi[jk]/rms_phi[jk]);
02556           mean_phi_ave->SetBinError(mean_phi_ave->FindBin(ieta), pow(double(rms_phi[jk]), -0.5));
02557         }
02558       }
02559       
02560       for (int ij=0; ij<nphimx; ij++) {
02561         if (rms_eta[ij] >0) {
02562           mean_eta_ave->Fill(ij+1, mean_eta[ij]/rms_eta[ij]);
02563           mean_eta_ave->SetBinError(mean_eta_ave->FindBin(ij+1), pow(double(rms_eta[ij]), -0.5));
02564         }
02565       }
02566       
02567       ysiz =450;
02568       gStyle->SetPadLeftMargin(0.13);
02569       gStyle->SetPadRightMargin(0.03);
02570 
02571 
02572       TCanvas *c2y = new TCanvas("c2", "Avearge signal in eta and phi", xsiz, ysiz);
02573       c2y->Divide(2,1);
02574       mean_eta_ave->GetXaxis()->SetTitle("#phi");
02575       mean_eta_ave->GetXaxis()->SetTitleSize(0.085);
02576       mean_eta_ave->GetXaxis()->SetTitleOffset(0.65); 
02577       mean_eta_ave->GetXaxis()->CenterTitle();
02578       mean_eta_ave->GetXaxis()->SetLabelSize(0.05);
02579       mean_eta_ave->GetXaxis()->SetLabelOffset(0.001); 
02580       
02581       mean_eta_ave->GetYaxis()->SetTitle("Signal (GeV)/MIP");
02582       mean_eta_ave->GetYaxis()->SetTitleSize(0.055);
02583       mean_eta_ave->GetYaxis()->SetTitleOffset(1.3); 
02584       mean_eta_ave->GetYaxis()->CenterTitle();
02585       mean_eta_ave->GetYaxis()->SetLabelSize(0.045);
02586       mean_eta_ave->GetYaxis()->SetLabelOffset(0.01); 
02587       mean_eta_ave->SetMarkerSize(0.60);
02588       mean_eta_ave->SetMarkerColor(2);
02589       mean_eta_ave->SetMarkerStyle(20);
02590 
02591       c2y->cd(1); mean_eta_ave->Draw();
02592       
02593       mean_phi_ave->GetXaxis()->SetTitle("#eta");
02594       mean_phi_ave->GetXaxis()->SetTitleSize(0.085);
02595       mean_phi_ave->GetXaxis()->SetTitleOffset(0.65); //55); 
02596       mean_phi_ave->GetXaxis()->CenterTitle();
02597       mean_phi_ave->GetXaxis()->SetLabelSize(0.05);
02598       mean_phi_ave->GetXaxis()->SetLabelOffset(0.001); 
02599       
02600       mean_phi_ave->GetYaxis()->SetTitle("Signal (GeV)/MIP");
02601       mean_phi_ave->GetYaxis()->SetTitleSize(0.055);
02602       mean_phi_ave->GetYaxis()->SetTitleOffset(1.3); 
02603       mean_phi_ave->GetYaxis()->CenterTitle();
02604       mean_phi_ave->GetYaxis()->SetLabelSize(0.045);
02605       mean_phi_ave->GetYaxis()->SetLabelOffset(0.01); 
02606       mean_phi_ave->SetMarkerSize(0.60);
02607       mean_phi_ave->SetMarkerColor(2);
02608       mean_phi_ave->SetMarkerStyle(20);
02609       
02610       c2y->cd(2); mean_phi_ave->Draw();
02611       
02612       sprintf(out_file, "high_hoaverage_eta_phi_%i.jpg",irunold); 
02613       c2y->SaveAs(out_file); 
02614       
02615       delete c2y;
02616       //      } else { //m_combined
02617 
02618       xsiz = 800;
02619       ysiz = 450;
02620       TCanvas *c3 = new TCanvas("c3", "Avearge signal in eta and phi", xsiz, ysiz);
02621       c3->Divide(2,1);
02622       mean_phi_hst->GetXaxis()->SetTitle("#eta");
02623       mean_phi_hst->GetXaxis()->SetTitleSize(0.065);
02624       mean_phi_hst->GetXaxis()->SetTitleOffset(0.9); 
02625       mean_phi_hst->GetXaxis()->CenterTitle();
02626       mean_phi_hst->GetXaxis()->SetLabelSize(0.065);
02627       mean_phi_hst->GetXaxis()->SetLabelOffset(0.001); 
02628       
02629       mean_phi_hst->GetYaxis()->SetTitle("GeV/MIP");
02630       mean_phi_hst->GetYaxis()->SetTitleSize(0.055);
02631       mean_phi_hst->GetYaxis()->SetTitleOffset(0.9); 
02632       mean_phi_hst->GetYaxis()->CenterTitle();
02633       mean_phi_hst->GetYaxis()->SetLabelSize(0.065);
02634       mean_phi_hst->GetYaxis()->SetLabelOffset(0.01); 
02635       
02636       mean_phi_hst->SetMarkerColor(4);
02637       mean_phi_hst->SetMarkerSize(0.8);
02638       mean_phi_hst->SetMarkerStyle(20);
02639       mean_phi_hst->Draw();
02640       
02641       sprintf(out_file, "low_mean_phi_hst_%i.jpg",irunold); 
02642       c3->SaveAs(out_file); 
02643       
02644       delete c3;
02645       
02646       //      } //m_combined
02647 
02648 
02649       gStyle->SetOptLogy(1);
02650       gStyle->SetPadTopMargin(.1);
02651       gStyle->SetPadLeftMargin(.15);
02652       xsiz = 800;
02653       ysiz = 500;
02654       TCanvas *c0x = new TCanvas("c0x", "Signal in each ring", xsiz, ysiz);
02655 
02656       c0x->Divide(3,2);
02657       for (int ij=0; ij<ringmx; ij++) {
02658         int iread = (ij==2) ? routmx : rout12mx;
02659         if (m_digiInput) {
02660           com_sigrsg[ij][iread]->GetXaxis()->SetTitle("Signal/ped (fC)");
02661         } else {
02662           com_sigrsg[ij][iread]->GetXaxis()->SetTitle("Signal/ped (GeV)");
02663         }
02664         com_sigrsg[ij][iread]->GetXaxis()->SetTitleSize(0.060);
02665         com_sigrsg[ij][iread]->GetXaxis()->SetTitleOffset(1.05);
02666         com_sigrsg[ij][iread]->GetXaxis()->CenterTitle();
02667         com_sigrsg[ij][iread]->GetXaxis()->SetLabelSize(0.065);
02668         com_sigrsg[ij][iread]->GetXaxis()->SetLabelOffset(0.01);
02669         
02670         com_sigrsg[ij][iread]->GetYaxis()->SetLabelSize(0.065);
02671         com_sigrsg[ij][iread]->GetYaxis()->SetLabelOffset(0.01); 
02672 
02673 
02674         com_sigrsg[ij][iread]->SetLineWidth(3);  
02675         com_sigrsg[ij][iread]->SetLineColor(4);  
02676         
02677         c0x->cd(ij+1); com_sigrsg[ij][iread]->Draw();
02678         
02679         com_crossg[ij][iread]->SetLineWidth(2);  
02680         com_crossg[ij][iread]->SetLineColor(2);
02681         com_crossg[ij][iread]->Draw("same");
02682       }
02683       sprintf(out_file, "hosig_ring_%i.jpg",irunold); 
02684       c0x->SaveAs(out_file);
02685       delete c0x;
02686 
02687       gStyle->SetTitleFontSize(0.06);
02688       gStyle->SetOptStat(0);
02689       gStyle->SetOptLogy(0);
02690 
02691       TCanvas *c0 = new TCanvas("c0", "Signal in each ring", xsiz, ysiz);
02692       
02693       c0->Divide(3,2);
02694       for (int jk=0; jk<ringmx; jk++) {
02695         peak_hpdrm[jk]->GetXaxis()->SetTitle("RM #");
02696         peak_hpdrm[jk]->GetXaxis()->SetTitleSize(0.070);
02697         peak_hpdrm[jk]->GetXaxis()->SetTitleOffset(1.0); 
02698         peak_hpdrm[jk]->GetXaxis()->CenterTitle();
02699         peak_hpdrm[jk]->GetXaxis()->SetLabelSize(0.065);
02700         peak_hpdrm[jk]->GetXaxis()->SetLabelOffset(0.01);
02701 
02702         peak_hpdrm[jk]->GetYaxis()->SetTitle("Peak(GeV)/MIP");
02703 
02704         peak_hpdrm[jk]->GetYaxis()->SetTitleSize(0.07);
02705         peak_hpdrm[jk]->GetYaxis()->SetTitleOffset(1.3); 
02706         peak_hpdrm[jk]->GetYaxis()->CenterTitle();
02707         peak_hpdrm[jk]->GetYaxis()->SetLabelSize(0.065);
02708         peak_hpdrm[jk]->GetYaxis()->SetLabelOffset(0.01); 
02709         //      peak_hpdrm[jk]->SetLineWidth(3);  
02710         //      peak_hpdrm[jk]->SetLineColor(4);  
02711         peak_hpdrm[jk]->SetMarkerSize(0.60);
02712         peak_hpdrm[jk]->SetMarkerColor(2);
02713         peak_hpdrm[jk]->SetMarkerStyle(20);
02714 
02715         
02716         c0->cd(jk+1); peak_hpdrm[jk]->Draw();
02717       }
02718       sprintf(out_file, "comb_peak_hpdrm_%i.jpg",irunold); 
02719       c0->SaveAs(out_file);
02720       
02721       delete c0;
02722 
02723       TCanvas *c1y = new TCanvas("c1y", "Signal in each ring", xsiz, ysiz);
02724       
02725       c1y->Divide(3,2);
02726       for (int jk=0; jk<ringmx; jk++) {
02727         const_hpdrm[jk]->GetXaxis()->SetTitle("RM #");
02728         const_hpdrm[jk]->GetXaxis()->SetTitleSize(0.070);
02729         const_hpdrm[jk]->GetXaxis()->SetTitleOffset(1.3); 
02730         const_hpdrm[jk]->GetXaxis()->CenterTitle();
02731         const_hpdrm[jk]->GetXaxis()->SetLabelSize(0.065);
02732         const_hpdrm[jk]->GetXaxis()->SetLabelOffset(0.01);
02733 
02734         if (m_digiInput) {
02735           const_hpdrm[jk]->GetYaxis()->SetTitle("Peak(fC)");
02736         } else {
02737           const_hpdrm[jk]->GetYaxis()->SetTitle("Peak(GeV)");
02738         }
02739         const_hpdrm[jk]->GetYaxis()->SetTitleSize(0.065);
02740         const_hpdrm[jk]->GetYaxis()->SetTitleOffset(1.0); 
02741         const_hpdrm[jk]->GetYaxis()->CenterTitle();
02742         const_hpdrm[jk]->GetYaxis()->SetLabelSize(0.065);
02743         const_hpdrm[jk]->GetYaxis()->SetLabelOffset(0.01); 
02744         //      const_hpdrm[jk]->SetLineWidth(3);  
02745         //      const_hpdrm[jk]->SetLineColor(4);  
02746         const_hpdrm[jk]->SetMarkerSize(0.60);
02747         const_hpdrm[jk]->SetMarkerColor(2);
02748         const_hpdrm[jk]->SetMarkerStyle(20);
02749 
02750         c1y->cd(jk+1); const_hpdrm[jk]->Draw();
02751       }
02752 
02753       sprintf(out_file, "comb_const_hpdrm_%i.jpg",irunold); 
02754       c1y->SaveAs(out_file);
02755       
02756       delete c1y;
02757 
02758     } //if (m_figure) {
02759 
02760     //    ps.Close();
02761     //    file_out.close();
02762 
02763   }// if (m_constant){ 
02764 
02765 
02766   if (m_figure) {
02767     for (int ij=0; ij<nphimx; ij++) {
02768       for (int jk=0; jk<netamx; jk++) {
02769         stat_eta[jk]->Fill(ij+1,sigrsg[jk][ij]->GetEntries());
02770         statmn_eta[jk]->Fill(ij+1,sigrsg[jk][ij]->GetMean());
02771       }
02772     }
02773     
02774     xsiz = 700;
02775     ysiz = 450;
02776     gStyle->SetTitleFontSize(0.09);
02777     gStyle->SetPadBottomMargin(0.14);
02778     gStyle->SetPadLeftMargin(0.17);
02779     gStyle->SetPadRightMargin(0.01);
02780     gStyle->SetNdivisions(303,"XY");
02781     gStyle->SetOptLogy(1);
02782     
02783     TCanvas *c2x = new TCanvas("c2x", "runfile", xsiz, ysiz); 
02784     c2x->Divide(5,3);
02785     for (int side=0; side <2; side++) {
02786       int nmn = 0;
02787       int nmx = netamx/2;
02788       if (side==1) {
02789         nmn = netamx/2;
02790         nmx = netamx;
02791       }
02792       int nzone = 0;
02793       char name[200];
02794 
02795       for (int ij=nmn; ij<nmx; ij++) {
02796         int ieta = (ij<15) ? ij+1 : 14-ij;
02797         c2x->cd(nzone+1); 
02798         if (m_digiInput) {
02799           sprintf(name,"fC(#eta=%i)",ieta);
02800         } else {
02801           sprintf(name,"GeV(#eta=%i)",ieta);
02802         }
02803         sigrsg[ij][nphimx]->GetXaxis()->SetTitle(name);
02804         sigrsg[ij][nphimx]->GetXaxis()->SetTitleSize(.08);
02805         sigrsg[ij][nphimx]->GetXaxis()->CenterTitle(); 
02806         sigrsg[ij][nphimx]->GetXaxis()->SetTitleOffset(0.90);
02807         sigrsg[ij][nphimx]->GetXaxis()->SetLabelSize(.08);
02808         sigrsg[ij][nphimx]->GetXaxis()->SetLabelOffset(.01);
02809         
02810         sigrsg[ij][nphimx]->GetYaxis()->SetLabelSize(.08);
02811         sigrsg[ij][nphimx]->GetYaxis()->SetLabelOffset(.01);    
02812         sigrsg[ij][nphimx]->SetLineWidth(2);
02813         sigrsg[ij][nphimx]->SetLineColor(4);
02814         sigrsg[ij][nphimx]->Draw();
02815         crossg[ij][nphimx]->SetLineWidth(2);
02816         crossg[ij][nphimx]->SetLineColor(2);
02817         crossg[ij][nphimx]->Draw("same");
02818         nzone++;
02819       } 
02820       
02821       sprintf(out_file, "sig_ho_%i_side%i.eps", irunold, side);
02822       c2x->SaveAs(out_file);
02823       
02824       sprintf(out_file, "sig_ho_%i_side%i.jpg", irunold, side);
02825       c2x->SaveAs(out_file);
02826     }
02827 
02828     gStyle->SetOptLogy(0);
02829     c2x = new TCanvas("c2x", "runfile", xsiz, ysiz); 
02830     c2x->Divide(5,3);
02831     for (int side=0; side <2; side++) {
02832       int nmn = 0;
02833       int nmx = netamx/2;
02834       if (side==1) {
02835         nmn = netamx/2;
02836         nmx = netamx;
02837       }
02838       int nzone = 0;
02839 
02840       nzone = 0;
02841       for (int ij=nmn; ij<nmx; ij++) {
02842         c2x->cd(nzone+1); 
02843         statmn_eta[ij]->SetLineWidth(2);  
02844         statmn_eta[ij]->SetLineColor(4);  
02845         statmn_eta[ij]->GetXaxis()->SetTitle("#phi index");     
02846         statmn_eta[ij]->GetXaxis()->SetTitleSize(.08);
02847         statmn_eta[ij]->GetXaxis()->CenterTitle(); 
02848         statmn_eta[ij]->GetXaxis()->SetTitleOffset(0.9);
02849         statmn_eta[ij]->GetYaxis()->SetLabelSize(.08);
02850         statmn_eta[ij]->GetYaxis()->SetLabelOffset(.01);        
02851         statmn_eta[ij]->GetXaxis()->SetLabelSize(.08);
02852         statmn_eta[ij]->GetXaxis()->SetLabelOffset(.01);
02853         if (m_digiInput) {
02854           statmn_eta[ij]->GetYaxis()->SetTitle("fC");
02855         } else {
02856           statmn_eta[ij]->GetYaxis()->SetTitle("GeV");
02857         }
02858         statmn_eta[ij]->GetYaxis()->SetTitleSize(.075);
02859         statmn_eta[ij]->GetYaxis()->CenterTitle(); 
02860         statmn_eta[ij]->GetYaxis()->SetTitleOffset(1.30);
02861         
02862         statmn_eta[ij]->Draw();
02863         nzone++;
02864       } 
02865       
02866       sprintf(out_file, "statmnho_%i_side%i.eps", irunold, side);
02867       c2x->SaveAs(out_file);
02868       
02869       sprintf(out_file, "statmnho_%i_side%i.jpg", irunold, side);
02870       c2x->SaveAs(out_file);
02871       
02872       gStyle->SetOptLogy(1);
02873       gStyle->SetNdivisions(203,"XY");
02874       
02875       nzone = 0;
02876       for (int ij=nmn; ij<nmx; ij++) {
02877         c2x->cd(nzone+1); 
02878         stat_eta[ij]->SetLineWidth(2);  
02879         stat_eta[ij]->SetLineColor(4);  
02880         stat_eta[ij]->GetXaxis()->SetTitle("#phi index");       
02881         stat_eta[ij]->GetXaxis()->SetTitleSize(.08);
02882         stat_eta[ij]->GetXaxis()->CenterTitle(); 
02883         stat_eta[ij]->GetXaxis()->SetTitleOffset(0.80);
02884         stat_eta[ij]->GetXaxis()->SetLabelSize(.08);
02885         stat_eta[ij]->GetXaxis()->SetLabelOffset(.01);  
02886         stat_eta[ij]->GetYaxis()->SetLabelSize(.08);
02887         stat_eta[ij]->GetYaxis()->SetLabelOffset(.01);
02888         
02889         stat_eta[ij]->Draw();
02890         nzone++;
02891       } 
02892       
02893       sprintf(out_file, "statho_%i_side%i.eps", irunold, side);
02894       c2x->SaveAs(out_file);
02895       
02896       sprintf(out_file, "statho_%i_side%i.jpg", irunold, side);
02897       c2x->SaveAs(out_file);
02898     }
02899     delete c2x;
02900 
02901   } //if (m_figure) {
02902 
02903   if (!m_constant) { //m_constant
02904     for (int j=0; j<netamx; j++) {
02905       for (int i=0; i<nphimx; i++) {
02906         if (crossg[j][i]) { delete crossg[j][i];}
02907         if (sigrsg[j][i]) { delete sigrsg[j][i];}
02908       }
02909     }
02910   }
02911 }
02912 
02913 //define this as a plug-in
02914 DEFINE_FWK_MODULE(HOCalibAnalyzer);
02915 
02916 
02917