Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-02-18 08:23:09

0001 //April 2015 : Removal of itrg1, itrg2, but addition of isect2, same is true in HOCalibVariables.h
0002 // -*- C++ -*-
0003 //
0004 // Package:    HOCalibAnalyzer
0005 // Class:      HOCalibAnalyzer
0006 //
0007 /**\class HOCalibAnalyzer HOCalibAnalyzer.cc Calibration/HOCalibAnalyzer/src/HOCalibAnalyzer.cc
0008 
0009  Description: <one line class summary>
0010 
0011  Implementation:
0012      <Notes on implementation>
0013 
0014 April 2015
0015 // Addition of these variables, ilumi (analyser), inslumi (analyser), nprim
0016 
0017 
0018 */
0019 //
0020 // Original Author:  Gobinda Majumder
0021 //         Created:  Sat Jul  7 09:51:31 CEST 2007
0022 //
0023 //
0024 
0025 // system include files
0026 #include <memory>
0027 
0028 // user include files
0029 #include "FWCore/Framework/interface/Frameworkfwd.h"
0030 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0031 
0032 #include "FWCore/Framework/interface/Event.h"
0033 #include "FWCore/Framework/interface/MakerMacros.h"
0034 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0035 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0036 
0037 #include "FWCore/ServiceRegistry/interface/Service.h"
0038 #include "FWCore/Utilities/interface/InputTag.h"
0039 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0040 #include "DataFormats/HcalCalibObjects/interface/HOCalibVariables.h"
0041 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0042 
0043 #include "TFile.h"
0044 #include "TH1F.h"
0045 #include "TH2F.h"
0046 #include "TTree.h"
0047 #include "TProfile.h"
0048 
0049 #include <string>
0050 
0051 #include <iostream>
0052 #include <fstream>
0053 #include <iomanip>
0054 
0055 using namespace std;
0056 using namespace edm;
0057 
0058 //#define EDM_ML_DEBUG
0059 //
0060 // class decleration
0061 //
0062 
0063 class HOCalibAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0064 public:
0065   explicit HOCalibAnalyzer(const edm::ParameterSet&);
0066   ~HOCalibAnalyzer() override;
0067 
0068 private:
0069   void beginJob() override {}
0070   void analyze(const edm::Event&, const edm::EventSetup&) override;
0071   void endJob() override {}
0072 
0073   static constexpr int netamx = 30;
0074   static constexpr int nphimx = 72;
0075   static constexpr int ringmx = 5;
0076   static constexpr int ncut = 14;
0077 
0078   const char* varcrit[3] = {"All", "steps", "n-1"};  // or opposite
0079 
0080   const double elosfact = (14.9 + 0.96 * fabs(log(8 * 2.8)) + 0.033 * 8 * (1.0 - pow(8, -0.33)));
0081 
0082   int getHOieta(int ij) { return (ij < netamx / 2) ? -netamx / 2 + ij : -netamx / 2 + ij + 1; }
0083   int invert_HOieta(int ieta) { return (ieta < 0) ? netamx / 2 + ieta : netamx / 2 + ieta - 1; }
0084 
0085   int mypow_2[31];
0086 
0087   const bool m_cosmic;
0088   const bool m_zeroField;
0089   const int m_bins;
0090   const double m_low;
0091   const double m_ahigh;
0092   const bool m_histFill;  //Stored signals of individual HO towers with default selection criteria
0093   const bool m_treeFill;  //Store rootuple without almost any selection criteria except a quality on muon
0094 
0095   int ipass;
0096 
0097   TTree* T1;
0098 
0099   TH1F* ho_indenergy[netamx][nphimx];
0100 
0101   TH1F* muonnm;
0102   TH1F* muonmm;
0103   TH1F* muonth;
0104   TH1F* muonph;
0105   TH1F* muonch;
0106 
0107   TH1F* sel_muonnm;
0108   TH1F* sel_muonmm;
0109   TH1F* sel_muonth;
0110   TH1F* sel_muonph;
0111   TH1F* sel_muonch;
0112 
0113   //  TProfile* sigvsevt[15][ncut];
0114   TH2F* sig_eta_evt[3 * netamx][ncut];  //For individual eta
0115   TH2F* sigvsevt[3 * netamx][ncut];
0116   TH1F* variab[3 * netamx][ncut];
0117 
0118   TH2F* mu_projection[ncut + 1];
0119 
0120   unsigned ievt, hoflag;
0121   int irun, ilumi, nprim, isect, isect2, ndof, nmuon;
0122 
0123   float inslumi, trkdr, trkdz, trkvx, trkvy, trkvz, trkmm, trkth, trkph, chisq, therr, pherr, hodx, hody, hoang, htime,
0124       hosig[9], hocorsig[18], hocro, hbhesig[9], caloen[3];
0125   float momatho, tkpt03, ecal03, hcal03;
0126   float tmphoang;
0127 
0128   int nevents[10];
0129 
0130   float ncount[ringmx][ncut + 10];
0131 
0132   edm::InputTag hoCalibVariableCollectionTag;
0133   const edm::EDGetTokenT<HOCalibVariableCollection> tok_ho_;
0134   const edm::EDGetTokenT<HORecHitCollection> tok_allho_;
0135   // ----------member data ---------------------------
0136 };
0137 
0138 //
0139 // constants, enums and typedefs
0140 //
0141 
0142 //
0143 // static data member definitions
0144 //
0145 
0146 //
0147 // constructors and destructor
0148 //
0149 
0150 HOCalibAnalyzer::HOCalibAnalyzer(const edm::ParameterSet& iConfig)
0151     : m_cosmic(iConfig.getUntrackedParameter<bool>("cosmic", true)),
0152       m_zeroField(iConfig.getUntrackedParameter<bool>("zeroField", false)),
0153       m_bins(iConfig.getUntrackedParameter<int>("HOSignalBins", 120)),
0154       m_low(iConfig.getUntrackedParameter<double>("lowerRange", -1.0)),
0155       m_ahigh(iConfig.getUntrackedParameter<double>("upperRange", 29.0)),
0156       m_histFill(iConfig.getUntrackedParameter<bool>("histFill", true)),
0157       m_treeFill(iConfig.getUntrackedParameter<bool>("treeFill", false)),
0158       tok_ho_(consumes<HOCalibVariableCollection>(iConfig.getParameter<edm::InputTag>("hoCalibVariableCollectionTag"))),
0159       tok_allho_(consumes<HORecHitCollection>(iConfig.getParameter<edm::InputTag>("hoInputTag"))) {
0160   // It is very likely you want the following in your configuration
0161   // hoCalibVariableCollectionTag = cms.InputTag('hoCalibProducer', 'HOCalibVariableCollection')
0162 
0163   usesResource(TFileService::kSharedResource);
0164 
0165   //now do what ever initialization is needed
0166   ipass = 0;
0167   for (int ij = 0; ij < 10; ij++) {
0168     nevents[ij] = 0;
0169   }
0170 
0171   edm::Service<TFileService> fs;
0172 
0173   T1 = fs->make<TTree>("T1", "HOSignal");
0174 
0175   T1->Branch("irun", &irun, "irun/I");
0176   T1->Branch("ievt", &ievt, "ievt/i");
0177 
0178   T1->Branch("isect", &isect, "isect/I");
0179   T1->Branch("isect2", &isect2, "isect2/I");
0180   T1->Branch("ndof", &ndof, "ndof/I");
0181   T1->Branch("nmuon", &nmuon, "nmuon/I");
0182 
0183   T1->Branch("ilumi", &ilumi, "ilumi/I");
0184   if (!m_cosmic) {
0185     T1->Branch("inslumi", &inslumi, "inslumi/F");
0186     T1->Branch("nprim", &nprim, "nprim/I");
0187     T1->Branch("tkpt03", &tkpt03, " tkpt03/F");
0188     T1->Branch("ecal03", &ecal03, " ecal03/F");
0189     T1->Branch("hcal03", &hcal03, " hcal03/F");
0190   }
0191 
0192   T1->Branch("trkdr", &trkdr, "trkdr/F");
0193   T1->Branch("trkdz", &trkdz, "trkdz/F");
0194 
0195   T1->Branch("trkvx", &trkvx, "trkvx/F");
0196   T1->Branch("trkvy", &trkvy, "trkvy/F");
0197   T1->Branch("trkvz", &trkvz, "trkvz/F");
0198   T1->Branch("trkmm", &trkmm, "trkmm/F");
0199   T1->Branch("trkth", &trkth, "trkth/F");
0200   T1->Branch("trkph", &trkph, "trkph/F");
0201 
0202   T1->Branch("chisq", &chisq, "chisq/F");
0203   T1->Branch("therr", &therr, "therr/F");
0204   T1->Branch("pherr", &pherr, "pherr/F");
0205   T1->Branch("hodx", &hodx, "hodx/F");
0206   T1->Branch("hody", &hody, "hody/F");
0207   T1->Branch("hoang", &hoang, "hoang/F");
0208 
0209   T1->Branch("momatho", &momatho, "momatho/F");
0210   T1->Branch("hoflag", &hoflag, "hoflag/i");
0211   T1->Branch("htime", &htime, "htime/F");
0212   T1->Branch("hosig", hosig, "hosig[9]/F");
0213   T1->Branch("hocro", &hocro, "hocro/F");
0214   T1->Branch("hocorsig", hocorsig, "hocorsig[18]/F");
0215   T1->Branch("caloen", caloen, "caloen[3]/F");
0216 
0217   char name[200];
0218   char title[200];
0219 
0220   if (m_histFill) {
0221     for (int ij = 0; ij < netamx; ij++) {
0222       int ieta = getHOieta(ij);
0223       for (int jk = 0; jk < nphimx; jk++) {
0224         sprintf(name, "ho_indenergy_%i_%i", ij, jk);
0225         sprintf(title, "ho IndEnergy (GeV) i#eta=%i i#phi=%i", ieta, jk + 1);
0226         ho_indenergy[ij][jk] = fs->make<TH1F>(name, title, 1200, m_low, m_ahigh);
0227       }
0228     }
0229   }
0230 
0231   muonnm = fs->make<TH1F>("muonnm", "No of muon", 10, -0.5, 9.5);
0232   muonmm = fs->make<TH1F>("muonmm", "P_{mu}", 200, -100., 100.);
0233   muonth = fs->make<TH1F>("muonth", "{Theta}_{mu}", 180, 0., 180.);
0234   muonph = fs->make<TH1F>("muonph", "{Phi}_{mu}", 180, -180., 180.);
0235   muonch = fs->make<TH1F>("muonch", "{chi^2}/ndf", 100, 0., 1000.);
0236 
0237   sel_muonnm = fs->make<TH1F>("sel_muonnm", "No of muon(sel)", 10, -0.5, 9.5);
0238   sel_muonmm = fs->make<TH1F>("sel_muonmm", "P_{mu}(sel)", 200, -100., 100.);
0239   sel_muonth = fs->make<TH1F>("sel_muonth", "{Theta}_{mu}(sel)", 180, 0., 180.);
0240   sel_muonph = fs->make<TH1F>("sel_muonph", "{Phi}_{mu}(sel)", 180, -180., 180.);
0241   sel_muonch = fs->make<TH1F>("sel_muonch", "{chi^2}/ndf(sel)", 100, 0., 1000.);
0242 
0243   float pival = acos(-1.);
0244 
0245   //if change order, change in iselect_wotime also and other efficiency numbers
0246   const char* varnam[ncut] = {"ndof",
0247                               "chisq",
0248                               "th",
0249                               "ph",
0250                               "therr",
0251                               "pherr",
0252                               "dircos",
0253                               "trkmm",
0254                               "nmuon",
0255                               "calo",
0256                               "trkiso",
0257                               "#phi-dir",
0258                               "#eta-dir",
0259                               "time"};
0260   int nbinxx[ncut] = {25, 60, 60, 60, 60, 60, 60, 120, 6, 60, 60, 120, 120, 60};
0261   double alowxx[ncut] = {5.5, 0., 0., -pival, 0.0, 0.0, 0.0, 0., 0.5, 0.0, 0.0, -20., -32., -45.0};
0262   double ahghxx[ncut] = {30.5, 40., pival, pival, 0.8, 0.02, 0.5, 300., 6.5, 10.0, 24.0, 20.0, 32.0, 45.0};
0263 
0264   for (int kl = 0; kl < ncut; kl++) {
0265     for (int jk = 0; jk < 3; jk++) {
0266       for (int ij = 0; ij < netamx; ij++) {
0267         sprintf(name, "sigeta_%i_%i_%i", kl, jk, ij);
0268         sprintf(title, "sigeta %s %s i#eta=%i", varnam[kl], varcrit[jk], getHOieta(ij));
0269         sig_eta_evt[netamx * jk + ij][kl] =
0270             fs->make<TH2F>(name, title, nbinxx[kl], alowxx[kl], ahghxx[kl], m_bins, m_low, m_ahigh);
0271       }
0272     }
0273   }
0274 
0275   for (int kl = 0; kl < ncut; kl++) {
0276     for (int ij = 0; ij < ringmx * 3; ij++) {
0277       int iring = ij % ringmx - 2;
0278       int iset = ij / ringmx;
0279       sprintf(name, "sigring_%i_%i", kl, ij);
0280       sprintf(title, "Signal %s %s Ring%i", varnam[kl], varcrit[iset], iring);
0281       sigvsevt[ij][kl] = fs->make<TH2F>(name, title, nbinxx[kl], alowxx[kl], ahghxx[kl], m_bins, m_low, m_ahigh);
0282     }
0283   }
0284 
0285   for (int kl = 0; kl < ncut; kl++) {
0286     for (int ij = 0; ij < ringmx * 3; ij++) {
0287       int iring = ij % ringmx - 2;
0288       int iset = ij / ringmx;
0289       sprintf(name, "varring_%i_%i", kl, ij);
0290       sprintf(title, "%s %s Ring%i", varnam[kl], varcrit[iset], iring);
0291       variab[ij][kl] = fs->make<TH1F>(name, title, nbinxx[kl], alowxx[kl], ahghxx[kl]);
0292     }
0293   }
0294 
0295   for (int ij = 0; ij <= ncut; ij++) {
0296     sprintf(name, "mu_projection_%i", ij);
0297     if (ij == 0) {
0298       sprintf(title, "All projected muon");
0299     } else {
0300       sprintf(title, "Projected muon with selection %s", varnam[ij - 1]);
0301     }
0302     mu_projection[ij] =
0303         fs->make<TH2F>(name, title, netamx + 1, -netamx / 2 - 0.5, netamx / 2 + 0.5, nphimx, 0.5, nphimx + 0.5);
0304   }
0305 
0306   for (int ij = 0; ij < 31; ij++) {
0307     mypow_2[ij] = pow(2, ij);
0308   }
0309   for (int ij = 0; ij < ringmx; ij++) {
0310     for (int jk = 0; jk < ncut + 10; jk++) {
0311       ncount[ij][jk] = 0.0;
0312     }
0313   }
0314 }
0315 
0316 HOCalibAnalyzer::~HOCalibAnalyzer() {
0317   // do anything here that needs to be done at desctruction time
0318   // (e.g. close files, deallocate resources etc.)
0319 
0320   edm::LogVerbatim("HOCalibAnalyzer") << " Total events = " << setw(7) << nevents[0] << " " << setw(7) << nevents[1]
0321                                       << " " << setw(7) << nevents[2] << " " << setw(7) << nevents[3] << " " << setw(7)
0322                                       << nevents[4] << " " << setw(7) << nevents[5] << " Selected events # is "
0323                                       << ipass;
0324 }
0325 
0326 //
0327 // member functions
0328 //
0329 
0330 // ------------ method called to for each event  ------------
0331 void HOCalibAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0332   nevents[0]++;
0333 
0334   using namespace edm;
0335 
0336   float pival = acos(-1.);
0337 
0338   ievt = iEvent.id().event();
0339   ilumi = iEvent.luminosityBlock();
0340 
0341   const edm::Handle<HOCalibVariableCollection>& HOCalib = iEvent.getHandle(tok_ho_);
0342 
0343   if (nevents[0] % 20000 == 1) {
0344     edm::LogVerbatim("HOCalibAnalyzer") << "nmuon event # " << setw(7) << nevents[0] << " " << setw(7) << nevents[1]
0345                                         << " " << setw(7) << nevents[2] << " " << setw(7) << nevents[3] << " "
0346                                         << setw(7) << nevents[4] << " " << setw(7) << nevents[5];
0347     edm::LogVerbatim("HOCalibAnalyzer") << " Run # " << iEvent.id().run() << " Evt # " << iEvent.id().event() << " "
0348                                         << int(HOCalib.isValid()) << " " << ipass;
0349   }
0350 
0351   if (HOCalib.isValid()) {
0352     nevents[1]++;
0353     nmuon = (*HOCalib).size();
0354 
0355     for (HOCalibVariableCollection::const_iterator hoC = (*HOCalib).begin(); hoC != (*HOCalib).end(); hoC++) {
0356       trkdr = (*hoC).trkdr;
0357       trkdz = (*hoC).trkdz;
0358 
0359       trkvx = (*hoC).trkvx;
0360       trkvy = (*hoC).trkvy;
0361       trkvz = (*hoC).trkvz;
0362 
0363       trkmm = (*hoC).trkmm;
0364       trkth = (*hoC).trkth;
0365       trkph = (*hoC).trkph;
0366 
0367       ndof = (int)(*hoC).ndof;
0368       chisq = (*hoC).chisq;
0369       momatho = (*hoC).momatho;
0370 
0371       therr = (*hoC).therr;
0372       pherr = (*hoC).pherr;
0373       trkph = (*hoC).trkph;
0374 
0375       if (!m_cosmic) {
0376         nprim = (*hoC).nprim;
0377         inslumi = (*hoC).inslumi;
0378         tkpt03 = (*hoC).tkpt03;
0379         ecal03 = (*hoC).ecal03;
0380         hcal03 = (*hoC).hcal03;
0381       }
0382 
0383       isect = (*hoC).isect;
0384       isect2 = (*hoC).isect2;
0385       hodx = (*hoC).hodx;
0386       hody = (*hoC).hody;
0387       hoang = (*hoC).hoang;
0388 
0389       tmphoang = sin(trkth) - hoang;
0390 
0391       htime = (*hoC).htime;
0392       hoflag = (*hoC).hoflag;
0393       for (int ij = 0; ij < 9; ij++) {
0394         hosig[ij] = (*hoC).hosig[ij];
0395 #ifdef EDM_ML_DEBUG
0396         edm::LogVerbatim("HOCalibAnalyzer") << "hosig " << ij << " " << hosig[ij];
0397 #endif
0398       }
0399       for (int ij = 0; ij < 18; ij++) {
0400         hocorsig[ij] = (*hoC).hocorsig[ij];
0401 #ifdef EDM_ML_DEBUG
0402         edm::LogVerbatim("HOCalibAnalyzer") << "hocorsig " << ij << " " << hocorsig[ij];
0403 #endif
0404       }
0405       hocro = (*hoC).hocro;
0406       for (int ij = 0; ij < 3; ij++) {
0407         caloen[ij] = (*hoC).caloen[ij];
0408       }
0409 
0410       int ipsall = 0;
0411       int ips0 = 0;
0412       int ips1 = 0;
0413       int ips2 = 0;
0414       int ips3 = 0;
0415       int ips4 = 0;
0416       int ips5 = 0;
0417       int ips6 = 0;
0418       int ips7 = 0;
0419       int ips8 = 0;
0420       int ips9 = 0;
0421       int ips10 = 0;
0422       int ips11 = 0;
0423       int ips12 = 0;
0424       int ips13 = 0;
0425 
0426       nevents[2]++;
0427       bool isZSps = (hosig[4] < -99.0) ? false : true;
0428 
0429       if ((!m_cosmic) && fabs(trkmm) < momatho)
0430         continue;
0431 
0432       nevents[3]++;
0433       if (fabs(trkth - pival / 2) < 0.000001)
0434         continue;  //22OCT07
0435       nevents[4]++;
0436 
0437       int ieta = int((abs(isect) % 10000) / 100.) - 50;  //an offset to acodate -ve eta values
0438       if (abs(ieta) >= 16)
0439         continue;
0440       nevents[5]++;
0441       int iphi = abs(isect) % 100;
0442 
0443       int iring = 0;
0444 
0445       int iring2 = iring + 2;
0446 
0447       double abshoang = (m_cosmic) ? fabs(hoang) : hoang;
0448 
0449       double elos = 1.0 / TMath::Max(0.1, abs(1.0 * hoang));
0450 
0451       if (!m_zeroField)
0452         elos *= ((14.9 + 0.96 * fabs(log(momatho * 2.8)) + 0.033 * momatho * (1.0 - pow(momatho, -0.33))) / elosfact);
0453 
0454       if (m_cosmic) {
0455         if (abs(ndof) >= 20 && abs(ndof) < 55) {
0456           ips0 = mypow_2[0];
0457           ipsall += ips0;
0458         }
0459         if (chisq > 0 && chisq < 12) {
0460           ips1 = mypow_2[1];
0461           ipsall += ips1;
0462         }  //18Jan2008
0463 
0464         if (trkth > 0.3 && trkth < pival - 0.3) {
0465           ips2 = mypow_2[2];
0466           ipsall += ips2;
0467         }  //No nead for pp evt
0468         if (trkph > -pival + 0.1 && trkph < -0.1) {
0469           ips3 = mypow_2[3];
0470           ipsall += ips3;
0471         }  //No nead for pp evt
0472 
0473         if (therr < 0.02) {
0474           ips4 = mypow_2[4];
0475           ipsall += ips4;
0476         }
0477         if (pherr < 0.0002) {
0478           ips5 = mypow_2[5];
0479           ipsall += ips5;
0480         }
0481         if (abshoang > 0.60 && abshoang < 1.0) {
0482           ips6 = mypow_2[6];
0483           ipsall += ips6;
0484         }
0485 
0486         if (m_zeroField || (fabs(momatho) > 5.0 && fabs(momatho) < 2000.0)) {
0487           ips7 = mypow_2[7];
0488           ipsall += ips7;
0489         }
0490 
0491         if (nmuon >= 1 && nmuon <= 3) {
0492           ips8 = mypow_2[8];
0493           ipsall += ips8;
0494         }
0495 
0496         //  if (hodx>0 && hody>0) { }
0497         ips9 = mypow_2[9];
0498         ipsall += ips9;
0499 
0500         ips10 = mypow_2[10];
0501         ipsall += ips10;
0502 
0503         if (iring2 == 2) {
0504           if (fabs(hodx) < 100 && fabs(hodx) > 2 && fabs(hocorsig[8]) < 40 && fabs(hocorsig[8]) > 2) {
0505             ips11 = mypow_2[11];
0506             ipsall += ips11;
0507           }
0508 
0509           if (fabs(hody) < 100 && fabs(hody) > 2 && fabs(hocorsig[9]) < 40 && fabs(hocorsig[9]) > 2) {
0510             ips12 = mypow_2[12];
0511             ipsall += ips12;
0512           }
0513 
0514         } else {
0515           if (fabs(hodx) < 100 && fabs(hodx) > 2) {
0516             ips11 = mypow_2[11];
0517             ipsall += ips11;
0518           }
0519 
0520           if (fabs(hody) < 100 && fabs(hody) > 2) {
0521             ips12 = mypow_2[12];
0522             ipsall += ips12;
0523           }
0524         }
0525 
0526         if (m_zeroField) {
0527           if (iring2 == 0) {
0528             if (htime > -60 && htime < 60) {
0529               ips13 = mypow_2[13];
0530               ipsall += ips13;
0531             }
0532           }
0533           if (iring2 == 1) {
0534             if (htime > -60 && htime < 60) {
0535               ips13 = mypow_2[13];
0536               ipsall += ips13;
0537             }
0538           }
0539           if (iring2 == 2) {
0540             if (htime > -60 && htime < 60) {
0541               ips13 = mypow_2[13];
0542               ipsall += ips13;
0543             }
0544           }
0545           if (iring2 == 3) {
0546             if (htime > -60 && htime < 60) {
0547               ips13 = mypow_2[13];
0548               ipsall += ips13;
0549             }
0550           }
0551           if (iring2 == 4) {
0552             if (htime > -60 && htime < 60) {
0553               ips13 = mypow_2[13];
0554               ipsall += ips13;
0555             }
0556           }
0557         } else {
0558           if (htime > -100 && htime < 100) {
0559             ips13 = mypow_2[13];
0560             ipsall += ips13;
0561           }
0562         }
0563       } else {
0564         if (abs(ndof) >= 10 && abs(ndof) < 25) {
0565           ips0 = mypow_2[0];
0566           ipsall += ips0;
0567         }
0568         if (chisq > 0 && chisq < 10) {
0569           ips1 = mypow_2[1];
0570           ipsall += ips1;
0571         }  //18Jan2008
0572 
0573         if (fabs(trkth - pival / 2) < 21.5) {
0574           ips2 = mypow_2[2];
0575           ipsall += ips2;
0576         }  //No nead for pp evt
0577         if (fabs(trkph + pival / 2) < 21.5) {
0578           ips3 = mypow_2[3];
0579           ipsall += ips3;
0580         }  //No nead for pp evt
0581 
0582         if (therr < 0.00002) {
0583           ips4 = mypow_2[4];
0584           ipsall += ips4;
0585         }
0586         if (pherr < 0.000002) {
0587           ips5 = mypow_2[5];
0588           ipsall += ips5;
0589         }
0590         //  if (abshoang >0.40 && abshoang <1.0) {ips6 = mypow_2[6];  ipsall +=ips6;}
0591         if (tmphoang < 0.065) {
0592           ips6 = mypow_2[6];
0593           ipsall += ips6;
0594         }
0595 
0596         if (fabs(momatho) < 250.0 && fabs(momatho) > 15.0) {
0597           if (iring2 == 2) {
0598             ips7 = mypow_2[7];
0599             ipsall += ips7;
0600           }
0601           if ((iring2 == 1 || iring2 == 3) && fabs(momatho) > 17.0) {
0602             ips7 = mypow_2[7];
0603             ipsall += ips7;
0604           }
0605           if ((iring2 == 0 || iring2 == 4) && fabs(momatho) > 20.0) {
0606             ips7 = mypow_2[7];
0607             ipsall += ips7;
0608           }
0609         }
0610 
0611         if (nmuon >= 1 && nmuon <= 3) {
0612           ips8 = mypow_2[8];
0613           ipsall += ips8;
0614         }
0615 
0616         if (ndof > 0 && caloen[0] < 15.0) {
0617           ips9 = mypow_2[9];
0618           ipsall += ips9;
0619         }  //5.0
0620         if (tkpt03 < 5.0) {
0621           ips10 = mypow_2[10];
0622           ipsall += ips10;
0623         }  //4.0
0624 
0625         if (iring2 == 2) {
0626           if (fabs(hodx) < 100 && fabs(hodx) > 2 && fabs(hocorsig[8]) < 40 && fabs(hocorsig[8]) > 2) {
0627             ips11 = mypow_2[11];
0628             ipsall += ips11;
0629           }
0630 
0631           if (fabs(hody) < 100 && fabs(hody) > 2 && fabs(hocorsig[9]) < 40 && fabs(hocorsig[9]) > 2) {
0632             ips12 = mypow_2[12];
0633             ipsall += ips12;
0634           }
0635 
0636         } else {
0637           if (fabs(hodx) < 100 && fabs(hodx) > 2) {
0638             ips11 = mypow_2[11];
0639             ipsall += ips11;
0640           }
0641 
0642           if (fabs(hody) < 100 && fabs(hody) > 2) {
0643             ips12 = mypow_2[12];
0644             ipsall += ips12;
0645           }
0646         }
0647 
0648         if (iring2 == 0) {
0649           if (htime > -20 && htime < 20) {
0650             ips13 = mypow_2[13];
0651             ipsall += ips13;
0652           }
0653         }
0654         if (iring2 == 1) {
0655           if (htime > -20 && htime < 20) {
0656             ips13 = mypow_2[13];
0657             ipsall += ips13;
0658           }
0659         }
0660         if (iring2 == 2) {
0661           if (htime > -30 && htime < 20) {
0662             ips13 = mypow_2[13];
0663             ipsall += ips13;
0664           }
0665         }
0666         if (iring2 == 3) {
0667           if (htime > -20 && htime < 20) {
0668             ips13 = mypow_2[13];
0669             ipsall += ips13;
0670           }
0671         }
0672         if (iring2 == 4) {
0673           if (htime > -20 && htime < 20) {
0674             ips13 = mypow_2[13];
0675             ipsall += ips13;
0676           }
0677         }
0678       }
0679       int tmpxet = invert_HOieta(ieta);
0680       double nomHOSig = hosig[4] / elos;
0681 
0682       if (ipsall - ips0 == mypow_2[ncut] - mypow_2[0] - 1) {
0683         if (isZSps) {
0684           sigvsevt[iring2][0]->Fill(abs(ndof), nomHOSig);
0685           sig_eta_evt[tmpxet][0]->Fill(abs(ndof), nomHOSig);
0686         }
0687         variab[iring2][0]->Fill(abs(ndof));
0688       }
0689       if (ipsall - ips1 == mypow_2[ncut] - mypow_2[1] - 1) {
0690         if (isZSps) {
0691           sigvsevt[iring2][1]->Fill(chisq, nomHOSig);
0692           sig_eta_evt[tmpxet][1]->Fill(chisq, nomHOSig);
0693         }
0694         variab[iring2][1]->Fill(chisq);
0695       }
0696       if (ipsall - ips2 == mypow_2[ncut] - mypow_2[2] - 1) {
0697         if (isZSps) {
0698           sigvsevt[iring2][2]->Fill(trkth, nomHOSig);
0699           sig_eta_evt[tmpxet][2]->Fill(trkth, nomHOSig);
0700         }
0701         variab[iring2][2]->Fill(trkth);
0702       }
0703       if (ipsall - ips3 == mypow_2[ncut] - mypow_2[3] - 1) {
0704         if (isZSps) {
0705           sigvsevt[iring2][3]->Fill(trkph, nomHOSig);
0706           sig_eta_evt[tmpxet][3]->Fill(trkph, nomHOSig);
0707         }
0708         variab[iring2][3]->Fill(trkph);
0709       }
0710       if (ipsall - ips4 == mypow_2[ncut] - mypow_2[4] - 1) {
0711         if (isZSps) {
0712           sigvsevt[iring2][4]->Fill(1000 * therr, nomHOSig);
0713           sig_eta_evt[tmpxet][4]->Fill(1000 * therr, nomHOSig);
0714         }
0715         variab[iring2][4]->Fill(1000 * therr);
0716       }
0717       if (ipsall - ips5 == mypow_2[ncut] - mypow_2[5] - 1) {
0718         if (isZSps) {
0719           sigvsevt[iring2][5]->Fill(1000 * pherr, nomHOSig);
0720           sig_eta_evt[tmpxet][5]->Fill(1000 * pherr, nomHOSig);
0721         }
0722         variab[iring2][5]->Fill(1000 * pherr);
0723       }
0724       if (ipsall - ips6 == mypow_2[ncut] - mypow_2[6] - 1) {
0725         if (isZSps) {
0726           sigvsevt[iring2][6]->Fill(tmphoang, (nomHOSig)*abshoang);
0727           sig_eta_evt[tmpxet][6]->Fill(tmphoang, (nomHOSig)*abshoang);
0728         }
0729         variab[iring2][6]->Fill(tmphoang);
0730       }
0731       if (ipsall - ips7 == mypow_2[ncut] - mypow_2[7] - 1) {
0732         if (isZSps) {
0733           sigvsevt[iring2][7]->Fill(fabs(trkmm), nomHOSig);
0734           sig_eta_evt[tmpxet][7]->Fill(fabs(trkmm), nomHOSig);
0735         }
0736         variab[iring2][7]->Fill(fabs(trkmm));
0737       }
0738       if (ipsall - ips8 == mypow_2[ncut] - mypow_2[8] - 1) {
0739         if (isZSps) {
0740           sigvsevt[iring2][8]->Fill(nmuon, nomHOSig);
0741           sig_eta_evt[tmpxet][8]->Fill(nmuon, nomHOSig);
0742         }
0743         variab[iring2][8]->Fill(nmuon);
0744       }
0745       if (!m_cosmic) {
0746         if (ipsall - ips9 == mypow_2[ncut] - mypow_2[9] - 1) {
0747           if (isZSps) {
0748             sigvsevt[iring2][9]->Fill(caloen[0], nomHOSig);
0749             sig_eta_evt[tmpxet][9]->Fill(caloen[0], nomHOSig);
0750           }
0751           variab[iring2][9]->Fill(caloen[0]);
0752         }
0753       }
0754       if (ipsall - ips10 == mypow_2[ncut] - mypow_2[10] - 1) {
0755         if (isZSps) {
0756           sigvsevt[iring2][10]->Fill(tkpt03, nomHOSig);
0757           sig_eta_evt[tmpxet][10]->Fill(tkpt03, nomHOSig);
0758         }
0759         variab[iring2][10]->Fill(tkpt03);
0760       }
0761       if (ipsall - ips11 == mypow_2[ncut] - mypow_2[11] - 1) {
0762         if (isZSps) {
0763           sigvsevt[iring2][11]->Fill(hodx, nomHOSig);
0764           sig_eta_evt[tmpxet][11]->Fill(hodx, nomHOSig);
0765         }
0766         variab[iring2][11]->Fill(hodx);
0767       }
0768       if (ipsall - ips12 == mypow_2[ncut] - mypow_2[12] - 1) {
0769         if (isZSps) {
0770           sigvsevt[iring2][12]->Fill(hody, nomHOSig);
0771           sig_eta_evt[tmpxet][12]->Fill(hody, nomHOSig);
0772         }
0773         variab[iring2][12]->Fill(hody);
0774       }
0775 
0776       if (ipsall - ips13 == mypow_2[ncut] - mypow_2[13] - 1) {
0777         if (isZSps) {
0778           sigvsevt[iring2][13]->Fill(htime, nomHOSig);
0779           sig_eta_evt[tmpxet][13]->Fill(htime, nomHOSig);
0780         }
0781         variab[iring2][13]->Fill(htime);
0782       }
0783 
0784       if (isZSps) {
0785         sigvsevt[iring2 + ringmx][0]->Fill(abs(ndof), nomHOSig);
0786         sig_eta_evt[netamx + tmpxet][0]->Fill(abs(ndof), nomHOSig);
0787       }
0788       variab[iring2 + 5][0]->Fill(abs(ndof));
0789 
0790       ncount[iring2][0]++;
0791       if (isZSps) {
0792         ncount[iring2][1]++;
0793       }
0794       if (ips0 > 0) {
0795         if (isZSps) {
0796           ncount[iring2][10]++;
0797           sigvsevt[iring2 + ringmx][1]->Fill(chisq, nomHOSig);
0798           sig_eta_evt[netamx + tmpxet][1]->Fill(chisq, nomHOSig);
0799         }
0800         variab[iring2 + ringmx][1]->Fill(chisq);
0801         mu_projection[1]->Fill(ieta, iphi);
0802         if (ips1 > 0) {
0803           if (isZSps) {
0804             ncount[iring2][11]++;
0805             sigvsevt[iring2 + ringmx][2]->Fill(trkth, nomHOSig);
0806             sig_eta_evt[netamx + tmpxet][2]->Fill(trkth, nomHOSig);
0807           }
0808           variab[iring2 + ringmx][2]->Fill(trkth);
0809           mu_projection[2]->Fill(ieta, iphi);
0810           if (ips2 > 0) {
0811             if (isZSps) {
0812               ncount[iring2][12]++;
0813               sigvsevt[iring2 + ringmx][3]->Fill(trkph, nomHOSig);
0814               sig_eta_evt[netamx + tmpxet][3]->Fill(trkph, nomHOSig);
0815             }
0816             variab[iring2 + ringmx][3]->Fill(trkph);
0817             mu_projection[3]->Fill(ieta, iphi);
0818             if (ips3 > 0) {
0819               if (isZSps) {
0820                 ncount[iring2][13]++;
0821                 sigvsevt[iring2 + ringmx][4]->Fill(1000 * therr, nomHOSig);
0822                 sig_eta_evt[netamx + tmpxet][4]->Fill(1000 * therr, nomHOSig);
0823               }
0824               variab[iring2 + ringmx][4]->Fill(1000 * therr);
0825               mu_projection[4]->Fill(ieta, iphi);
0826               if (ips4 > 0) {
0827                 if (isZSps) {
0828                   ncount[iring2][14]++;
0829                   sigvsevt[iring2 + ringmx][5]->Fill(1000 * pherr, nomHOSig);
0830                   sig_eta_evt[netamx + tmpxet][5]->Fill(1000 * pherr, nomHOSig);
0831                 }
0832                 variab[iring2 + ringmx][5]->Fill(1000 * pherr);
0833                 mu_projection[5]->Fill(ieta, iphi);
0834                 if (ips5 > 0) {
0835                   if (isZSps) {
0836                     ncount[iring2][15]++;
0837                     sigvsevt[iring2 + ringmx][6]->Fill(tmphoang, (nomHOSig)*abshoang);
0838                     sig_eta_evt[netamx + tmpxet][6]->Fill(tmphoang, (nomHOSig)*abshoang);
0839                   }
0840                   variab[iring2 + ringmx][6]->Fill(tmphoang);
0841                   mu_projection[6]->Fill(ieta, iphi);
0842                   if (ips6 > 0) {
0843                     if (isZSps) {
0844                       ncount[iring2][16]++;
0845                       sigvsevt[iring2 + ringmx][7]->Fill(fabs(trkmm), nomHOSig);
0846                       sig_eta_evt[netamx + tmpxet][7]->Fill(fabs(trkmm), nomHOSig);
0847                     }
0848                     variab[iring2 + ringmx][7]->Fill(fabs(trkmm));
0849                     mu_projection[7]->Fill(ieta, iphi);
0850                     if (ips7 > 0) {
0851                       ncount[iring2][4]++;  //Efficiency of Muon detection
0852                       if (isZSps) {
0853                         ncount[iring2][17]++;
0854                         sigvsevt[iring2 + ringmx][8]->Fill(nmuon, nomHOSig);
0855                         sig_eta_evt[netamx + tmpxet][8]->Fill(nmuon, nomHOSig);
0856                       }
0857                       variab[iring2 + ringmx][8]->Fill(nmuon);
0858                       mu_projection[8]->Fill(ieta, iphi);
0859                       if (ips8 > 0) {
0860                         if (!m_cosmic) {
0861                           if (isZSps) {
0862                             ncount[iring2][18]++;
0863                             sigvsevt[iring2 + ringmx][9]->Fill(caloen[0], nomHOSig);
0864                             sig_eta_evt[netamx + tmpxet][9]->Fill(caloen[0], nomHOSig);
0865                           }
0866                           variab[iring2 + ringmx][9]->Fill(caloen[0]);
0867                           mu_projection[9]->Fill(ieta, iphi);
0868                         }
0869                         if (ips9 > 0) {
0870                           if (isZSps) {
0871                             ncount[iring2][19]++;
0872                             sigvsevt[iring2 + ringmx][10]->Fill(tkpt03, nomHOSig);
0873                             sig_eta_evt[netamx + tmpxet][10]->Fill(tkpt03, nomHOSig);
0874                           }
0875                           variab[iring2 + ringmx][10]->Fill(tkpt03);
0876                           mu_projection[10]->Fill(ieta, iphi);
0877                           if (ips10 > 0) {
0878                             ncount[iring2][3]++;  //Efficiency of Muon detection
0879                             if (isZSps) {
0880                               ncount[iring2][20]++;
0881                               sigvsevt[iring2 + ringmx][11]->Fill(hodx, nomHOSig);
0882                               sig_eta_evt[netamx + tmpxet][11]->Fill(hodx, nomHOSig);
0883                             }
0884                             variab[iring2 + ringmx][11]->Fill(hodx);
0885                             mu_projection[11]->Fill(ieta, iphi);
0886 
0887                             if (ips11 > 0) {
0888                               if (isZSps) {
0889                                 ncount[iring2][21]++;
0890                                 sigvsevt[iring2 + ringmx][12]->Fill(hody, nomHOSig);
0891                                 sig_eta_evt[netamx + tmpxet][12]->Fill(hody, nomHOSig);
0892                               }
0893                               variab[iring2 + ringmx][12]->Fill(hody);
0894                               mu_projection[12]->Fill(ieta, iphi);
0895 
0896                               if (ips12 > 0) {
0897                                 ncount[iring2][2]++;  //Efficiency of Muon detection
0898                                 if (isZSps) {
0899                                   ncount[iring2][22]++;
0900                                   sigvsevt[iring2 + ringmx][13]->Fill(htime, nomHOSig);
0901                                   sig_eta_evt[tmpxet + ringmx][13]->Fill(htime, nomHOSig);
0902                                 }
0903                                 variab[iring2 + ringmx][13]->Fill(htime);
0904                                 mu_projection[13]->Fill(ieta, iphi);
0905 
0906                                 if (ips13 > 0) {
0907                                   if (isZSps) {
0908                                     ncount[iring2][23]++;
0909                                     mu_projection[14]->Fill(ieta, iphi);
0910                                   }
0911                                 }
0912                               }
0913                             }
0914                           }
0915                         }
0916                       }
0917                     }
0918                   }
0919                 }
0920               }
0921             }
0922           }
0923         }
0924       }
0925       if (isZSps) {
0926         sigvsevt[iring2 + 2 * ringmx][0]->Fill(abs(ndof), nomHOSig);
0927         sigvsevt[iring2 + 2 * ringmx][1]->Fill(chisq, nomHOSig);
0928         sigvsevt[iring2 + 2 * ringmx][2]->Fill(trkth, nomHOSig);
0929         sigvsevt[iring2 + 2 * ringmx][3]->Fill(trkph, nomHOSig);
0930         sigvsevt[iring2 + 2 * ringmx][4]->Fill(1000 * therr, nomHOSig);
0931         sigvsevt[iring2 + 2 * ringmx][5]->Fill(1000 * pherr, nomHOSig);
0932         if (abshoang > 0.01) {
0933           sigvsevt[iring2 + 2 * ringmx][6]->Fill(tmphoang, (nomHOSig)*abshoang);
0934         }
0935         sigvsevt[iring2 + 2 * ringmx][7]->Fill(fabs(trkmm), nomHOSig);
0936         sigvsevt[iring2 + 2 * ringmx][8]->Fill(nmuon, nomHOSig);
0937         if (!m_cosmic)
0938           sigvsevt[iring2 + 2 * ringmx][9]->Fill(caloen[0], nomHOSig);
0939         sigvsevt[iring2 + 2 * ringmx][10]->Fill(tkpt03, nomHOSig);
0940         sigvsevt[iring2 + 2 * ringmx][11]->Fill(hodx, nomHOSig);
0941         sigvsevt[iring2 + 2 * ringmx][12]->Fill(hody, nomHOSig);
0942         sigvsevt[iring2 + 2 * ringmx][13]->Fill(htime, nomHOSig);
0943 
0944         sig_eta_evt[2 * netamx + tmpxet][0]->Fill(abs(ndof), nomHOSig);
0945         sig_eta_evt[2 * netamx + tmpxet][1]->Fill(chisq, nomHOSig);
0946         sig_eta_evt[2 * netamx + tmpxet][2]->Fill(trkth, nomHOSig);
0947         sig_eta_evt[2 * netamx + tmpxet][3]->Fill(trkph, nomHOSig);
0948         sig_eta_evt[2 * netamx + tmpxet][4]->Fill(1000 * therr, nomHOSig);
0949         sig_eta_evt[2 * netamx + tmpxet][5]->Fill(1000 * pherr, nomHOSig);
0950         if (abshoang > 0.01) {
0951           sig_eta_evt[2 * netamx + tmpxet][6]->Fill(tmphoang, (nomHOSig)*abshoang);
0952         }
0953         sig_eta_evt[2 * netamx + tmpxet][7]->Fill(fabs(trkmm), nomHOSig);
0954         sig_eta_evt[2 * netamx + tmpxet][8]->Fill(nmuon, nomHOSig);
0955         if (!m_cosmic)
0956           sig_eta_evt[2 * netamx + tmpxet][9]->Fill(caloen[0], nomHOSig);
0957         sig_eta_evt[2 * netamx + tmpxet][10]->Fill(tkpt03, nomHOSig);
0958         sig_eta_evt[2 * netamx + tmpxet][11]->Fill(hodx, nomHOSig);
0959         sig_eta_evt[2 * netamx + tmpxet][12]->Fill(hody, nomHOSig);
0960         sig_eta_evt[2 * netamx + tmpxet][13]->Fill(htime, nomHOSig);
0961       }
0962 
0963       variab[iring2 + 2 * ringmx][0]->Fill(abs(ndof));
0964       variab[iring2 + 2 * ringmx][1]->Fill(chisq);
0965       variab[iring2 + 2 * ringmx][2]->Fill(trkth);
0966       variab[iring2 + 2 * ringmx][3]->Fill(trkph);
0967       variab[iring2 + 2 * ringmx][4]->Fill(1000 * therr);
0968       variab[iring2 + 2 * ringmx][5]->Fill(1000 * pherr);
0969       variab[iring2 + 2 * ringmx][6]->Fill(tmphoang);
0970       variab[iring2 + 2 * ringmx][7]->Fill(fabs(trkmm));
0971       variab[iring2 + 2 * ringmx][8]->Fill(nmuon);
0972       if (!m_cosmic)
0973         variab[iring2 + 2 * ringmx][9]->Fill(caloen[0]);
0974       variab[iring2 + 2 * ringmx][10]->Fill(tkpt03);
0975       variab[iring2 + 2 * ringmx][11]->Fill(hodx);
0976       variab[iring2 + 2 * ringmx][12]->Fill(hody);
0977       variab[iring2 + 2 * ringmx][13]->Fill(htime);
0978 
0979       muonnm->Fill(nmuon);
0980       muonmm->Fill(trkmm);
0981       muonth->Fill(trkth * 180 / pival);
0982       muonph->Fill(trkph * 180 / pival);
0983       muonch->Fill(chisq);
0984 
0985       int iselect = (ipsall == mypow_2[ncut] - 1) ? 1 : 0;
0986 
0987       if (iselect == 1) {
0988         ipass++;
0989         sel_muonnm->Fill(nmuon);
0990         sel_muonmm->Fill(trkmm);
0991         sel_muonth->Fill(trkth * 180 / pival);
0992         sel_muonph->Fill(trkph * 180 / pival);
0993         sel_muonch->Fill(chisq);
0994         if (m_histFill && tmpxet >= 0 && tmpxet < netamx && iphi >= 0 && iphi < nphimx) {
0995           ho_indenergy[tmpxet][iphi - 1]->Fill(nomHOSig);
0996         }
0997         if (m_treeFill) {
0998           T1->Fill();
0999         }
1000       }
1001     }  //for (HOCalibVariableCollection::const_iterator hoC=(*HOCalib).begin(); hoC!=(*HOCalib).end(); hoC++){
1002   }    //if (isCosMu)
1003 }
1004 
1005 //define this as a plug-in
1006 DEFINE_FWK_MODULE(HOCalibAnalyzer);