Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 11:59:40

0001 // -*- C++ -*-
0002 // Package:    SiStripChannelGain
0003 // Class:      SiStripGainCosmicCalculator
0004 // Original Author:  G. Bruno, D. Kcira
0005 //         Created:  Mon May 20 10:04:31 CET 2007
0006 #include "CalibTracker/SiStripChannelGain/plugins/SiStripGainCosmicCalculator.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0010 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0011 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0012 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0013 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetType.h"
0014 #include "DataFormats/FEDRawData/interface/FEDNumbering.h"
0015 #include "CLHEP/Random/RandFlat.h"
0016 #include "CLHEP/Random/RandGauss.h"
0017 #include "Geometry/CommonDetUnit/interface/PixelGeomDetType.h"
0018 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0019 #include "DataFormats/TrackReco/interface/Track.h"
0020 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0021 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
0022 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
0023 #include "DataFormats/TrackerCommon/interface/SiStripSubStructure.h"
0024 #include "DataFormats/DetId/interface/DetId.h"
0025 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0026 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0027 //#include "DQM/SiStripCommon/interface/SiStripGenerateKey.h"
0028 
0029 //---------------------------------------------------------------------------------------------------------
0030 SiStripGainCosmicCalculator::SiStripGainCosmicCalculator(const edm::ParameterSet& iConfig)
0031     : ConditionDBWriter<SiStripApvGain>(iConfig) {
0032   edm::LogInfo("SiStripGainCosmicCalculator::SiStripGainCosmicCalculator");
0033   ExpectedChargeDeposition = 200.;
0034   edm::LogInfo("SiStripApvGainCalculator::SiStripApvGainCalculator")
0035       << "ExpectedChargeDeposition=" << ExpectedChargeDeposition;
0036 
0037   TrackProducer = iConfig.getParameter<std::string>("TrackProducer");
0038   TrackLabel = iConfig.getParameter<std::string>("TrackLabel");
0039 
0040   detModulesToBeExcluded.clear();
0041   detModulesToBeExcluded = iConfig.getParameter<std::vector<unsigned> >("detModulesToBeExcluded");
0042   MinNrEntries = iConfig.getUntrackedParameter<unsigned>("minNrEntries", 20);
0043   MaxChi2OverNDF = iConfig.getUntrackedParameter<double>("maxChi2OverNDF", 5.);
0044 
0045   outputHistogramsInRootFile = iConfig.getParameter<bool>("OutputHistogramsInRootFile");
0046   outputFileName = iConfig.getParameter<std::string>("OutputFileName");
0047 
0048   edm::LogInfo("SiStripApvGainCalculator")
0049       << "Clusters from " << detModulesToBeExcluded.size() << " modules will be ignored in the calibration:";
0050   edm::LogInfo("SiStripApvGainCalculator") << "The calibration for these DetIds will be set to a default value";
0051   for (std::vector<uint32_t>::const_iterator imod = detModulesToBeExcluded.begin();
0052        imod != detModulesToBeExcluded.end();
0053        imod++) {
0054     edm::LogInfo("SiStripApvGainCalculator") << "exclude detid = " << *imod;
0055   }
0056 
0057   printdebug_ = iConfig.getUntrackedParameter<bool>("printDebug", false);
0058 
0059   tTopoToken_ = esConsumes<edm::Transition::BeginRun>();
0060   detCablingToken_ = esConsumes<edm::Transition::BeginRun>();
0061   tkGeomToken_ = esConsumes<edm::Transition::BeginRun>();
0062 }
0063 
0064 SiStripGainCosmicCalculator::~SiStripGainCosmicCalculator() {
0065   edm::LogInfo("SiStripGainCosmicCalculator::~SiStripGainCosmicCalculator");
0066 }
0067 
0068 void SiStripGainCosmicCalculator::algoEndJob() {}
0069 
0070 void SiStripGainCosmicCalculator::algoBeginJob(const edm::EventSetup& iSetup) {
0071   tTopo_ = &iSetup.getData(tTopoToken_);
0072   siStripDetCabling_ = &iSetup.getData(detCablingToken_);
0073   tkGeom_ = &iSetup.getData(tkGeomToken_);
0074 
0075   std::cout << "SiStripGainCosmicCalculator::algoBeginJob called" << std::endl;
0076   total_nr_of_events = 0;
0077   HlistAPVPairs = new TObjArray();
0078   HlistOtherHistos = new TObjArray();
0079   //
0080   HlistOtherHistos->Add(new TH1F(Form("APVPairCorrections"), Form("APVPairCorrections"), 50, -1., 4.));
0081   HlistOtherHistos->Add(new TH1F(Form("APVPairCorrectionsTIB1mono"), Form("APVPairCorrectionsTIB1mono"), 50, -1., 4.));
0082   HlistOtherHistos->Add(
0083       new TH1F(Form("APVPairCorrectionsTIB1stereo"), Form("APVPairCorrectionsTIB1stereo"), 50, -1., 4.));
0084   HlistOtherHistos->Add(new TH1F(Form("APVPairCorrectionsTIB2"), Form("APVPairCorrectionsTIB2"), 50, -1., 4.));
0085   HlistOtherHistos->Add(new TH1F(Form("APVPairCorrectionsTOB1"), Form("APVPairCorrectionsTOB1"), 50, -1., 4.));
0086   HlistOtherHistos->Add(new TH1F(Form("APVPairCorrectionsTOB2"), Form("APVPairCorrectionsTOB2"), 50, -1., 4.));
0087   HlistOtherHistos->Add(new TH1F(Form("LocalAngle"), Form("LocalAngle"), 70, -0.1, 3.4));
0088   HlistOtherHistos->Add(new TH1F(Form("LocalAngleAbsoluteCosine"), Form("LocalAngleAbsoluteCosine"), 48, -0.1, 1.1));
0089   HlistOtherHistos->Add(new TH1F(Form("LocalPosition_cm"), Form("LocalPosition_cm"), 100, -5., 5.));
0090   HlistOtherHistos->Add(new TH1F(Form("LocalPosition_normalized"), Form("LocalPosition_normalized"), 100, -1.1, 1.1));
0091   TH1F* local_histo = new TH1F(Form("SiStripRecHitType"), Form("SiStripRecHitType"), 2, 0.5, 2.5);
0092   HlistOtherHistos->Add(local_histo);
0093   local_histo->GetXaxis()->SetBinLabel(1, "simple");
0094   local_histo->GetXaxis()->SetBinLabel(2, "matched");
0095 
0096   // get cabling and find out list of active detectors
0097   std::vector<uint32_t> activeDets;
0098   activeDets.clear();
0099   SelectedDetIds.clear();
0100   siStripDetCabling_->addActiveDetectorsRawIds(activeDets);
0101   //    SelectedDetIds = activeDets; // all active detector modules
0102   // use SiStripSubStructure for selecting certain regions
0103   SiStripSubStructure::getTIBDetectors(
0104       activeDets, SelectedDetIds, tTopo_, 0, 0, 0, 0);  // this adds rawDetIds to SelectedDetIds
0105   SiStripSubStructure::getTOBDetectors(
0106       activeDets, SelectedDetIds, tTopo_, 0, 0, 0);  // this adds rawDetIds to SelectedDetIds
0107   // get tracker geometry and find nr. of apv pairs for each active detector
0108   for (auto det : tkGeom_->dets()) {  // loop over detector modules
0109     if (dynamic_cast<const StripGeomDetUnit*>(det) != nullptr) {
0110       uint32_t detid = det->geographicalId().rawId();
0111       // get thickness for all detector modules, not just for active, this is strange
0112       double module_thickness =
0113           det->surface()
0114               .bounds()
0115               .thickness();  // get thickness of detector from GeomDet (DetContainer == vector<GeomDet*>)
0116       thickness_map.insert(std::make_pair(detid, module_thickness));
0117       //
0118       const bool is_active_detector =
0119           std::end(SelectedDetIds) != std::find(std::begin(SelectedDetIds), std::end(SelectedDetIds), detid);
0120       //
0121       const bool exclude_this_detid =
0122           std::end(detModulesToBeExcluded) !=
0123           std::find(std::begin(detModulesToBeExcluded), std::end(detModulesToBeExcluded), detid);
0124       //
0125       if (is_active_detector &&
0126           (!exclude_this_detid)) {  // check whether is active detector and that should not be excluded
0127         const StripTopology& p = dynamic_cast<const StripGeomDetUnit*>(det)->specificTopology();
0128         unsigned short NAPVPairs = p.nstrips() / 256;
0129         if (NAPVPairs < 2 || NAPVPairs > 3) {
0130           edm::LogError("SiStripGainCosmicCalculator")
0131               << "Problem with Number of strips in detector: " << p.nstrips() << " Exiting program";
0132           exit(1);
0133         }
0134         for (int iapp = 0; iapp < NAPVPairs; iapp++) {
0135           TString hid = Form("ChargeAPVPair_%i_%i", detid, iapp);
0136           HlistAPVPairs->Add(
0137               new TH1F(hid, hid, 45, 0., 1350.));  // multiply by 3 to take into account division by width
0138         }
0139       }
0140     }
0141   }
0142 }
0143 
0144 //---------------------------------------------------------------------------------------------------------
0145 void SiStripGainCosmicCalculator::algoAnalyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0146   using namespace edm;
0147   total_nr_of_events++;
0148 
0149   //TO BE RESTORED
0150   //  anglefinder_->init(event,iSetup);
0151 
0152   // get seeds
0153   //  edm::Handle<TrajectorySeedCollection> seedcoll;
0154   //  event.getByType(seedcoll);
0155   // get tracks
0156   Handle<reco::TrackCollection> trackCollection;
0157   iEvent.getByLabel(TrackProducer, TrackLabel, trackCollection);
0158   const reco::TrackCollection* tracks = trackCollection.product();
0159 
0160   //  // get magnetic field
0161   //  edm::ESHandle<MagneticField> esmagfield;
0162   //  es.get<IdealMagneticFieldRecord>().get(esmagfield);
0163   //  magfield=&(*esmagfield);
0164   // loop over tracks
0165   for (reco::TrackCollection::const_iterator itr = tracks->begin(); itr != tracks->end();
0166        itr++) {  // looping over tracks
0167 
0168     //TO BE RESTORED
0169     //    std::vector<std::pair<const TrackingRecHit *,float> >hitangle =anglefinder_->findtrackangle((*(*seedcoll).begin()),*itr);
0170     std::vector<std::pair<const TrackingRecHit*, float> >
0171         hitangle;  // =anglefinder_->findtrackangle((*(*seedcoll).begin()),*itr);
0172 
0173     for (std::vector<std::pair<const TrackingRecHit*, float> >::const_iterator hitangle_iter = hitangle.begin();
0174          hitangle_iter != hitangle.end();
0175          hitangle_iter++) {
0176       const TrackingRecHit* trechit = hitangle_iter->first;
0177       float local_angle = hitangle_iter->second;
0178       LocalPoint local_position = trechit->localPosition();
0179       const SiStripRecHit2D* sistripsimplehit = dynamic_cast<const SiStripRecHit2D*>(trechit);
0180       const SiStripMatchedRecHit2D* sistripmatchedhit = dynamic_cast<const SiStripMatchedRecHit2D*>(trechit);
0181       //      std::cout<<" hit/matched "<<std::ios::hex<<sistripsimplehit<<" "<<sistripmatchedhit<<std::endl;
0182       ((TH1F*)HlistOtherHistos->FindObject("LocalAngle"))->Fill(local_angle);
0183       ((TH1F*)HlistOtherHistos->FindObject("LocalAngleAbsoluteCosine"))->Fill(fabs(cos(local_angle)));
0184       if (sistripsimplehit) {
0185         ((TH1F*)HlistOtherHistos->FindObject("SiStripRecHitType"))->Fill(1.);
0186         const SiStripRecHit2D::ClusterRef& cluster = sistripsimplehit->cluster();
0187         const auto& ampls = cluster->amplitudes();
0188         uint32_t thedetid = 0;  // is zero since long time cluster->geographicalId();
0189         double module_width = moduleWidth(thedetid);
0190         ((TH1F*)HlistOtherHistos->FindObject("LocalPosition_cm"))->Fill(local_position.x());
0191         ((TH1F*)HlistOtherHistos->FindObject("LocalPosition_normalized"))->Fill(local_position.x() / module_width);
0192         double module_thickness = moduleThickness(thedetid);
0193         int ifirststrip = cluster->firstStrip();
0194         int theapvpairid = int(float(ifirststrip) / 256.);
0195         TH1F* histopointer = (TH1F*)HlistAPVPairs->FindObject(Form("ChargeAPVPair_%i_%i", thedetid, theapvpairid));
0196         if (histopointer) {
0197           short cCharge = 0;
0198           for (unsigned int iampl = 0; iampl < ampls.size(); iampl++) {
0199             cCharge += ampls[iampl];
0200           }
0201           double cluster_charge_over_path = ((double)cCharge) * fabs(cos(local_angle)) / (10. * module_thickness);
0202           histopointer->Fill(cluster_charge_over_path);
0203         }
0204       } else {
0205         if (sistripmatchedhit)
0206           ((TH1F*)HlistOtherHistos->FindObject("SiStripRecHitType"))->Fill(2.);
0207       }
0208     }
0209   }
0210 }
0211 
0212 //---------------------------------------------------------------------------------------------------------
0213 std::pair<double, double> SiStripGainCosmicCalculator::getPeakOfLandau(
0214     TH1F* inputHisto) {  // automated fitting with finding of the appropriate nr. of ADCs
0215   // set some default dummy value and return if no entries
0216   double adcs = -0.5;
0217   double error = 0.;
0218   double nr_of_entries = inputHisto->GetEntries();
0219   if (nr_of_entries < MinNrEntries) {
0220     return std::make_pair(adcs, error);
0221   }
0222   //
0223   //  // fit with initial setting of  parameter values
0224   //  double rms_of_histogram = inputHisto->GetRMS();
0225   //  TF1 *landaufit = new TF1("landaufit","landau",0.,450.);
0226   //  landaufit->SetParameters(nr_of_entries,mean_of_histogram,rms_of_histogram);
0227   //  inputHisto->Fit("landaufit","0Q+");
0228   //  delete landaufit;
0229   //
0230   // perform fit with standard landau
0231   // make our own copy to avoid problems with threads
0232   std::unique_ptr<TF1> fitfunction(new TF1("landaufit", "landau"));
0233   inputHisto->Fit(fitfunction.get(), "0Q");
0234   adcs = fitfunction->GetParameter("MPV");
0235   error = fitfunction->GetParError(1);  // MPV is parameter 1 (0=constant, 1=MPV, 2=Sigma)
0236   double chi2 = fitfunction->GetChisquare();
0237   double ndf = fitfunction->GetNDF();
0238   double chi2overndf = chi2 / ndf;
0239   // in case things went wrong, try to refit in smaller range
0240   if (adcs < 2. || (error / adcs) > 1.8) {
0241     inputHisto->Fit(fitfunction.get(), "0Q", nullptr, 0., 400.);
0242     std::cout << "refitting landau for histogram " << inputHisto->GetTitle() << std::endl;
0243     std::cout << "initial error/adcs =" << error << " / " << adcs << std::endl;
0244     std::cout << "new     error/adcs =" << fitfunction->GetParError(1) << " / " << fitfunction->GetParameter("MPV")
0245               << std::endl;
0246     adcs = fitfunction->GetParameter("MPV");
0247     error = fitfunction->GetParError(1);  // MPV is parameter 1 (0=constant, 1=MPV, 2=Sigma)
0248     chi2 = fitfunction->GetChisquare();
0249     ndf = fitfunction->GetNDF();
0250     chi2overndf = chi2 / ndf;
0251   }
0252   // if still wrong, give up
0253   if (adcs < 2. || chi2overndf > MaxChi2OverNDF) {
0254     adcs = -0.5;
0255     error = 0.;
0256   }
0257   return std::make_pair(adcs, error);
0258 }
0259 
0260 //---------------------------------------------------------------------------------------------------------
0261 double SiStripGainCosmicCalculator::moduleWidth(const uint32_t detid)  // get width of the module detid
0262 {  //dk: copied from A. Giammanco and hacked,  module_width values : 10.49 12.03 6.144 7.14 9.3696
0263   double module_width = 0.;
0264   const GeomDetUnit* it = tkGeom_->idToDetUnit(DetId(detid));
0265   if (dynamic_cast<const StripGeomDetUnit*>(it) == nullptr && dynamic_cast<const PixelGeomDetUnit*>(it) == nullptr) {
0266     std::cout << "this detID doesn't seem to belong to the Tracker" << std::endl;
0267   } else {
0268     module_width = it->surface().bounds().width();
0269   }
0270   return module_width;
0271 }
0272 
0273 //---------------------------------------------------------------------------------------------------------
0274 double SiStripGainCosmicCalculator::moduleThickness(const uint32_t detid)  // get thickness of the module detid
0275 {                                                                          //dk: copied from A. Giammanco and hacked
0276   double module_thickness = 0.;
0277   const GeomDetUnit* it = tkGeom_->idToDetUnit(DetId(detid));
0278   if (dynamic_cast<const StripGeomDetUnit*>(it) == nullptr && dynamic_cast<const PixelGeomDetUnit*>(it) == nullptr) {
0279     std::cout << "this detID doesn't seem to belong to the Tracker" << std::endl;
0280   } else {
0281     module_thickness = it->surface().bounds().thickness();
0282   }
0283   return module_thickness;
0284 }
0285 
0286 //---------------------------------------------------------------------------------------------------------
0287 std::unique_ptr<SiStripApvGain> SiStripGainCosmicCalculator::getNewObject() {
0288   std::cout << "SiStripGainCosmicCalculator::getNewObject called" << std::endl;
0289 
0290   std::cout << "total_nr_of_events=" << total_nr_of_events << std::endl;
0291   // book some more histograms
0292   TH1F* ChargeOfEachAPVPair = new TH1F("ChargeOfEachAPVPair", "ChargeOfEachAPVPair", 1, 0, 1);
0293   ChargeOfEachAPVPair->SetCanExtend(TH1::kAllAxes);
0294   TH1F* EntriesApvPairs = new TH1F("EntriesApvPairs", "EntriesApvPairs", 1, 0, 1);
0295   EntriesApvPairs->SetCanExtend(TH1::kAllAxes);
0296   TH1F* NrOfEntries =
0297       new TH1F("NrOfEntries", "NrOfEntries", 351, -0.5, 350.5);  // NrOfEntries->SetCanExtend(TH1::kAllAxes);
0298   TH1F* ModuleThickness = new TH1F("ModuleThickness", "ModuleThickness", 2, 0.5, 2.5);
0299   HlistOtherHistos->Add(ModuleThickness);
0300   ModuleThickness->GetXaxis()->SetBinLabel(1, "320mu");
0301   ModuleThickness->GetXaxis()->SetBinLabel(2, "500mu");
0302   ModuleThickness->SetYTitle("Nr APVPairs");
0303   TH1F* ModuleWidth = new TH1F("ModuleWidth", "ModuleWidth", 5, 0.5, 5.5);
0304   HlistOtherHistos->Add(ModuleWidth);
0305   ModuleWidth->GetXaxis()->SetBinLabel(1, "6.144cm");
0306   ModuleWidth->GetXaxis()->SetBinLabel(2, "7.14cm");
0307   ModuleWidth->GetXaxis()->SetBinLabel(3, "9.3696cm");
0308   ModuleWidth->GetXaxis()->SetBinLabel(4, "10.49cm");
0309   ModuleWidth->GetXaxis()->SetBinLabel(5, "12.03cm");
0310   ModuleWidth->SetYTitle("Nr APVPairs");
0311   // loop over single histograms and extract peak value of charge
0312   HlistAPVPairs->Sort();  // sort alfabetically
0313   TIter hiterator(HlistAPVPairs);
0314   double MeanCharge = 0.;
0315   double NrOfApvPairs = 0.;
0316   TH1F* MyHisto = (TH1F*)hiterator();
0317   while (MyHisto) {
0318     TString histo_title = MyHisto->GetTitle();
0319     if (histo_title.Contains("ChargeAPVPair_")) {
0320       std::pair<double, double> two_values = getPeakOfLandau(MyHisto);
0321       double local_nrofadcs = two_values.first;
0322       double local_sigma = two_values.second;
0323       ChargeOfEachAPVPair->Fill(histo_title, local_nrofadcs);
0324       int ichbin = ChargeOfEachAPVPair->GetXaxis()->FindBin(histo_title.Data());
0325       ChargeOfEachAPVPair->SetBinError(ichbin, local_sigma);
0326       EntriesApvPairs->Fill(histo_title, MyHisto->GetEntries());
0327       NrOfEntries->Fill(MyHisto->GetEntries());
0328       if (local_nrofadcs > 0) {  // if nr of adcs is negative, the fitting routine could not extract meaningfull numbers
0329         MeanCharge += local_nrofadcs;
0330         NrOfApvPairs += 1.;  // count nr of apv pairs since do not know whether nr of bins of histogram is the same
0331       }
0332     }
0333     MyHisto = (TH1F*)hiterator();
0334   }
0335   ChargeOfEachAPVPair->LabelsDeflate("X");
0336   EntriesApvPairs->LabelsDeflate("X");  // trim nr. of bins to match active labels
0337   HlistOtherHistos->Add(ChargeOfEachAPVPair);
0338   HlistOtherHistos->Add(EntriesApvPairs);
0339   HlistOtherHistos->Add(NrOfEntries);
0340   MeanCharge = MeanCharge / NrOfApvPairs;
0341   // calculate correction
0342   TH1F* CorrectionOfEachAPVPair = (TH1F*)ChargeOfEachAPVPair->Clone("CorrectionOfEachAPVPair");
0343   TH1F* ChargeOfEachAPVPairControlView =
0344       new TH1F("ChargeOfEachAPVPairControlView", "ChargeOfEachAPVPairControlView", 1, 0, 1);
0345   ChargeOfEachAPVPairControlView->SetCanExtend(TH1::kAllAxes);
0346   TH1F* CorrectionOfEachAPVPairControlView =
0347       new TH1F("CorrectionOfEachAPVPairControlView", "CorrectionOfEachAPVPairControlView", 1, 0, 1);
0348   CorrectionOfEachAPVPairControlView->SetCanExtend(TH1::kAllAxes);
0349   std::ofstream APVPairTextOutput("apvpair_corrections.txt");
0350   APVPairTextOutput << "# MeanCharge = " << MeanCharge << std::endl;
0351   APVPairTextOutput << "# Nr. of APVPairs = " << NrOfApvPairs << std::endl;
0352   for (int ibin = 1; ibin <= ChargeOfEachAPVPair->GetNbinsX(); ibin++) {
0353     TString local_bin_label = ChargeOfEachAPVPair->GetXaxis()->GetBinLabel(ibin);
0354     double local_charge_over_path = ChargeOfEachAPVPair->GetBinContent(ibin);
0355     if (local_bin_label.Contains("ChargeAPVPair_") &&
0356         local_charge_over_path > 0.0000001) {  // calculate correction only for meaningful numbers
0357       uint32_t extracted_detid;
0358       std::istringstream read_label((local_bin_label(14, 9)).Data());
0359       read_label >> extracted_detid;
0360       unsigned short extracted_apvpairid;
0361       std::istringstream read_apvpair((local_bin_label(24, 1)).Data());
0362       read_apvpair >> extracted_apvpairid;
0363       double local_error_of_charge = ChargeOfEachAPVPair->GetBinError(ibin);
0364       double local_correction = -0.5;
0365       double local_error_correction = 0.;
0366       local_correction =
0367           MeanCharge / local_charge_over_path;  // later use ExpectedChargeDeposition instead of MeanCharge
0368       local_error_correction = local_correction * local_error_of_charge / local_charge_over_path;
0369       if (local_error_correction > 1.8) {  // understand why error too large sometimes
0370         std::cout << "too large error " << local_error_correction << " for histogram " << local_bin_label << std::endl;
0371       }
0372       double nr_of_entries = EntriesApvPairs->GetBinContent(ibin);
0373       APVPairTextOutput << local_bin_label << " " << local_correction << " " << local_charge_over_path << " "
0374                         << nr_of_entries << std::endl;
0375       CorrectionOfEachAPVPair->SetBinContent(ibin, local_correction);
0376       CorrectionOfEachAPVPair->SetBinError(ibin, local_error_correction);
0377       ((TH1F*)HlistOtherHistos->FindObject("APVPairCorrections"))->Fill(local_correction);
0378       DetId thedetId = DetId(extracted_detid);
0379       unsigned int generalized_layer = 0;
0380       // calculate generalized_layer:  31,32 = TIB1, 33 = TIB2, 33 = TIB3, 51 = TOB1, 52 = TOB2, 60 = TEC
0381       if (thedetId.subdetId() == StripSubdetector::TIB) {
0382         generalized_layer =
0383             10 * thedetId.subdetId() + tTopo_->tibLayer(thedetId.rawId()) + tTopo_->tibStereo(thedetId.rawId());
0384         if (tTopo_->tibLayer(thedetId.rawId()) == 2) {
0385           generalized_layer++;
0386           if (tTopo_->tibGlued(thedetId.rawId()))
0387             edm::LogError("ClusterMTCCFilter") << "WRONGGGG" << std::endl;
0388         }
0389       } else {
0390         generalized_layer = 10 * thedetId.subdetId();
0391         if (thedetId.subdetId() == StripSubdetector::TOB) {
0392           generalized_layer += tTopo_->tobLayer(thedetId.rawId());
0393         }
0394       }
0395       if (generalized_layer == 31) {
0396         ((TH1F*)HlistOtherHistos->FindObject("APVPairCorrectionsTIB1mono"))->Fill(local_correction);
0397       }
0398       if (generalized_layer == 32) {
0399         ((TH1F*)HlistOtherHistos->FindObject("APVPairCorrectionsTIB1stereo"))->Fill(local_correction);
0400       }
0401       if (generalized_layer == 33) {
0402         ((TH1F*)HlistOtherHistos->FindObject("APVPairCorrectionsTIB2"))->Fill(local_correction);
0403       }
0404       if (generalized_layer == 51) {
0405         ((TH1F*)HlistOtherHistos->FindObject("APVPairCorrectionsTOB1"))->Fill(local_correction);
0406       }
0407       if (generalized_layer == 52) {
0408         ((TH1F*)HlistOtherHistos->FindObject("APVPairCorrectionsTOB2"))->Fill(local_correction);
0409       }
0410       // control view
0411       const FedChannelConnection& fedchannelconnection =
0412           siStripDetCabling_->getConnection(extracted_detid, extracted_apvpairid);
0413       std::ostringstream local_key;
0414       // in S. Mersi's analysis the APVPair id seems to be used instead of the lldChannel, hence use the same here
0415       local_key << "fecCrate" << fedchannelconnection.fecCrate() << "_fecSlot" << fedchannelconnection.fecSlot()
0416                 << "_fecRing" << fedchannelconnection.fecRing() << "_ccuAddr" << fedchannelconnection.ccuAddr()
0417                 << "_ccuChan" << fedchannelconnection.ccuChan() << "_apvPair" << extracted_apvpairid;
0418       TString control_key = local_key.str();
0419       ChargeOfEachAPVPairControlView->Fill(control_key, local_charge_over_path);
0420       int ibin1 = ChargeOfEachAPVPairControlView->GetXaxis()->FindBin(control_key);
0421       ChargeOfEachAPVPairControlView->SetBinError(ibin1, local_error_of_charge);
0422       CorrectionOfEachAPVPairControlView->Fill(control_key, local_correction);
0423       int ibin2 = CorrectionOfEachAPVPairControlView->GetXaxis()->FindBin(control_key);
0424       CorrectionOfEachAPVPairControlView->SetBinError(ibin2, local_error_correction);
0425       // thickness of each module
0426       double module_thickness = moduleThickness(extracted_detid);
0427       if (fabs(module_thickness - 0.032) < 0.001)
0428         ModuleThickness->Fill(1);
0429       if (fabs(module_thickness - 0.05) < 0.001)
0430         ModuleThickness->Fill(2);
0431       // width of each module
0432       double module_width = moduleWidth(extracted_detid);
0433       if (fabs(module_width - 6.144) < 0.01)
0434         ModuleWidth->Fill(1);
0435       if (fabs(module_width - 7.14) < 0.01)
0436         ModuleWidth->Fill(2);
0437       if (fabs(module_width - 9.3696) < 0.01)
0438         ModuleWidth->Fill(3);
0439       if (fabs(module_width - 10.49) < 0.01)
0440         ModuleWidth->Fill(4);
0441       if (fabs(module_width - 12.03) < 0.01)
0442         ModuleWidth->Fill(5);
0443     }
0444   }
0445   HlistOtherHistos->Add(CorrectionOfEachAPVPair);
0446   ChargeOfEachAPVPairControlView->LabelsDeflate("X");
0447   CorrectionOfEachAPVPairControlView->LabelsDeflate("X");
0448   HlistOtherHistos->Add(ChargeOfEachAPVPairControlView);
0449   HlistOtherHistos->Add(CorrectionOfEachAPVPairControlView);
0450   // output histograms to file
0451 
0452   if (outputHistogramsInRootFile) {
0453     TFile* outputfile = new TFile(outputFileName, "RECREATE");
0454     HlistAPVPairs->Write();
0455     HlistOtherHistos->Write();
0456     outputfile->Close();
0457   }
0458 
0459   auto obj = std::make_unique<SiStripApvGain>();
0460 
0461   //   for(std::map<uint32_t,OptoScanAnalysis*>::const_iterator it = analyses.begin(); it != analyses.end(); it++){
0462   //     //Generate Gain for det detid
0463   //     std::vector<float> theSiStripVector;
0464   //     for(unsigned short j=0; j<it->second; j++){
0465   //       float gain;
0466 
0467   //       //      if(sigmaGain_/meanGain_ < 0.00001) gain = meanGain_;
0468   //       //      else{
0469   //       gain = CLHEP::RandGauss::shoot(meanGain_, sigmaGain_);
0470   //       if(gain<=minimumPosValue_) gain=minimumPosValue_;
0471   //       //      }
0472 
0473   //       if (printdebug_)
0474   //    edm::LogInfo("SiStripGainCalculator") << "detid " << it->first << " \t"
0475   //                          << " apv " << j << " \t"
0476   //                          << gain    << " \t"
0477   //                          << std::endl;
0478   //       theSiStripVector.push_back(gain);
0479   //     }
0480   //     SiStripApvGain::Range range(theSiStripVector.begin(),theSiStripVector.end());
0481   //     if ( ! obj->put(it->first,range) )
0482   //       edm::LogError("SiStripGainCalculator")<<"[SiStripGainCalculator::beginJob] detid already exists"<<std::endl;
0483   //   }
0484 
0485   return obj;
0486 }