Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:35:30 

0001 #ifndef HcalObjRepresent_h
0002 #define HcalObjRepresent_h
0003 
0004 #include "CondFormats/HcalObjects/interface/HcalGains.h"
0005 
0006 //#include "CondCore/Utilities/interface/PayLoadInspector.h"
0007 //#include "CondCore/Utilities/interface/InspectorPythonWrapper.h"
0008 #include "CondFormats/HcalObjects/interface/HcalDetIdRelationship.h"
0009 #include "CondFormats/HcalObjects/interface/HcalCondObjectContainer.h"
0010 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0011 #include "DataFormats/HcalDetId/interface/HcalOtherDetId.h"
0012 #include "DataFormats/HcalDetId/interface/HcalCastorDetId.h"
0013 #include "DataFormats/HcalDetId/interface/HcalZDCDetId.h"
0014 #include "FWCore/Utilities/interface/Exception.h"
0015 
0016 #include <string>
0017 #include <fstream>
0018 #include <sstream>
0019 #include <vector>
0020 
0021 #include "TH1F.h"
0022 #include "TH2F.h"
0023 #include "DataFormats/DetId/interface/DetId.h"
0024 
0025 #include "TROOT.h"
0026 #include "TCanvas.h"
0027 #include "TStyle.h"
0028 #include "TColor.h"
0029 #include "TLine.h"
0030 #include "TLatex.h"
0031 #include "TProfile.h"
0032 #include "TPaveLabel.h"
0033 
0034 //functions for correct representation of data in summary and plot
0035 namespace HcalObjRepresent {
0036 
0037   //used to produce all display objects for payload inspector
0038   template <class Items, class Item>
0039   class HcalDataContainer {
0040   public:
0041     HcalDataContainer(std::shared_ptr<Items> payload, unsigned int run) : payload_(payload), run_(run) {
0042       PlotMode_ = "Map";
0043     }
0044 
0045     virtual ~HcalDataContainer() {}
0046     // For easier channel mapping
0047     typedef std::tuple<int, int, int> Coord;
0048     typedef std::map<Coord, Item> tHcalValCont;
0049     // mapping of pair of subdetector name (e.g, "HE") and depth number (e.g. 3) to Histogram of data for that subdetector/depth pair
0050     typedef std::map<std::pair<std::string, int>, TH2F*> DepthMap;
0051 
0052     ///////////////// public Get functions  /////////////////
0053     unsigned int GetRun() { return run_; }
0054 
0055     std::string GetTopoMode() { return TopoMode_; }
0056 
0057     std::map<std::string, int> GetSubDetDepths() { return subDetDepths_; }
0058 
0059     DepthMap GetDepths() {
0060       fillValConts();
0061       //std::cout << "Got depths with run number = " << std::to_string(GetRun()) << std::endl;
0062       return depths_;
0063     }
0064 
0065     ///////////////// setting object fields /////////////////
0066 
0067     // Fills a tHcalValCont for each subdetector, setting Topology Mode along the way
0068     void fillValConts() {
0069       if (!depths_.empty())
0070         return;
0071       int iphi, ieta, depth;
0072       HcalDetId detId;
0073 
0074       std::string subDetName;
0075       std::vector<Item> itemsVec;
0076       std::pair<std::string, int> depthKey;
0077       const char* histLabel;
0078       for (std::pair<std::string, std::vector<Item> > cont : (*payload_).getAllContainers()) {
0079         subDetName = std::get<0>(cont);
0080         itemsVec = std::get<1>(cont);
0081 
0082         auto valContainer = getContFromString(subDetName);
0083 
0084         for (Item item : itemsVec) {
0085           detId = HcalDetId(item.rawId());
0086           iphi = detId.iphi();
0087           ieta = detId.ieta();
0088           depth = detId.depth();
0089           Coord coord = std::make_tuple(depth, ieta, iphi);
0090           //Add hist if it's not there, AND if it is one we care about: HO,HB,HE,HF, AND it's not an empty entry (depth = 0 in HO)
0091           if (subDetName[0] == 'H' && depth != 0) {
0092             valContainer->insert(std::pair<std::tuple<int, int, int>, Item>(coord, item));
0093 
0094             depthKey = std::make_pair(subDetName, depth);
0095 
0096             auto depthVal = depths_.find(depthKey);
0097             if (depthVal == depths_.end()) {
0098               histLabel = ("run" + std::to_string(run_) + "_" + subDetName + "_d" + std::to_string(depth)).c_str();
0099               depths_.insert(std::make_pair(std::make_pair(subDetName, depth),
0100                                             new TH2F(histLabel, histLabel, 83, -42.5, 41.5, 71, 0.5, 71.5)));
0101             }
0102             depths_[depthKey]->Fill(ieta, iphi, getValue(&item));
0103           }
0104         }
0105       }
0106 
0107       //Still need to know which hists to take when done; decide now
0108       setTopoModeFromValConts();
0109     }
0110 
0111     const std::vector<std::pair<std::string, std::vector<Item> > > getAllItems() {
0112       return (*payload_).getAllContainers();
0113     }
0114 
0115     ////NOTE to be implemented in PayloadInspector classes
0116     virtual float getValue(const Item* item) {
0117       throw cms::Exception("Value definition not found") << "getValue definition not found for " << payload_->myname();
0118     }
0119 
0120     //Gets Hcal Object at given coordinate
0121     //Currently unused but remains as a potentially useful function
0122     Item* getItemFromValCont(std::string subDetName, int depth, int ieta, int iphi, bool throwOnFail) {
0123       Item* cell = nullptr;
0124       Coord coord = std::make_tuple(depth, ieta, iphi);
0125       tHcalValCont* valContainer = getContFromString(subDetName);
0126 
0127       auto it = valContainer->find(coord);
0128       if (it != valContainer->end())
0129         cell = &it->second;
0130 
0131       if ((!cell)) {
0132         if (throwOnFail) {
0133           throw cms::Exception("Conditions not found")
0134               << "Unavailable Conditions of type " << payload_->myname() << " for cell " << subDetName << " (" << depth
0135               << ", " << ieta << ", " << iphi << ")";
0136         }
0137       }
0138       return cell;
0139     }
0140 
0141     /////////////////// Building Graphics //////////////////////
0142 
0143     //TODO: remove zero entries from doing divide and subtract
0144 
0145     // To generate Ratios of two IOVs
0146     void Divide(HcalDataContainer* dataCont2) {
0147       //TODO: Do something like looping over this and that depths setting every empty bin (which I think means content=0) to -999. Then replacing the divide call with another manual loop that sets all bins with content -999 to 0
0148 
0149       PlotMode_ = "Ratio";
0150       DepthMap::iterator depth1;
0151       std::pair<std::string, int> key;
0152       DepthMap myDepths = this->GetDepths();
0153       DepthMap theirDepths = dataCont2->GetDepths();
0154       for (depth1 = myDepths.begin(); depth1 != myDepths.end(); depth1++) {
0155         key = depth1->first;
0156         if (theirDepths.count(key) != 0) {
0157           myDepths.at(key)->Divide((const TH1*)theirDepths.at(key));
0158         } else {
0159           throw cms::Exception("Unaligned Conditions")
0160               << "trying to plot conditions for " << payload_->myname() << "; found value for " << std::get<0>(key)
0161               << " depth " << std::to_string(std::get<1>(key)) << " in run " << GetRun() << " but not in run "
0162               << dataCont2->GetRun();
0163         }
0164       }
0165     }
0166 
0167     // To generate Diffs of two IOVs
0168     void Subtract(HcalDataContainer* dataCont2) {
0169       PlotMode_ = "Diff";
0170       DepthMap::iterator depth1;
0171       std::pair<std::string, int> key;
0172       DepthMap myDepths = this->GetDepths();
0173       DepthMap theirDepths = dataCont2->GetDepths();
0174       for (depth1 = myDepths.begin(); depth1 != myDepths.end(); depth1++) {
0175         key = depth1->first;
0176         if (theirDepths.count(key) != 0) {
0177           myDepths.at(key)->Add(myDepths.at(key), theirDepths.at(key), 1, -1);
0178         } else {
0179           throw cms::Exception("Unaligned Conditions")
0180               << "trying to plot conditions for " << payload_->myname() << "; found value for " << std::get<0>(key)
0181               << " depth " << std::to_string(std::get<1>(key)) << " in run " << GetRun() << " but not in run "
0182               << dataCont2->GetRun();
0183         }
0184       }
0185     }
0186 
0187     // wisely determines what range to set histogram axis to
0188     std::pair<float, float> GetRange(TH1* hist) {
0189       if (PlotMode_ == "Ratio") {
0190         float amp;
0191         Double_t adjustMin = 1;
0192         Double_t tempMin;
0193         int nBinsX = hist->GetXaxis()->GetNbins();
0194         int nBinsY = hist->GetYaxis()->GetNbins();
0195         for (int i = 0; i < nBinsX; i++) {
0196           for (int j = 0; j < nBinsY; j++) {
0197             tempMin = hist->GetBinContent(i, j);
0198             if ((tempMin != 0) && (tempMin < adjustMin))
0199               adjustMin = tempMin;
0200           }
0201         }
0202         amp = std::max((1 - adjustMin), (hist->GetMaximum() - 1));
0203         //amp = std::max((1 - hist->GetMinimum()),(hist->GetMaximum() - 1) );
0204         return std::make_pair(1 - amp, 1 + amp);
0205       } else if (PlotMode_ == "Diff") {
0206         float amp;
0207         amp = std::max((0 - hist->GetMinimum()), hist->GetMaximum());
0208         return std::make_pair((-1 * amp), amp);
0209       } else {
0210         Double_t adjustMin = 10000;
0211         Double_t tempMin;
0212         int nBinsX = hist->GetXaxis()->GetNbins();
0213         int nBinsY = hist->GetYaxis()->GetNbins();
0214         for (int i = 0; i < nBinsX; i++) {
0215           for (int j = 0; j < nBinsY; j++) {
0216             tempMin = hist->GetBinContent(i, j);
0217             if ((tempMin != 0) && (tempMin < adjustMin))
0218               adjustMin = tempMin;
0219           }
0220         }
0221         return std::make_pair(((adjustMin == 10000) ? hist->GetMinimum() : adjustMin), hist->GetMaximum());
0222       }
0223     }
0224 
0225     // set style
0226     void initGraphics() {
0227       gStyle->SetOptStat(0);
0228       gStyle->SetPalette(1);
0229       gStyle->SetOptFit(0);
0230       gStyle->SetLabelFont(42);
0231       gStyle->SetLabelFont(42);
0232       gStyle->SetTitleFont(42);
0233       gStyle->SetTitleFont(42);
0234       gStyle->SetMarkerSize(0);
0235       gStyle->SetTitleOffset(1.3, "Y");
0236       gStyle->SetTitleOffset(1.0, "X");
0237       gStyle->SetNdivisions(510);
0238       gStyle->SetStatH(0.11);
0239       gStyle->SetStatW(0.33);
0240       gStyle->SetTitleW(0.4);
0241       gStyle->SetTitleX(0.13);
0242       gStyle->SetPadTickX(1);
0243       gStyle->SetPadTickY(1);
0244     }
0245 
0246     TH1D* GetProjection(TH2F* hist, std::string plotType, const char* newName, std::string subDetName, int depth) {
0247       //TODO: Also want average for standard projection of 2DHist (not ratio or diff)?
0248       //if (PlotMode_ != "Ratio") return (plotType=="EtaProfile") ? ((TH2F*)(hist->Clone("temp")))->ProjectionX(newName) : ((TH2F*)(hist->Clone("temp")))->ProjectionY(newName);
0249 
0250       //TH1D* projection = ((TH2F*)(depths_[std::make_pair(subDetName,depth)]->Clone("temp")))->ProjectionX(newName);
0251 
0252       int xBins = (plotType == "EtaProfile") ? 83 : 71;
0253       int etaMin = -42, etaMax = 42, phiMin = 1, phiMax = 72;
0254       int xMin = (plotType == "EtaProfile") ? etaMin : phiMin;
0255       int xMax = (plotType == "EtaProfile") ? etaMax : phiMax;
0256       int otherMin = (plotType == "EtaProfile") ? phiMin : etaMin;
0257       int otherMax = (plotType == "EtaProfile") ? phiMax : etaMax;
0258       TH1D* retHist = new TH1D(newName, newName, xBins, xMin, xMax);
0259       int numChannels;
0260       Double_t sumVal;
0261       Double_t channelVal;
0262       int ieta, iphi;
0263       int bin = 0;
0264       for (int i = xMin; i <= xMax; i++ && bin++) {
0265         numChannels = 0;
0266         sumVal = 0;
0267         for (int j = otherMin; j <= otherMax; j++) {
0268           ieta = (plotType == "EtaProfile") ? i : j;
0269           ieta += 42;
0270           iphi = (plotType == "EtaProfile") ? j : i;
0271           iphi += -1;
0272           channelVal = hist->GetBinContent(ieta, iphi);
0273           //std::cout << "(ieta, iphi)    :       (" << std::to_string(ieta) << ", " << std::to_string(iphi) << ")" << std::endl;
0274           if (channelVal != 0) {
0275             sumVal += channelVal;
0276             numChannels++;
0277           }
0278         }
0279         //if(sumVal !=0) projection->SetBinContent(i, sumVal/((Double_t)numChannels));//retHist->Fill(i,sumVal/((Double_t)numChannels));
0280         if (sumVal != 0)
0281           retHist->Fill(i, sumVal / ((Double_t)numChannels));
0282       }
0283       return retHist;
0284       //return projection;
0285     }
0286 
0287     // fills a canvas with given subdetector information, plotting all depths
0288     void FillCanv(TCanvas* canvas,
0289                   std::string subDetName,
0290                   int startDepth = 1,
0291                   int startCanv = 1,
0292                   std::string plotForm = "2DHist") {
0293       const char* newName;
0294       std::pair<float, float> range;
0295       int padNum;
0296       int maxDepth = (subDetName == "HO") ? 4 : subDetDepths_[subDetName];
0297       TH1D* projection;
0298       TLatex label;
0299       for (int i = startDepth; i <= maxDepth; i++) {
0300         //skip if data not obtained; TODO: Maybe add text on plot saying data not found?
0301         if (depths_.count(std::make_pair(subDetName, i)) == 0) {
0302           return;
0303         }
0304 
0305         padNum = i + startCanv - 1;
0306         if (subDetName == "HO")
0307           padNum = padNum - 3;
0308         canvas->cd(padNum);
0309         canvas->GetPad(padNum)->SetGridx(1);
0310         canvas->GetPad(padNum)->SetGridy(1);
0311 
0312         if (plotForm == "2DHist") {
0313           canvas->GetPad(padNum)->SetRightMargin(0.13);
0314           range = GetRange(depths_[std::make_pair(subDetName, i)]);
0315           depths_[std::make_pair(subDetName, i)]->Draw("colz");
0316           depths_[std::make_pair(subDetName, i)]->SetContour(100);
0317           depths_[std::make_pair(subDetName, i)]->GetXaxis()->SetTitle("ieta");
0318           depths_[std::make_pair(subDetName, i)]->GetYaxis()->SetTitle("iphi");
0319           depths_[std::make_pair(subDetName, i)]->GetXaxis()->CenterTitle();
0320           depths_[std::make_pair(subDetName, i)]->GetYaxis()->CenterTitle();
0321           depths_[std::make_pair(subDetName, i)]->GetZaxis()->SetRangeUser(std::get<0>(range), std::get<1>(range));
0322           depths_[std::make_pair(subDetName, i)]->GetYaxis()->SetTitleSize(0.06);
0323           depths_[std::make_pair(subDetName, i)]->GetYaxis()->SetTitleOffset(0.80);
0324           depths_[std::make_pair(subDetName, i)]->GetXaxis()->SetTitleSize(0.06);
0325           depths_[std::make_pair(subDetName, i)]->GetXaxis()->SetTitleOffset(0.80);
0326           depths_[std::make_pair(subDetName, i)]->GetYaxis()->SetLabelSize(0.055);
0327           depths_[std::make_pair(subDetName, i)]->GetXaxis()->SetLabelSize(0.055);
0328         } else {
0329           canvas->GetPad(padNum)->SetLeftMargin(0.152);
0330           canvas->GetPad(padNum)->SetRightMargin(0.02);
0331           //gStyle->SetTitleOffset(1.6,"Y");
0332           newName = ("run_" + std::to_string(run_) + "_" + subDetName + "_d" + std::to_string(i) + "_" +
0333                      (plotForm == "EtaProfile" ? "ieta" : "iphi"))
0334                         .c_str();
0335           //projection = ((TH2F*)(depths_[std::make_pair(subDetName,i)]->Clone("temp")))->ProjectionX(newName);
0336           projection = GetProjection(depths_[std::make_pair(subDetName, i)], plotForm, newName, subDetName, i);
0337           range = GetRange(projection);
0338           projection->Draw("hist");
0339           projection->GetXaxis()->SetTitle((plotForm == "EtaProfile" ? "ieta" : "iphi"));
0340           projection->GetXaxis()->CenterTitle();
0341           projection->GetYaxis()->SetTitle(
0342               (payload_->myname() + " " + (PlotMode_ == "Map" ? "" : PlotMode_) + " " + GetUnit(payload_->myname()))
0343                   .c_str());
0344           label.SetNDC();
0345           label.SetTextAlign(26);
0346           label.SetTextSize(0.05);
0347           label.DrawLatex(
0348               0.3, 0.95, ("run_" + std::to_string(run_) + "_" + subDetName + "_d" + std::to_string(i)).c_str());
0349           projection->GetYaxis()->CenterTitle();
0350           projection->GetXaxis()->SetTitleSize(0.06);
0351           projection->GetYaxis()->SetTitleSize(0.06);
0352           projection->GetXaxis()->SetTitleOffset(0.80);
0353           projection->GetXaxis()->SetLabelSize(0.055);
0354           projection->GetYaxis()->SetTitleOffset(1.34);
0355           projection->GetYaxis()->SetLabelSize(0.055);
0356         }
0357       }
0358     }
0359 
0360     //// functions called in Payload Inspector classes to import final canvases to be plotted
0361 
0362     // profile = "EtaProfile" || "PhiProfile"
0363     TCanvas* getCanvasAll(std::string profile = "2DHist") {
0364       fillValConts();
0365       initGraphics();
0366       TCanvas* HAll = new TCanvas("HAll", "HAll", 1680, (GetTopoMode() == "2015/2016") ? 1680 : 2500);
0367       HAll->Divide(3, (GetTopoMode() == "2015/2016") ? 3 : 6, 0.02, 0.01);
0368       FillCanv(HAll, "HB", 1, 1, profile);
0369       FillCanv(HAll, "HO", 4, 3, profile);
0370       FillCanv(HAll, "HF", 1, 4, profile);
0371       FillCanv(HAll, "HE", 1, (GetTopoMode() == "2015/2016") ? 7 : 10, profile);
0372       return HAll;
0373     }
0374 
0375     TCanvas* getCanvasHF() {
0376       fillValConts();
0377       initGraphics();
0378       TCanvas* HF = new TCanvas("HF", "HF", 1600, 1000);
0379       HF->Divide(3, 2, 0.02, 0.01);
0380       FillCanv(HF, "HF");
0381       return HF;
0382     }
0383     TCanvas* getCanvasHE() {
0384       fillValConts();
0385       initGraphics();
0386       TCanvas* HE = new TCanvas("HE", "HE", 1680, 1680);
0387       HE->Divide(3, 3, 0.02, 0.01);
0388       FillCanv(HE, "HE");
0389       return HE;
0390     }
0391     TCanvas* getCanvasHBHO() {
0392       fillValConts();
0393       initGraphics();
0394       TCanvas* HBHO = new TCanvas("HBHO", "HBHO", 1680, 1680);
0395       HBHO->Divide(3, 3, 0.02, 0.01);
0396       FillCanv(HBHO, "HB");
0397       FillCanv(HBHO, "HO", 4, 3);
0398       FillCanv(HBHO, "HB", 1, 4, "EtaProfile");
0399       FillCanv(HBHO, "HO", 4, 6, "EtaProfile");
0400       FillCanv(HBHO, "HB", 1, 7, "PhiProfile");
0401       FillCanv(HBHO, "HO", 4, 9, "PhiProfile");
0402       return HBHO;
0403     }
0404 
0405     std::string GetUnit(std::string type) {
0406       std::string unit = units_[type];
0407       if (unit.empty())
0408         return "";
0409       else
0410         return "(" + unit + ")";
0411     }
0412 
0413   private:
0414     DepthMap depths_;
0415     std::shared_ptr<Items> payload_;
0416     unsigned int run_;
0417     std::string TopoMode_;
0418     // "Map", "Ratio", or "Diff"
0419     std::string PlotMode_;
0420 
0421     tHcalValCont HBvalContainer;
0422     tHcalValCont HEvalContainer;
0423     tHcalValCont HOvalContainer;
0424     tHcalValCont HFvalContainer;
0425     tHcalValCont HTvalContainer;
0426     tHcalValCont ZDCvalContainer;
0427     tHcalValCont CALIBvalContainer;
0428     tHcalValCont CASTORvalContainer;
0429     std::map<std::string, int> subDetDepths_;
0430     std::map<std::string, std::string> units_ = {{"HcalPedestals", "ADC"},
0431                                                  {"HcalGains", ""},             //dimensionless TODO: verify
0432                                                  {"HcalL1TriggerObjects", ""},  //dimensionless TODO: Verify
0433                                                  {"HcalPedestalWidths", "ADC"},
0434                                                  {"HcalRespCorrs", ""},  //dimensionless TODO: verify
0435                                                  {"Dark Current", ""},
0436                                                  {"fcByPE", ""},
0437                                                  {"crossTalk", ""},
0438                                                  {"parLin", ""}};
0439 
0440     tHcalValCont* getContFromString(std::string subDetString) {
0441       if (subDetString == "HB")
0442         return &HBvalContainer;
0443       else if (subDetString == "HE")
0444         return &HEvalContainer;
0445       else if (subDetString == "HF")
0446         return &HFvalContainer;
0447       else if (subDetString == "HO")
0448         return &HOvalContainer;
0449       else if (subDetString == "HT")
0450         return &HTvalContainer;
0451       else if (subDetString == "CALIB")
0452         return &CALIBvalContainer;
0453       else if (subDetString == "CASTOR")
0454         return &CASTORvalContainer;
0455       //else return &ZDCvalContainer;
0456       else if (subDetString == "ZDC_EM" || subDetString == "ZDC" || subDetString == "ZDC_HAD" ||
0457                subDetString == "ZDC_LUM")
0458         return &ZDCvalContainer;
0459       else
0460         throw cms::Exception("subDetString " + subDetString + " not found in Item");
0461     }
0462 
0463     void setTopoModeFromValConts(bool throwOnFail = false) {
0464       // Check HEP17 alternate channel for 2017, just by checking if the 7th depth is there, or if HF has 4 depths
0465       if (depths_.count(std::make_pair("HF", 4)) != 0 || depths_.count(std::make_pair("HE", 7)) != 0)
0466         TopoMode_ = "2017";
0467       // Check endcap depth unique to 2018
0468       else if (HEvalContainer.count(std::make_tuple(7, -26, 63)) != 0)
0469         TopoMode_ = "2018";
0470       // if not 2017 or 2018, 2015 and 2016 are the same
0471       else
0472         TopoMode_ = "2015/2016";
0473       //NOTE: HO's one depth is labeled depth 4
0474       if (TopoMode_ == "2018" || TopoMode_ == "2017") {
0475         subDetDepths_.insert(std::pair<std::string, int>("HB", 2));
0476         subDetDepths_.insert(std::pair<std::string, int>("HE", 7));
0477         subDetDepths_.insert(std::pair<std::string, int>("HF", 4));
0478         subDetDepths_.insert(std::pair<std::string, int>("HO", 1));
0479       } else if (TopoMode_ == "2015/2016") {
0480         subDetDepths_.insert(std::pair<std::string, int>("HB", 2));
0481         subDetDepths_.insert(std::pair<std::string, int>("HE", 3));
0482         subDetDepths_.insert(std::pair<std::string, int>("HF", 2));
0483         subDetDepths_.insert(std::pair<std::string, int>("HO", 1));
0484       }
0485     }
0486   };
0487 
0488   inline void drawTable(int nbRows, int nbColumns) {
0489     TLine* l = new TLine;
0490     l->SetLineWidth(1);
0491     for (int i = 1; i < nbRows; i++) {
0492       double y = (double)i;
0493       l = new TLine(0., y, nbColumns, y);
0494       l->Draw();
0495     }
0496 
0497     for (int i = 1; i < nbColumns; i++) {
0498       double x = (double)i;
0499       double y = (double)nbRows;
0500       l = new TLine(x, 0., x, y);
0501       l->Draw();
0502     }
0503   }
0504 
0505   inline std::string SciNotatStr(float num) {
0506     // Create an output string stream
0507     std::ostringstream streamObj2;
0508 
0509     if (num == -1)
0510       return "NOT FOUND";
0511 
0512     //Add double to stream
0513     streamObj2 << num;
0514     // Get string from output string stream
0515     std::string strObj2 = streamObj2.str();
0516 
0517     return strObj2;
0518   }
0519 
0520   inline std::string IntToBinary(unsigned int number) {
0521     std::stringstream ss;
0522     unsigned int mask = 1 << 31;
0523     for (unsigned short int i = 0; i < 32; ++i) {
0524       //if (!(i % 4))
0525       //    ss << "_";
0526       if (mask & number)
0527         ss << "1";
0528       else
0529         ss << "0";
0530       mask = mask >> 1;
0531     }
0532     return ss.str();
0533   }
0534 
0535   inline const bool isBitSet(unsigned int bitnumber, unsigned int status) {
0536     unsigned int statadd = 0x1 << (bitnumber);
0537     return (status & statadd) ? (true) : (false);
0538   }
0539 
0540   inline std::string getBitsSummary(uint32_t bits, std::string statusBitArray[], short unsigned int bitMap[]) {
0541     std::stringstream ss;
0542     for (unsigned int i = 0; i < 9; ++i) {
0543       if (isBitSet(bitMap[i], bits)) {
0544         ss << "[" << bitMap[i] << "]" << statusBitArray[bitMap[i]] << "; ";
0545       }
0546     }
0547     ss << std::endl;
0548     return ss.str();
0549   }
0550 
0551   //functions for making plot:
0552   inline void setBinLabels(std::vector<TH2F>& depth) {
0553     // Set labels for all depth histograms
0554     for (unsigned int i = 0; i < depth.size(); ++i) {
0555       depth[i].SetXTitle("i#eta");
0556       depth[i].SetYTitle("i#phi");
0557     }
0558 
0559     std::stringstream label;
0560 
0561     // set label on every other bin
0562     for (int i = -41; i <= -29; i = i + 2) {
0563       label << i;
0564       depth[0].GetXaxis()->SetBinLabel(i + 42, label.str().c_str());
0565       depth[1].GetXaxis()->SetBinLabel(i + 42, label.str().c_str());
0566       label.str("");
0567     }
0568     depth[0].GetXaxis()->SetBinLabel(14, "-29HE");
0569     depth[1].GetXaxis()->SetBinLabel(14, "-29HE");
0570 
0571     // offset by one for HE
0572     for (int i = -27; i <= 27; i = i + 2) {
0573       label << i;
0574       depth[0].GetXaxis()->SetBinLabel(i + 43, label.str().c_str());
0575       label.str("");
0576     }
0577     depth[0].GetXaxis()->SetBinLabel(72, "29HE");
0578     for (int i = 29; i <= 41; i = i + 2) {
0579       label << i;
0580       depth[0].GetXaxis()->SetBinLabel(i + 44, label.str().c_str());
0581       label.str("");
0582     }
0583     for (int i = 16; i <= 28; i = i + 2) {
0584       label << i - 43;
0585       depth[1].GetXaxis()->SetBinLabel(i, label.str().c_str());
0586       label.str("");
0587     }
0588     depth[1].GetXaxis()->SetBinLabel(29, "NULL");
0589     for (int i = 15; i <= 27; i = i + 2) {
0590       label << i;
0591       depth[1].GetXaxis()->SetBinLabel(i + 15, label.str().c_str());
0592       label.str("");
0593     }
0594 
0595     depth[1].GetXaxis()->SetBinLabel(44, "29HE");
0596     for (int i = 29; i <= 41; i = i + 2) {
0597       label << i;
0598       depth[1].GetXaxis()->SetBinLabel(i + 16, label.str().c_str());
0599       label.str("");
0600     }
0601 
0602     // HE depth 3 labels;
0603     depth[2].GetXaxis()->SetBinLabel(1, "-28");
0604     depth[2].GetXaxis()->SetBinLabel(2, "-27");
0605     depth[2].GetXaxis()->SetBinLabel(3, "Null");
0606     depth[2].GetXaxis()->SetBinLabel(4, "-16");
0607     depth[2].GetXaxis()->SetBinLabel(5, "Null");
0608     depth[2].GetXaxis()->SetBinLabel(6, "16");
0609     depth[2].GetXaxis()->SetBinLabel(7, "Null");
0610     depth[2].GetXaxis()->SetBinLabel(8, "27");
0611     depth[2].GetXaxis()->SetBinLabel(9, "28");
0612   }
0613   // Now define functions that can be used in conjunction with EtaPhi histograms
0614 
0615   // This arrays the eta binning for depth 2 histograms (with a gap between -15 -> +15)
0616   const int binmapd2[] = {-42, -41, -40, -39, -38, -37, -36, -35, -34, -33, -32, -31, -30, -29,   -28,
0617                           -27, -26, -25, -24, -23, -22, -21, -20, -19, -18, -17, -16, -15, -9999, 15,
0618                           16,  17,  18,  19,  20,  21,  22,  23,  24,  25,  26,  27,  28,  29,    30,
0619                           31,  32,  33,  34,  35,  36,  37,  38,  39,  40,  41,  42};
0620 
0621   // This stores eta binning in depth 3 (where HE is only present at a few ieta values)
0622 
0623   const int binmapd3[] = {-28, -27, -9999, -16, -9999, 16, -9999, 27, 28};
0624 
0625   inline int CalcEtaBin(int subdet, int ieta, int depth) {
0626     // This takes the eta value from a subdetector and return an eta counter value as used by eta-phi array
0627     // (ieta=-41 corresponds to bin 0, +41 to bin 85 -- there are two offsets to deal with the overlap at |ieta|=29).
0628     // For HO, ieta = -15 corresponds to bin 0, and ieta=15 is bin 30
0629     // For HE depth 3, things are more complicated, but feeding the ieta value will give back the corresponding counter eta value
0630 
0631     // The CalcEtaBin value is the value as used within our array counters, and thus starts at 0.
0632     // If you are using it with getBinContent or setBinContent, you will need to add +1 to the result of this function
0633 
0634     int etabin = -9999;  // default invalid value
0635 
0636     if (depth == 1) {
0637       // Depth 1 is fairly straightforward -- just shift HF-, HF+ by -/+1
0638       etabin = ieta + 42;
0639       if (subdet == HcalForward) {
0640         ieta < 0 ? etabin-- : etabin++;
0641       }
0642     }
0643 
0644     else if (depth == 2) {
0645       // Depth 2 is more complicated, given that there are no cells in the range |ieta|<15
0646       if (ieta < -14) {
0647         etabin = ieta + 42;
0648         if (subdet == HcalForward)
0649           etabin--;
0650       } else if (ieta > 14) {
0651         etabin = ieta + 14;
0652         if (subdet == HcalForward)
0653           etabin++;
0654       }
0655 
0656     }
0657     // HO is also straightforward; a simple offset to the ieta value is applied
0658     else if (subdet == HcalOuter && abs(ieta) < 16)
0659       etabin = ieta + 15;
0660     else if (subdet == HcalEndcap) {
0661       // HE depth 3 has spotty coverage; hard-code the bin response
0662       if (depth == 3) {
0663         if (ieta == -28)
0664           etabin = 0;
0665         else if (ieta == -27)
0666           etabin = 1;
0667         else if (ieta == -16)
0668           etabin = 3;
0669         else if (ieta == 16)
0670           etabin = 5;
0671         else if (ieta == 27)
0672           etabin = 7;
0673         else if (ieta == 28)
0674           etabin = 8;
0675       }
0676     }
0677     return etabin;
0678   }
0679 
0680   inline int CalcIeta(int subdet, int eta, int depth) {
0681     // This function returns the 'true' ieta value given subdet, eta, and depth
0682     // Here 'eta' is the index from our arrays (it starts at 0);
0683     // remember that histogram bins start with bin 1, so there's an offset of 1
0684     // to consider if using getBinContent(eta,phi)
0685 
0686     // eta runs from 0...X  (X depends on depth)
0687     int ieta = -9999;  // default value is nonsensical
0688     if (subdet == HcalBarrel) {
0689       if (depth == 1) {
0690         ieta = eta - 42;
0691         if (ieta == 0)
0692           return -9999;
0693         return ieta;
0694       } else if (depth == 2) {
0695         ieta = binmapd2[eta];
0696         if (ieta == 0)
0697           return -9999;
0698         if (ieta == 17 || ieta == -17)
0699           return -9999;  // no depth 2 cells at |ieta| = 17
0700         return ieta;
0701       } else
0702         return -9999;  // non-physical value
0703     } else if (subdet == HcalForward) {
0704       if (depth == 1) {
0705         ieta = eta - 42;
0706         if (eta < 13)
0707           ieta++;
0708         else if (eta > 71)
0709           ieta--;
0710         else
0711           return -9999;  // if outside forward range, return dummy
0712         return ieta;
0713       } else if (depth == 2) {
0714         ieta = binmapd2[eta];  // special map for depth 2
0715         if (ieta <= -30)
0716           ieta++;
0717         else if (ieta >= 30)
0718           ieta--;
0719         else
0720           return -9999;
0721         return ieta;
0722       } else
0723         return -9999;
0724     }
0725 
0726     else if (subdet == HcalEndcap) {
0727       if (depth == 1)
0728         ieta = eta - 42;
0729       else if (depth == 2) {
0730         ieta = binmapd2[eta];
0731         if (abs(ieta) > 29 || abs(ieta) < 18)
0732           return -9999;  // outside HE
0733         if (ieta == 0)
0734           return -9999;
0735         return ieta;
0736       } else if (depth == 3) {
0737         if (eta < 0 || eta > 8)
0738           return -9999;
0739         else
0740           ieta = binmapd3[eta];  // special map for depth 3
0741         if (ieta == 0)
0742           return -9999;
0743         return ieta;
0744       } else
0745         return -9999;
0746     }  // HcalEndcap
0747     else if (subdet == HcalOuter) {
0748       if (depth != 4)
0749         return -9999;
0750       else {
0751         ieta = eta - 15;  // bin 0 is ieta=-15, all bins increment normally from there
0752         if (abs(ieta) > 15)
0753           return -9999;
0754         if (ieta == 0)
0755           return -9999;
0756         return ieta;
0757       }
0758     }  // HcalOuter
0759     if (ieta == 0)
0760       return -9999;
0761     return ieta;
0762   }
0763 
0764   inline int CalcIeta(int eta, int depth) {
0765     // This version of CalcIeta does the same as the function above,
0766     // but does not require that 'subdet' be specified.
0767 
0768     // returns ieta value give an eta counter.
0769     // eta runs from 0...X  (X depends on depth)
0770     int ieta = -9999;
0771     if (eta < 0)
0772       return ieta;
0773     if (depth == 1) {
0774       ieta =
0775           eta -
0776           42;  // default shift: bin 0 corresponds to a histogram ieta of -42 (which is offset by 1 from true HF value of -41)
0777       if (eta < 13)
0778         ieta++;
0779       else if (eta > 71)
0780         ieta--;
0781       if (ieta == 0)
0782         ieta = -9999;
0783       return ieta;
0784     } else if (depth == 2) {
0785       if (eta > 57)
0786         return -9999;
0787       else {
0788         ieta = binmapd2[eta];
0789         if (ieta == -9999)
0790           return ieta;
0791         if (ieta == 0)
0792           return -9999;
0793         if (ieta == 17 || ieta == -17)
0794           return -9999;  // no depth 2 cells at |ieta| = 17
0795         else if (ieta <= -30)
0796           ieta++;
0797         else if (ieta >= 30)
0798           ieta--;
0799         return ieta;
0800       }
0801     } else if (depth == 3) {
0802       if (eta > 8)
0803         return -9999;
0804       else
0805         ieta = binmapd3[eta];
0806       if (ieta == 0)
0807         return -9999;
0808       return ieta;
0809     } else if (depth == 4) {
0810       ieta = eta - 15;  // bin 0 is ieta=-15, all bins increment normally from there
0811       if (abs(ieta) > 15)
0812         return -9999;
0813       if (ieta == 0)
0814         return -9999;
0815       return ieta;
0816     }
0817     return ieta;  // avoids compilation warning
0818   }
0819 
0820   // Functions to check whether a given (eta,depth) value is valid for a given subdetector
0821 
0822   inline std::vector<std::string> HcalEtaPhiHistNames() {
0823     std::vector<std::string> name;
0824     name.push_back("HB HE HF Depth 1 ");
0825     name.push_back("HB HE HF Depth 2 ");
0826     name.push_back("HE Depth 3 ");
0827     name.push_back("HO Depth 4 ");
0828     return name;
0829   }
0830 
0831   inline bool isHB(int etabin, int depth) {
0832     if (depth > 2)
0833       return false;
0834     else if (depth < 1)
0835       return false;
0836     else {
0837       int ieta = CalcIeta(etabin, depth);
0838       if (ieta == -9999)
0839         return false;
0840       if (depth == 1) {
0841         if (abs(ieta) <= 16)
0842           return true;
0843         else
0844           return false;
0845       } else if (depth == 2) {
0846         if (abs(ieta) == 15 || abs(ieta) == 16)
0847           return true;
0848         else
0849           return false;
0850       }
0851     }
0852     return false;
0853   }
0854 
0855   inline bool isHE(int etabin, int depth) {
0856     if (depth > 3)
0857       return false;
0858     else if (depth < 1)
0859       return false;
0860     else {
0861       int ieta = CalcIeta(etabin, depth);
0862       if (ieta == -9999)
0863         return false;
0864       if (depth == 1) {
0865         if (abs(ieta) >= 17 && abs(ieta) <= 28)
0866           return true;
0867         if (ieta == -29 && etabin == 13)
0868           return true;  // HE -29
0869         if (ieta == 29 && etabin == 71)
0870           return true;  // HE +29
0871       } else if (depth == 2) {
0872         if (abs(ieta) >= 17 && abs(ieta) <= 28)
0873           return true;
0874         if (ieta == -29 && etabin == 13)
0875           return true;  // HE -29
0876         if (ieta == 29 && etabin == 43)
0877           return true;  // HE +29
0878       } else if (depth == 3)
0879         return true;
0880     }
0881     return false;
0882   }
0883 
0884   inline bool isHF(int etabin, int depth) {
0885     if (depth > 2)
0886       return false;
0887     else if (depth < 1)
0888       return false;
0889     else {
0890       int ieta = CalcIeta(etabin, depth);
0891       if (ieta == -9999)
0892         return false;
0893       if (depth == 1) {
0894         if (ieta == -29 && etabin == 13)
0895           return false;  // HE -29
0896         else if (ieta == 29 && etabin == 71)
0897           return false;  // HE +29
0898         else if (abs(ieta) >= 29)
0899           return true;
0900       } else if (depth == 2) {
0901         if (ieta == -29 && etabin == 13)
0902           return false;  // HE -29
0903         else if (ieta == 29 && etabin == 43)
0904           return false;  // HE +29
0905         else if (abs(ieta) >= 29)
0906           return true;
0907       }
0908     }
0909     return false;
0910   }
0911 
0912   inline bool isHO(int etabin, int depth) {
0913     if (depth != 4)
0914       return false;
0915     int ieta = CalcIeta(etabin, depth);
0916     if (ieta != -9999)
0917       return true;
0918     return false;
0919   }
0920 
0921   // Checks whether HO region contains SiPM
0922 
0923   inline bool isSiPM(int ieta, int iphi, int depth) {
0924     if (depth != 4)
0925       return false;
0926     // HOP1
0927     if (ieta >= 5 && ieta <= 10 && iphi >= 47 && iphi <= 58)
0928       return true;
0929     // HOP2
0930     if (ieta >= 11 && ieta <= 15 && iphi >= 59 && iphi <= 70)
0931       return true;
0932     return false;
0933   }  // bool isSiPM
0934 
0935   // Checks whether (subdet, ieta, iphi, depth) value is a valid Hcal cell
0936 
0937   inline bool validDetId(HcalSubdetector sd, int ies, int ip, int dp) {
0938     // inputs are (subdetector, ieta, iphi, depth)
0939     // stolen from latest version of DataFormats/HcalDetId/src/HcalDetId.cc (not yet available in CMSSW_2_1_9)
0940 
0941     const int ie(abs(ies));
0942 
0943     return (
0944         (ip >= 1) && (ip <= 72) && (dp >= 1) && (ie >= 1) &&
0945         (((sd == HcalBarrel) && (((ie <= 14) && (dp == 1)) || (((ie == 15) || (ie == 16)) && (dp <= 2)))) ||
0946          ((sd == HcalEndcap) &&
0947           (((ie == 16) && (dp == 3)) || ((ie == 17) && (dp == 1)) || ((ie >= 18) && (ie <= 20) && (dp <= 2)) ||
0948            ((ie >= 21) && (ie <= 26) && (dp <= 2) && (ip % 2 == 1)) ||
0949            ((ie >= 27) && (ie <= 28) && (dp <= 3) && (ip % 2 == 1)) || ((ie == 29) && (dp <= 2) && (ip % 2 == 1)))) ||
0950          ((sd == HcalOuter) && (ie <= 15) && (dp == 4)) ||
0951          ((sd == HcalForward) && (dp <= 2) &&
0952           (((ie >= 29) && (ie <= 39) && (ip % 2 == 1)) || ((ie >= 40) && (ie <= 41) && (ip % 4 == 3))))));
0953 
0954   }  // bool validDetId(HcalSubdetector sd, int ies, int ip, int dp)
0955 
0956   // Sets eta, phi labels for 'summary' eta-phi plots (identical to Depth 1 Eta-Phi labelling)
0957 
0958   inline void SetEtaPhiLabels(TH2F& h) {
0959     std::stringstream label;
0960     for (int i = -41; i <= -29; i = i + 2) {
0961       label << i;
0962       h.GetXaxis()->SetBinLabel(i + 42, label.str().c_str());
0963       label.str("");
0964     }
0965     h.GetXaxis()->SetBinLabel(14, "-29HE");
0966 
0967     // offset by one for HE
0968     for (int i = -27; i <= 27; i = i + 2) {
0969       label << i;
0970       h.GetXaxis()->SetBinLabel(i + 43, label.str().c_str());
0971       label.str("");
0972     }
0973     h.GetXaxis()->SetBinLabel(72, "29HE");
0974     for (int i = 29; i <= 41; i = i + 2) {
0975       label << i;
0976       h.GetXaxis()->SetBinLabel(i + 44, label.str().c_str());
0977       label.str("");
0978     }
0979     return;
0980   }
0981 
0982   // Fill Unphysical bins in histograms
0983   inline void FillUnphysicalHEHFBins(std::vector<TH2F>& hh) {
0984     int ieta = 0;
0985     int iphi = 0;
0986     // First 2 depths have 5-10-20 degree corrections
0987     for (unsigned int d = 0; d < 3; ++d) {
0988       //BUG CAN BE HERE:
0989       //if (hh[d] != 0) continue;
0990 
0991       for (int eta = 0; eta < hh[d].GetNbinsX(); ++eta) {
0992         ieta = CalcIeta(eta, d + 1);
0993         if (ieta == -9999 || abs(ieta) < 21)
0994           continue;
0995         for (int phi = 0; phi < hh[d].GetNbinsY(); ++phi) {
0996           iphi = phi + 1;
0997           if (iphi % 2 == 1 && abs(ieta) < 40 && iphi < 73) {
0998             hh[d].SetBinContent(eta + 1, iphi + 1, hh[d].GetBinContent(eta + 1, iphi));
0999           }
1000           // last two eta strips span 20 degrees in phi
1001           // Fill the phi cell above iphi, and the 2 below it
1002           else if (abs(ieta) > 39 && iphi % 4 == 3 && iphi < 73) {
1003             //ieta=40, iphi=3 covers iphi 3,4,5,6
1004             hh[d].SetBinContent(eta + 1, (iphi) % 72 + 1, hh[d].GetBinContent(eta + 1, iphi));
1005             hh[d].SetBinContent(eta + 1, (iphi + 1) % 72 + 1, hh[d].GetBinContent(eta + 1, iphi));
1006             hh[d].SetBinContent(eta + 1, (iphi + 2) % 72 + 1, hh[d].GetBinContent(eta + 1, iphi));
1007           }
1008         }  // for (int phi...)
1009       }  // for (int eta...)
1010     }  // for (int d=0;...)
1011     // no corrections needed for HO (depth 4)
1012     return;
1013   }  // FillUnphysicalHEHFBins(MonitorElement* hh)
1014 
1015   //Fill unphysical bins for single ME
1016   inline void FillUnphysicalHEHFBins(TH2F& hh) {
1017     // Fills unphysical HE/HF bins for Summary Histogram
1018     // Summary Histogram is binned with the same binning as the Depth 1 EtaPhiHists
1019 
1020     //CAN BE BUG HERE:
1021     //if (hh==0) return;
1022 
1023     int ieta = 0;
1024     int iphi = 0;
1025     int etabins = hh.GetNbinsX();
1026     int phibins = hh.GetNbinsY();
1027     float binval = 0;
1028     for (int eta = 0; eta < etabins; ++eta)  // loop over eta bins
1029     {
1030       ieta = CalcIeta(eta, 1);
1031       if (ieta == -9999 || abs(ieta) < 21)
1032         continue;  // ignore etas that don't exist, or that have 5 degree phi binning
1033 
1034       for (int phi = 0; phi < phibins; ++phi) {
1035         iphi = phi + 1;
1036         if (iphi % 2 == 1 && abs(ieta) < 40 && iphi < 73)  // 10 degree phi binning condition
1037         {
1038           binval = hh.GetBinContent(eta + 1, iphi);
1039           hh.SetBinContent(eta + 1, iphi + 1, binval);
1040         }  // if (iphi%2==1...)
1041         else if (abs(ieta) > 39 && iphi % 4 == 3 && iphi < 73)  // 20 degree phi binning condition
1042         {
1043           // Set last two eta strips where each cell spans 20 degrees in phi
1044           // Set next phi cell above iphi, and 2 cells below the actual cell
1045           hh.SetBinContent(eta + 1, (iphi) % 72 + 1, hh.GetBinContent(eta + 1, iphi));
1046           hh.SetBinContent(eta + 1, (iphi + 1) % 72 + 1, hh.GetBinContent(eta + 1, iphi));
1047           hh.SetBinContent(eta + 1, (iphi + 2) % 72 + 1, hh.GetBinContent(eta + 1, iphi));
1048         }  // else if (abs(ieta)>39 ...)
1049       }  // for (int phi=0;phi<72;++phi)
1050 
1051     }  // for (int eta=0; eta< (etaBins_-2);++eta)
1052 
1053     return;
1054   }  // FillUnphysicalHEHFBins(std::vector<MonitorElement*> &hh)
1055 
1056   // special fill call based on detid -- eventually will need special treatment
1057   inline void Fill(HcalDetId& id, double val /*=1*/, std::vector<TH2F>& depth) {
1058     // If in HF, need to shift by 1 bin (-1 bin lower in -HF, +1 bin higher in +HF)
1059     if (id.subdet() == HcalForward)
1060       depth[id.depth() - 1].Fill(id.ieta() < 0 ? id.ieta() - 1 : id.ieta() + 1, id.iphi(), val);
1061     else
1062       depth[id.depth() - 1].Fill(id.ieta(), id.iphi(), val);
1063   }
1064 
1065   inline void Reset(std::vector<TH2F>& depth) {
1066     for (unsigned int d = 0; d < depth.size(); d++)
1067       //BUG CAN BE HERE:
1068       //if(depth[d])
1069       depth[d].Reset();
1070   }  // void Reset(void)
1071 
1072   inline void setup(std::vector<TH2F>& depth, std::string name, std::string units = "") {
1073     std::string unittitle, unitname;
1074     if (units.empty()) {
1075       unitname = units;
1076       unittitle = "No Units";
1077     } else {
1078       unitname = " " + units;
1079       unittitle = units;
1080     }
1081 
1082     // Push back depth plots
1083     ////////Create 4 plots:
1084     //1. create first plot
1085     depth.push_back(TH2F(("HB HE HF Depth 1 " + name + unitname).c_str(),
1086                          (name + " Depth 1 -- HB HE HF (" + unittitle + ")").c_str(),
1087                          85,
1088                          -42.5,
1089                          42.5,
1090                          72,
1091                          0.5,
1092                          72.5));
1093 
1094     //2.1 prepare second plot
1095     float ybins[73];
1096     for (int i = 0; i <= 72; i++)
1097       ybins[i] = (float)(i + 0.5);
1098     float xbinsd2[] = {-42.5, -41.5, -40.5, -39.5, -38.5, -37.5, -36.5, -35.5, -34.5, -33.5, -32.5, -31.5,
1099                        -30.5, -29.5, -28.5, -27.5, -26.5, -25.5, -24.5, -23.5, -22.5, -21.5, -20.5, -19.5,
1100                        -18.5, -17.5, -16.5, -15.5, -14.5, 14.5,  15.5,  16.5,  17.5,  18.5,  19.5,  20.5,
1101                        21.5,  22.5,  23.5,  24.5,  25.5,  26.5,  27.5,  28.5,  29.5,  30.5,  31.5,  32.5,
1102                        33.5,  34.5,  35.5,  36.5,  37.5,  38.5,  39.5,  40.5,  41.5,  42.5};
1103 
1104     //2.2 create second plot
1105     depth.push_back(TH2F(("HB HE HF Depth 2 " + name + unitname).c_str(),
1106                          (name + " Depth 2 -- HB HE HF (" + unittitle + ")").c_str(),
1107                          57,
1108                          xbinsd2,
1109                          72,
1110                          ybins));
1111 
1112     //3.1 Set up variable-sized bins for HE depth 3 (MonitorElement also requires phi bins to be entered in array format)
1113     float xbins[] = {-28.5, -27.5, -26.5, -16.5, -15.5, 15.5, 16.5, 26.5, 27.5, 28.5};
1114     //3.2
1115     depth.push_back(TH2F(("HE Depth 3 " + name + unitname).c_str(),
1116                          (name + " Depth 3 -- HE (" + unittitle + ")").c_str(),
1117                          // Use variable-sized eta bins
1118                          9,
1119                          xbins,
1120                          72,
1121                          ybins));
1122 
1123     //4.1 HO bins are fixed width, but cover a smaller eta range (-15 -> 15)
1124     depth.push_back(TH2F(("HO Depth 4 " + name + unitname).c_str(),
1125                          (name + " Depth 4 -- HO (" + unittitle + ")").c_str(),
1126                          31,
1127                          -15.5,
1128                          15.5,
1129                          72,
1130                          0.5,
1131                          72.5));
1132 
1133     for (unsigned int i = 0; i < depth.size(); ++i)
1134       depth[i].Draw("colz");
1135 
1136     setBinLabels(depth);  // set axis titles, special bins
1137   }
1138 
1139   inline void fillOneGain(std::vector<TH2F>& graphData,
1140                           HcalGains::tAllContWithNames& allContainers,
1141                           std::string name,
1142                           int id,
1143                           std::string units = "") {
1144     setup(graphData, name);
1145 
1146     std::stringstream x;
1147     // Change the titles of each individual histogram
1148     for (unsigned int d = 0; d < graphData.size(); ++d) {
1149       graphData[d].Reset();
1150       x << "Gain " << id << " for HCAL depth " << d + 1;
1151 
1152       //BUG CAN BE HERE:
1153       //if (ChannelStatus->depth[d])
1154       graphData[d].SetTitle(
1155           x.str()
1156               .c_str());  // replace "setTitle" with "SetTitle", since you are using TH2F objects instead of MonitorElements
1157       x.str("");
1158     }
1159 
1160     HcalDetId hcal_id;
1161     int ieta, depth, iphi;
1162 
1163     //main loop
1164     // get all containers with names
1165     //HcalGains::tAllContWithNames allContainers = object().getAllContainers();
1166 
1167     //ITERATORS AND VALUES:
1168     HcalGains::tAllContWithNames::const_iterator iter;
1169     std::vector<HcalGain>::const_iterator contIter;
1170     float gain = 0.0;
1171 
1172     //Run trough given id gain:
1173     int i = id;
1174 
1175     //run trough all pair containers
1176     for (iter = allContainers.begin(); iter != allContainers.end(); ++iter) {
1177       //Run trough all values:
1178       for (contIter = (*iter).second.begin(); contIter != (*iter).second.end(); ++contIter) {
1179         hcal_id = HcalDetId((uint32_t)(*contIter).rawId());
1180 
1181         depth = hcal_id.depth();
1182         if (depth < 1 || depth > 4)
1183           continue;
1184 
1185         ieta = hcal_id.ieta();
1186         iphi = hcal_id.iphi();
1187 
1188         if (hcal_id.subdet() == HcalForward)
1189           ieta > 0 ? ++ieta : --ieta;
1190 
1191         //GET VALUE:
1192         gain = (*contIter).getValue(i);
1193         //logstatus = log2(1.*channelBits)+1;
1194 
1195         //FILLING GOES HERE:
1196         graphData[depth - 1].Fill(ieta, iphi, gain);
1197 
1198         //FOR DEBUGGING:
1199         //std::cout << "ieta: " << ieta << "; iphi: " << iphi << "; logstatus: " << logstatus << "; channelBits: " << channelBits<< std::endl;
1200       }
1201     }
1202   }
1203 
1204   //FillUnphysicalHEHFBins(graphData);
1205   //return ("kasdasd");
1206 
1207   class ADataRepr  //Sample base class for c++ inheritance tutorial
1208   {
1209   public:
1210     ADataRepr(unsigned int d) : m_total(d) {}
1211     virtual ~ADataRepr() {}
1212     unsigned int nr, id;
1213     std::stringstream filename, rootname, plotname;
1214 
1215     void fillOneGain(std::vector<TH2F>& graphData, std::string units = "") {
1216       std::stringstream ss("");
1217 
1218       if (m_total == 1)
1219         ss << rootname.str() << " for HCAL depth ";
1220       else
1221         ss << rootname.str() << nr << " for HCAL depth ";
1222 
1223       setup(graphData, ss.str());
1224 
1225       // Change the titles of each individual histogram
1226       for (unsigned int d = 0; d < graphData.size(); ++d) {
1227         graphData[d].Reset();
1228 
1229         ss.str("");
1230         if (m_total == 1)
1231           ss << plotname.str() << " for HCAL depth " << d + 1;
1232         else
1233           ss << plotname.str() << nr << " for HCAL depth " << d + 1;
1234 
1235         //BUG CAN BE HERE:
1236         //if (ChannelStatus->depth[d])
1237         graphData[d].SetTitle(
1238             ss.str()
1239                 .c_str());  // replace "setTitle" with "SetTitle", since you are using TH2F objects instead of MonitorElements
1240         ss.str("");
1241       }
1242       //overload this function:
1243       doFillIn(graphData);
1244 
1245       FillUnphysicalHEHFBins(graphData);
1246 
1247       ss.str("");
1248       if (m_total == 1)
1249         ss << filename.str() << ".png";
1250       else
1251         ss << filename.str() << nr << ".png";
1252       draw(graphData, ss.str());
1253       //FOR DEBUGGING:
1254       //std::cout << "ieta: " << ieta << "; iphi: " << iphi << "; logstatus: " << logstatus << "; channelBits: " << channelBits<< std::endl;
1255     }
1256 
1257   protected:
1258     unsigned int m_total;
1259     HcalDetId hcal_id;
1260     int ieta, depth, iphi;
1261 
1262     virtual void doFillIn(std::vector<TH2F>& graphData) = 0;
1263 
1264   private:
1265     void draw(std::vector<TH2F>& graphData, std::string filename) {
1266       //Drawing...
1267       // use David's palette
1268       gStyle->SetPalette(1);
1269       const Int_t NCont = 999;
1270       gStyle->SetNumberContours(NCont);
1271       TCanvas canvas("CC map", "CC map", 840, 369 * 4);
1272 
1273       TPad pad1("pad1", "pad1", 0.0, 0.75, 1.0, 1.0);
1274       pad1.Draw();
1275       TPad pad2("pad2", "pad2", 0.0, 0.5, 1.0, 0.75);
1276       pad2.Draw();
1277       TPad pad3("pad3", "pad3", 0.0, 0.25, 1.0, 0.5);
1278       pad3.Draw();
1279       TPad pad4("pad4", "pad4", 0.0, 0.0, 1.0, 0.25);
1280       pad4.Draw();
1281 
1282       pad1.cd();
1283       graphData[0].SetStats(false);
1284       graphData[0].Draw("colz");
1285 
1286       pad2.cd();
1287       graphData[1].SetStats(false);
1288       graphData[1].Draw("colz");
1289 
1290       pad3.cd();
1291       graphData[2].SetStats(false);
1292       graphData[2].Draw("colz");
1293 
1294       pad4.cd();
1295       graphData[3].SetStats(false);
1296       graphData[3].Draw("colz");
1297 
1298       canvas.SaveAs(filename.c_str());
1299     }
1300 
1301     void setup(std::vector<TH2F>& depth, std::string name, std::string units = "") {
1302       std::string unittitle, unitname;
1303       if (units.empty()) {
1304         unitname = units;
1305         unittitle = "No Units";
1306       } else {
1307         unitname = " " + units;
1308         unittitle = units;
1309       }
1310 
1311       // Push back depth plots
1312       ////////Create 4 plots:
1313       //1. create first plot
1314       depth.push_back(TH2F(("HB HE HF Depth 1 " + name + unitname).c_str(),
1315                            (name + " Depth 1 -- HB HE HF (" + unittitle + ")").c_str(),
1316                            85,
1317                            -42.5,
1318                            42.5,
1319                            72,
1320                            0.5,
1321                            72.5));
1322 
1323       //2.1 prepare second plot
1324       float ybins[73];
1325       for (int i = 0; i <= 72; i++)
1326         ybins[i] = (float)(i + 0.5);
1327       float xbinsd2[] = {-42.5, -41.5, -40.5, -39.5, -38.5, -37.5, -36.5, -35.5, -34.5, -33.5, -32.5, -31.5,
1328                          -30.5, -29.5, -28.5, -27.5, -26.5, -25.5, -24.5, -23.5, -22.5, -21.5, -20.5, -19.5,
1329                          -18.5, -17.5, -16.5, -15.5, -14.5, 14.5,  15.5,  16.5,  17.5,  18.5,  19.5,  20.5,
1330                          21.5,  22.5,  23.5,  24.5,  25.5,  26.5,  27.5,  28.5,  29.5,  30.5,  31.5,  32.5,
1331                          33.5,  34.5,  35.5,  36.5,  37.5,  38.5,  39.5,  40.5,  41.5,  42.5};
1332 
1333       //2.2 create second plot
1334       depth.push_back(TH2F(("HB HE HF Depth 2 " + name + unitname).c_str(),
1335                            (name + " Depth 2 -- HB HE HF (" + unittitle + ")").c_str(),
1336                            57,
1337                            xbinsd2,
1338                            72,
1339                            ybins));
1340 
1341       //3.1 Set up variable-sized bins for HE depth 3 (MonitorElement also requires phi bins to be entered in array format)
1342       float xbins[] = {-28.5, -27.5, -26.5, -16.5, -15.5, 15.5, 16.5, 26.5, 27.5, 28.5};
1343       //3.2
1344       depth.push_back(TH2F(("HE Depth 3 " + name + unitname).c_str(),
1345                            (name + " Depth 3 -- HE (" + unittitle + ")").c_str(),
1346                            // Use variable-sized eta bins
1347                            9,
1348                            xbins,
1349                            72,
1350                            ybins));
1351 
1352       //4.1 HO bins are fixed width, but cover a smaller eta range (-15 -> 15)
1353       depth.push_back(TH2F(("HO Depth 4 " + name + unitname).c_str(),
1354                            (name + " Depth 4 -- HO (" + unittitle + ")").c_str(),
1355                            31,
1356                            -15.5,
1357                            15.5,
1358                            72,
1359                            0.5,
1360                            72.5));
1361 
1362       for (unsigned int i = 0; i < depth.size(); ++i)
1363         depth[i].Draw("colz");
1364 
1365       setBinLabels(depth);  // set axis titles, special bins
1366     }
1367 
1368     //functions for making plot:
1369     void setBinLabels(std::vector<TH2F>& depth) {
1370       // Set labels for all depth histograms
1371       for (unsigned int i = 0; i < depth.size(); ++i) {
1372         depth[i].SetXTitle("i#eta");
1373         depth[i].SetYTitle("i#phi");
1374       }
1375 
1376       std::stringstream label;
1377 
1378       // set label on every other bin
1379       for (int i = -41; i <= -29; i = i + 2) {
1380         label << i;
1381         depth[0].GetXaxis()->SetBinLabel(i + 42, label.str().c_str());
1382         depth[1].GetXaxis()->SetBinLabel(i + 42, label.str().c_str());
1383         label.str("");
1384       }
1385       depth[0].GetXaxis()->SetBinLabel(14, "-29HE");
1386       depth[1].GetXaxis()->SetBinLabel(14, "-29HE");
1387 
1388       // offset by one for HE
1389       for (int i = -27; i <= 27; i = i + 2) {
1390         label << i;
1391         depth[0].GetXaxis()->SetBinLabel(i + 43, label.str().c_str());
1392         label.str("");
1393       }
1394       depth[0].GetXaxis()->SetBinLabel(72, "29HE");
1395       for (int i = 29; i <= 41; i = i + 2) {
1396         label << i;
1397         depth[0].GetXaxis()->SetBinLabel(i + 44, label.str().c_str());
1398         label.str("");
1399       }
1400       for (int i = 16; i <= 28; i = i + 2) {
1401         label << i - 43;
1402         depth[1].GetXaxis()->SetBinLabel(i, label.str().c_str());
1403         label.str("");
1404       }
1405       depth[1].GetXaxis()->SetBinLabel(29, "NULL");
1406       for (int i = 15; i <= 27; i = i + 2) {
1407         label << i;
1408         depth[1].GetXaxis()->SetBinLabel(i + 15, label.str().c_str());
1409         label.str("");
1410       }
1411 
1412       depth[1].GetXaxis()->SetBinLabel(44, "29HE");
1413       for (int i = 29; i <= 41; i = i + 2) {
1414         label << i;
1415         depth[1].GetXaxis()->SetBinLabel(i + 16, label.str().c_str());
1416         label.str("");
1417       }
1418 
1419       // HE depth 3 labels;
1420       depth[2].GetXaxis()->SetBinLabel(1, "-28");
1421       depth[2].GetXaxis()->SetBinLabel(2, "-27");
1422       depth[2].GetXaxis()->SetBinLabel(3, "Null");
1423       depth[2].GetXaxis()->SetBinLabel(4, "-16");
1424       depth[2].GetXaxis()->SetBinLabel(5, "Null");
1425       depth[2].GetXaxis()->SetBinLabel(6, "16");
1426       depth[2].GetXaxis()->SetBinLabel(7, "Null");
1427       depth[2].GetXaxis()->SetBinLabel(8, "27");
1428       depth[2].GetXaxis()->SetBinLabel(9, "28");
1429     }
1430   };
1431 }  // namespace HcalObjRepresent
1432 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!