CSCObjects/test/CSCGainsStudy.cc

0001 /* EDAnalyzer to study property of gains.
0002  *
0003  * \author Dominique Fortin
0004  */
0005
0006 #include <memory>
0007 #include <iostream>
0008 #include <stdexcept>
0009
0010 // user include files
0011 #include <CondFormats/CSCObjects/test/CSCGainsStudy.h>
0012
0013 #include "FWCore/Framework/interface/ESHandle.h"
0014 #include "FWCore/Framework/interface/MakerMacros.h"
0015
0016 #include "FWCore/Framework/interface/EventSetup.h"
0017
0018 #include "TFile.h"
0019
0020 using namespace std;
0021 using namespace edm;
0022
0023 // Constructor
0024 CSCGainsStudy::CSCGainsStudy(const ParameterSet& pset) : gainsToken{esConsumes()} {
0025   // Get the various input parameters
0026   debug = pset.getUntrackedParameter<bool>("debug");
0027   rootFileName = pset.getUntrackedParameter<string>("rootFileName");
0028
0029   if (debug)
0030     cout << "[CSCGainsStudy] Constructor called" << endl;
0031
0032   // Create the root file
0033   theFile = new TFile(rootFileName.c_str(), "RECREATE");
0034   theFile->cd();
0035
0036   // Book the histograms
0037   All_CSC = new HCSCGains("All_CSC");
0038   // ME+1/1
0039   ME_11_27 = new HCSCGains("ME_11_27");
0040   ME_11_28 = new HCSCGains("ME_11_28");
0041   ME_11_29 = new HCSCGains("ME_11_29");
0042   ME_11_30 = new HCSCGains("ME_11_30");
0043   ME_11_31 = new HCSCGains("ME_11_31");
0044   ME_11_32 = new HCSCGains("ME_11_32");
0045   // ME+1/2
0046   ME_12_27 = new HCSCGains("ME_12_27");
0047   ME_12_28 = new HCSCGains("ME_12_28");
0048   ME_12_29 = new HCSCGains("ME_12_29");
0049   ME_12_30 = new HCSCGains("ME_12_30");
0050   ME_12_31 = new HCSCGains("ME_12_31");
0051   ME_12_32 = new HCSCGains("ME_12_32");
0052   // ME+1/3
0053   ME_13_27 = new HCSCGains("ME_13_27");
0054   ME_13_28 = new HCSCGains("ME_13_28");
0055   ME_13_29 = new HCSCGains("ME_13_29");
0056   ME_13_30 = new HCSCGains("ME_13_30");
0057   ME_13_31 = new HCSCGains("ME_13_31");
0058   ME_13_32 = new HCSCGains("ME_13_32");
0059   // ME+2/1
0060   ME_21_14 = new HCSCGains("ME_21_14");
0061   ME_21_15 = new HCSCGains("ME_21_15");
0062   ME_21_16 = new HCSCGains("ME_21_16");
0063   // ME+2/2
0064   ME_22_27 = new HCSCGains("ME_22_27");
0065   ME_22_28 = new HCSCGains("ME_22_28");
0066   ME_22_29 = new HCSCGains("ME_22_29");
0067   ME_22_30 = new HCSCGains("ME_22_30");
0068   ME_22_31 = new HCSCGains("ME_22_31");
0069   ME_22_32 = new HCSCGains("ME_22_32");
0070   // ME+3/1
0071   ME_31_14 = new HCSCGains("ME_31_14");
0072   ME_31_15 = new HCSCGains("ME_31_15");
0073   ME_31_16 = new HCSCGains("ME_31_16");
0074   // ME+3/2
0075   ME_32_27 = new HCSCGains("ME_32_27");
0076   ME_32_28 = new HCSCGains("ME_32_28");
0077   ME_32_29 = new HCSCGains("ME_32_29");
0078   ME_32_30 = new HCSCGains("ME_32_30");
0079   ME_32_31 = new HCSCGains("ME_32_31");
0080   ME_32_32 = new HCSCGains("ME_32_32");
0081 }
0082
0083 // Destructor
0084 CSCGainsStudy::~CSCGainsStudy() {
0085   if (debug)
0086     cout << "[CSCGainsStudy] Destructor called" << endl;
0087
0088   // Write the histos to file
0089   theFile->cd();
0090   All_CSC->Write();
0091   // ME+1/1
0092   ME_11_27->Write();
0093   ME_11_28->Write();
0094   ME_11_29->Write();
0095   ME_11_30->Write();
0096   ME_11_31->Write();
0097   ME_11_32->Write();
0098   // ME+1/2
0099   ME_12_27->Write();
0100   ME_12_28->Write();
0101   ME_12_29->Write();
0102   ME_12_30->Write();
0103   ME_12_31->Write();
0104   ME_12_32->Write();
0105   // ME+1/3
0106   ME_13_27->Write();
0107   ME_13_28->Write();
0108   ME_13_29->Write();
0109   ME_13_30->Write();
0110   ME_13_31->Write();
0111   ME_13_32->Write();
0112   // ME+2/1
0113   ME_21_14->Write();
0114   ME_21_15->Write();
0115   ME_21_16->Write();
0116   // ME+2/2
0117   ME_22_27->Write();
0118   ME_22_28->Write();
0119   ME_22_29->Write();
0120   ME_22_30->Write();
0121   ME_22_31->Write();
0122   ME_22_32->Write();
0123   // ME+3/1
0124   ME_31_14->Write();
0125   ME_31_15->Write();
0126   ME_31_16->Write();
0127   // ME+3/2
0128   ME_32_27->Write();
0129   ME_32_28->Write();
0130   ME_32_29->Write();
0131   ME_32_30->Write();
0132   ME_32_31->Write();
0133   ME_32_32->Write();
0134
0135   // Close file
0136   theFile->Close();
0137   if (debug)
0138     cout << "************* Finished writing histograms to file" << endl;
0139
0140   delete theFile;
0141 }
0142
0143 /* analyze
0144  *
0145  */
0146 void CSCGainsStudy::analyze(const Event& event, const EventSetup& eventSetup) {
0147   // Get the gains and compute global gain average to store for later use in strip calibration
0148   const CSCGains* hGains_ = &eventSetup.getData(gainsToken);
0149
0150   // Store so it can be used in member functions
0151   pGains = hGains_;
0152
0153   // Get global gain average:
0154   float AvgStripGain = getStripGainAvg();
0155
0156   if (debug)
0157     std::cout << "Global average gain is " << AvgStripGain << std::endl;
0158
0159   HCSCGains* histo = 0;
0160
0161   TString prefix = "ME_";
0162   float thegain = -10.;
0163   float weight = -10.;
0164   float weight0 = 1.;
0165
0166   // Build iterator which loops on all layer id:
0167   map<int, vector<CSCGains::Item> >::const_iterator it;
0168
0169   for (it = pGains->gains.begin(); it != pGains->gains.end(); ++it) {
0170     // Channel id used for retrieving gains from database is
0171     // chId=220000000 + ec*100000 + st*10000 + rg*1000 + ch*10 + la;
0172
0173     const unsigned ChId = it->first;
0174     unsigned offset = (ChId - 220000000);
0175     unsigned ec = offset / 100000;                                            // endcap
0176     unsigned st = (offset - ec * 100000) / 10000;                             // station
0177     unsigned rg = (offset - ec * 100000 - st * 10000) / 1000;                 // ring
0178     unsigned ch = (offset - ec * 100000 - st * 10000 - rg * 1000) / 10;       // chamber
0179     unsigned la = (offset - ec * 100000 - st * 10000 - rg * 1000 - ch * 10);  // layer
0180
0181     if (la == 0)
0182       continue;  // layer == 0 means whole chamber...
0183
0184     int channel = 1;
0185
0186     // Build iterator which loops on all channels:
0187     vector<CSCGains::Item>::const_iterator gain_i;
0188
0189     for (gain_i = it->second.begin(); gain_i != it->second.end(); ++gain_i) {
0190       thegain = gain_i->gain_slope;
0191       weight = AvgStripGain / thegain;
0192
0193       //      if (debug) std::cout << "the weight is " << weight << std::endl;
0194
0195       if (weight <= 0.)
0196         weight = -10.;  // ignore strips with no gain computed
0197       if (channel % 100 == 1)
0198         weight0 = weight;  // can't make comparison with first strip, so set diff to zero
0199
0200       // Fill corrections factor for all strip at once
0201       histo = All_CSC;
0202       histo->Fill(weight, weight0, channel, la);
0203
0204       // Now fill correction factor for given chamber
0205       // Get station
0206       if (st == 1) {  // station 1
0207
0208         if (rg == 1) {  // ring 1
0209           if (ch == 27) {
0210             histo = ME_11_27;
0211           } else if (ch == 28) {
0212             histo = ME_11_28;
0213           } else if (ch == 29) {
0214             histo = ME_11_29;
0215           } else if (ch == 30) {
0216             histo = ME_11_30;
0217           } else if (ch == 31) {
0218             histo = ME_11_31;
0219           } else {
0220             histo = ME_11_32;
0221           }
0222         } else if (rg == 2) {  // ring 2
0223           if (ch == 27) {
0224             histo = ME_12_27;
0225           } else if (ch == 28) {
0226             histo = ME_12_28;
0227           } else if (ch == 29) {
0228             histo = ME_12_29;
0229           } else if (ch == 30) {
0230             histo = ME_12_30;
0231           } else if (ch == 31) {
0232             histo = ME_12_31;
0233           } else {
0234             histo = ME_12_32;
0235           }
0236         } else {  // ring 3
0237           if (ch == 27) {
0238             histo = ME_13_27;
0239           } else if (ch == 28) {
0240             histo = ME_13_28;
0241           } else if (ch == 29) {
0242             histo = ME_13_29;
0243           } else if (ch == 30) {
0244             histo = ME_13_30;
0245           } else if (ch == 31) {
0246             histo = ME_13_31;
0247           } else {
0248             histo = ME_13_32;
0249           }
0250         }
0251
0252       } else if (st == 2) {  // station 2
0253
0254         if (rg == 1) {  // ring 1
0255           if (ch == 14) {
0256             histo = ME_21_14;
0257           } else if (ch == 15) {
0258             histo = ME_21_15;
0259           } else {
0260             histo = ME_21_16;
0261           }
0262         } else {  // ring 2
0263           if (ch == 27) {
0264             histo = ME_22_27;
0265           } else if (ch == 28) {
0266             histo = ME_22_28;
0267           } else if (ch == 29) {
0268             histo = ME_22_29;
0269           } else if (ch == 30) {
0270             histo = ME_22_30;
0271           } else if (ch == 31) {
0272             histo = ME_22_31;
0273           } else {
0274             histo = ME_22_32;
0275           }
0276         }
0277
0278       } else {  // station 3
0279
0280         if (rg == 1) {  // ring 1
0281           if (ch == 14) {
0282             histo = ME_31_14;
0283           } else if (ch == 15) {
0284             histo = ME_31_15;
0285           } else {
0286             histo = ME_31_16;
0287           }
0288         } else {  // ring 2
0289           if (ch == 27) {
0290             histo = ME_32_27;
0291           } else if (ch == 28) {
0292             histo = ME_32_28;
0293           } else if (ch == 29) {
0294             histo = ME_32_29;
0295           } else if (ch == 30) {
0296             histo = ME_32_30;
0297           } else if (ch == 31) {
0298             histo = ME_32_31;
0299           } else {
0300             histo = ME_32_32;
0301           }
0302         }
0303       }
0304       histo->Fill(weight, weight0, channel, la);
0305       weight0 = weight;
0306       channel++;
0307     }
0308   }
0309 }
0310
0311 /* getStripGainAvg
0312  *
0313  */
0314 float CSCGainsStudy::getStripGainAvg() {
0315   int n_strip = 0;
0316   float gain_tot = 0.;
0317   float gain_avg = -1.;
0318
0319   // Build iterator which loops on all layer id:
0320   map<int, vector<CSCGains::Item> >::const_iterator it;
0321
0322   for (it = pGains->gains.begin(); it != pGains->gains.end(); ++it) {
0323     const unsigned ChId = it->first;
0324     unsigned offset = (ChId - 220000000);
0325     unsigned ec = offset / 100000;                                            // endcap
0326     unsigned st = (offset - ec * 100000) / 10000;                             // station
0327     unsigned rg = (offset - ec * 100000 - st * 10000) / 1000;                 // ring
0328     unsigned ch = (offset - ec * 100000 - st * 10000 - rg * 1000) / 10;       // chamber
0329     unsigned la = (offset - ec * 100000 - st * 10000 - rg * 1000 - ch * 10);  // layer
0330
0331     if (la == 0)
0332       continue;  // layer == 0 means whole chamber...
0333
0334     // Build iterator which loops on all channels:
0335     vector<CSCGains::Item>::const_iterator gain_i;
0336
0337     for (gain_i = it->second.begin(); gain_i != it->second.end(); ++gain_i) {
0338       // Make sure channel isn't dead, otherwise don't include in average !
0339       if (gain_i->gain_slope > 0.) {
0340         gain_tot += gain_i->gain_slope;
0341         n_strip++;
0342       }
0343     }
0344   }
0345   // Average gain
0346   if (n_strip > 0)
0347     gain_avg = gain_tot / n_strip;
0348
0349   return gain_avg;
0350 }
0351
0352 DEFINE_FWK_MODULE(CSCGainsStudy);