#include <SiStripGainCosmicCalculator.h>

Inheritance diagram for SiStripGainCosmicCalculator:

Public Member Functions
	SiStripGainCosmicCalculator (const edm::ParameterSet &)

	~SiStripGainCosmicCalculator () override

Public Member Functions inherited from ConditionDBWriter< SiStripApvGain >
	ConditionDBWriter (const edm::ParameterSet &iConfig)

	~ConditionDBWriter () override

Public Member Functions inherited from edm::one::EDAnalyzer< edm::one::WatchRuns, edm::one::WatchLuminosityBlocks, edm::one::SharedResources >
	EDAnalyzer ()=default

	EDAnalyzer (const EDAnalyzer &)=delete

SerialTaskQueue *	globalLuminosityBlocksQueue () final

SerialTaskQueue *	globalRunsQueue () final

const EDAnalyzer &	operator= (const EDAnalyzer &)=delete

bool	wantsGlobalLuminosityBlocks () const noexcept final

bool	wantsGlobalRuns () const noexcept final

bool	wantsInputProcessBlocks () const noexcept final

bool	wantsProcessBlocks () const noexcept final

Public Member Functions inherited from edm::one::EDAnalyzerBase
void	callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)

	EDAnalyzerBase ()

ModuleDescription const &	moduleDescription () const

bool	wantsStreamLuminosityBlocks () const noexcept

bool	wantsStreamRuns () const noexcept

	~EDAnalyzerBase () override

Public Member Functions inherited from edm::EDConsumerBase
std::vector< ConsumesInfo >	consumesInfo () const

void	convertCurrentProcessAlias (std::string const &processName)
	Convert "@currentProcess" in InputTag process names to the actual current process name. More...

	EDConsumerBase ()

	EDConsumerBase (EDConsumerBase const &)=delete

	EDConsumerBase (EDConsumerBase &&)=default

ESResolverIndex const *	esGetTokenIndices (edm::Transition iTrans) const

std::vector< ESResolverIndex > const &	esGetTokenIndicesVector (edm::Transition iTrans) const

std::vector< ESRecordIndex > const &	esGetTokenRecordIndicesVector (edm::Transition iTrans) const

ProductResolverIndexAndSkipBit	indexFrom (EDGetToken, BranchType, TypeID const &) const

void	itemsMayGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const

void	itemsToGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const

std::vector< ProductResolverIndexAndSkipBit > const &	itemsToGetFrom (BranchType iType) const

void	labelsForToken (EDGetToken iToken, Labels &oLabels) const

void	modulesWhoseProductsAreConsumed (std::array< std::vector< ModuleDescription const > , NumBranchTypes > &modulesAll, std::vector< ModuleProcessName > &modulesInPreviousProcesses, ProductRegistry const &preg, std::map< std::string, ModuleDescription const *> const &labelsToDesc, std::string const &processName) const

EDConsumerBase const &	operator= (EDConsumerBase const &)=delete

EDConsumerBase &	operator= (EDConsumerBase &&)=default

bool	registeredToConsume (ProductResolverIndex, bool, BranchType) const

void	selectInputProcessBlocks (ProductRegistry const &productRegistry, ProcessBlockHelperBase const &processBlockHelperBase)

ProductResolverIndexAndSkipBit	uncheckedIndexFrom (EDGetToken) const

void	updateLookup (BranchType iBranchType, ProductResolverIndexHelper const &, bool iPrefetchMayGet)

void	updateLookup (eventsetup::ESRecordsToProductResolverIndices const &)

virtual	~EDConsumerBase () noexcept(false)

Private Member Functions
void	algoAnalyze (const edm::Event &, const edm::EventSetup &) override

void	algoBeginJob (const edm::EventSetup &) override

void	algoEndJob () override

std::unique_ptr< SiStripApvGain >	getNewObject () override

std::pair< double, double >	getPeakOfLandau (TH1F *inputHisto)

double	moduleThickness (const uint32_t detid)

double	moduleWidth (const uint32_t detid)

Private Attributes
edm::ESGetToken< SiStripDetCabling, SiStripDetCablingRcd >	detCablingToken_

std::vector< uint32_t >	detModulesToBeExcluded

double	ExpectedChargeDeposition

TObjArray *	HlistAPVPairs

TObjArray *	HlistOtherHistos

double	MaxChi2OverNDF

unsigned int	MinNrEntries

TString	outputFileName

bool	outputHistogramsInRootFile

bool	printdebug_

std::vector< uint32_t >	SelectedDetIds

const SiStripDetCabling *	siStripDetCabling_ = nullptr

std::map< uint32_t, double >	thickness_map

const TrackerGeometry *	tkGeom_ = nullptr

edm::ESGetToken< TrackerGeometry, TrackerDigiGeometryRecord >	tkGeomToken_

uint32_t	total_nr_of_events

std::string	TrackLabel

std::string	TrackProducer

const TrackerTopology *	tTopo_ = nullptr

edm::ESGetToken< TrackerTopology, TrackerTopologyRcd >	tTopoToken_

Additional Inherited Members
Public Types inherited from edm::one::EDAnalyzerBase
typedef EDAnalyzerBase	ModuleType

Public Types inherited from edm::EDConsumerBase
typedef ProductLabels	Labels

Static Public Member Functions inherited from ConditionDBWriter< SiStripApvGain >
static void	fillPSetDescription (edm::ParameterSetDescription &desc)

Static Public Member Functions inherited from edm::one::EDAnalyzerBase
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &descriptions)

Protected Member Functions inherited from ConditionDBWriter< SiStripApvGain >
void	setDoStore (const bool doStore)
	When set to false the payload will not be written to the db. More...

void	storeOnDbNow ()

cond::Time_t	timeOfLastIOV ()

Protected Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	consumes (edm::InputTag const &tag)

template<BranchType B = InEvent>
EDConsumerBaseAdaptor< B >	consumes (edm::InputTag tag) noexcept

EDGetToken	consumes (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	consumes (TypeToGet const &id, edm::InputTag const &tag)

ConsumesCollector	consumesCollector ()
	Use a ConsumesCollector to gather consumes information from helper functions. More...

template<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto	esConsumes ()

template<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto	esConsumes (ESInputTag const &tag)

template<Transition Tr = Transition::Event>
constexpr auto	esConsumes ()

template<Transition Tr = Transition::Event>
auto	esConsumes (ESInputTag tag)

template<Transition Tr = Transition::Event>
ESGetTokenGeneric	esConsumes (eventsetup::EventSetupRecordKey const &iRecord, eventsetup::DataKey const &iKey)
	Used with EventSetupRecord::doGet. More...

template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	mayConsume (edm::InputTag const &tag)

EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

void	resetItemsToGetFrom (BranchType iType)

Detailed Description

Definition at line 32 of file SiStripGainCosmicCalculator.h.

Constructor & Destructor Documentation

◆ SiStripGainCosmicCalculator()

SiStripGainCosmicCalculator::SiStripGainCosmicCalculator ( const edm::ParameterSet & iConfig )

explicit

Definition at line 30 of file SiStripGainCosmicCalculator.cc.

References detCablingToken_, detModulesToBeExcluded, ExpectedChargeDeposition, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), MaxChi2OverNDF, MinNrEntries, outputFileName, outputHistogramsInRootFile, printdebug_, AlCaHLTBitMon_QueryRunRegistry::string, tkGeomToken_, TrackLabel, and tTopoToken_.

     : ConditionDBWriter<SiStripApvGain>(iConfig) {
   edm::LogInfo("SiStripGainCosmicCalculator::SiStripGainCosmicCalculator");
   ExpectedChargeDeposition = 200.;
   edm::LogInfo("SiStripApvGainCalculator::SiStripApvGainCalculator")
       << "ExpectedChargeDeposition=" << ExpectedChargeDeposition;
 
   TrackProducer = iConfig.getParameter<std::string>("TrackProducer");
   TrackLabel = iConfig.getParameter<std::string>("TrackLabel");
 
   detModulesToBeExcluded.clear();
   detModulesToBeExcluded = iConfig.getParameter<std::vector<unsigned> >("detModulesToBeExcluded");
   MinNrEntries = iConfig.getUntrackedParameter<unsigned>("minNrEntries", 20);
   MaxChi2OverNDF = iConfig.getUntrackedParameter<double>("maxChi2OverNDF", 5.);
 
   outputHistogramsInRootFile = iConfig.getParameter<bool>("OutputHistogramsInRootFile");
   outputFileName = iConfig.getParameter<std::string>("OutputFileName");
 
   edm::LogInfo("SiStripApvGainCalculator")
       << "Clusters from " << detModulesToBeExcluded.size() << " modules will be ignored in the calibration:";
   edm::LogInfo("SiStripApvGainCalculator") << "The calibration for these DetIds will be set to a default value";
   for (std::vector<uint32_t>::const_iterator imod = detModulesToBeExcluded.begin();
        imod != detModulesToBeExcluded.end();
        imod++) {
     edm::LogInfo("SiStripApvGainCalculator") << "exclude detid = " << *imod;
   }
 
   printdebug_ = iConfig.getUntrackedParameter<bool>("printDebug", false);
 
   tTopoToken_ = esConsumes<edm::Transition::BeginRun>();
   detCablingToken_ = esConsumes<edm::Transition::BeginRun>();
   tkGeomToken_ = esConsumes<edm::Transition::BeginRun>();
 }

◆ ~SiStripGainCosmicCalculator()

SiStripGainCosmicCalculator::~SiStripGainCosmicCalculator ( )

override

Definition at line 64 of file SiStripGainCosmicCalculator.cc.

                                                           {
   edm::LogInfo("SiStripGainCosmicCalculator::~SiStripGainCosmicCalculator");
 }

Member Function Documentation

◆ algoAnalyze()

void SiStripGainCosmicCalculator::algoAnalyze	(	const edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)

overrideprivatevirtual

Reimplemented from ConditionDBWriter< SiStripApvGain >.

Definition at line 145 of file SiStripGainCosmicCalculator.cc.

References SiStripRecHit2D::cluster(), funct::cos(), HcalObjRepresent::Fill(), HlistAPVPairs, HlistOtherHistos, iEvent, createfilelist::int, TrackingRecHit::localPosition(), moduleThickness(), moduleWidth(), total_nr_of_events, JetHT_cfg::trackCollection, TrackLabel, DiMuonV_cfg::tracks, trackerHitRTTI::vector, and PV3DBase< T, PVType, FrameType >::x().

                                                                                                  {
   using namespace edm;
   total_nr_of_events++;
 
   //TO BE RESTORED
   //  anglefinder_->init(event,iSetup);
 
   // get seeds
   //  edm::Handle<TrajectorySeedCollection> seedcoll;
   //  event.getByType(seedcoll);
   // get tracks
   Handle<reco::TrackCollection> trackCollection;
   iEvent.getByLabel(TrackProducer, TrackLabel, trackCollection);
   const reco::TrackCollection* tracks = trackCollection.product();
 
   //  // get magnetic field
   //  edm::ESHandle<MagneticField> esmagfield;
   //  es.get<IdealMagneticFieldRecord>().get(esmagfield);
   //  magfield=&(*esmagfield);
   // loop over tracks
   for (reco::TrackCollection::const_iterator itr = tracks->begin(); itr != tracks->end();
        itr++) {  // looping over tracks
 
     //TO BE RESTORED
     //    std::vector<std::pair<const TrackingRecHit *,float> >hitangle =anglefinder_->findtrackangle((*(*seedcoll).begin()),*itr);
     std::vector<std::pair<const TrackingRecHit*, float> >
         hitangle;  // =anglefinder_->findtrackangle((*(*seedcoll).begin()),*itr);
 
     for (std::vector<std::pair<const TrackingRecHit*, float> >::const_iterator hitangle_iter = hitangle.begin();
          hitangle_iter != hitangle.end();
          hitangle_iter++) {
       const TrackingRecHit* trechit = hitangle_iter->first;
       float local_angle = hitangle_iter->second;
       LocalPoint local_position = trechit->localPosition();
       const SiStripRecHit2D* sistripsimplehit = dynamic_cast<const SiStripRecHit2D*>(trechit);
       const SiStripMatchedRecHit2D* sistripmatchedhit = dynamic_cast<const SiStripMatchedRecHit2D*>(trechit);
       //      std::cout<<" hit/matched "<<std::ios::hex<<sistripsimplehit<<" "<<sistripmatchedhit<<std::endl;
       ((TH1F*)HlistOtherHistos->FindObject("LocalAngle"))->Fill(local_angle);
       ((TH1F*)HlistOtherHistos->FindObject("LocalAngleAbsoluteCosine"))->Fill(fabs(cos(local_angle)));
       if (sistripsimplehit) {
         ((TH1F*)HlistOtherHistos->FindObject("SiStripRecHitType"))->Fill(1.);
         const SiStripRecHit2D::ClusterRef& cluster = sistripsimplehit->cluster();
         const auto& ampls = cluster->amplitudes();
         uint32_t thedetid = 0;  // is zero since long time cluster->geographicalId();
         double module_width = moduleWidth(thedetid);
         ((TH1F*)HlistOtherHistos->FindObject("LocalPosition_cm"))->Fill(local_position.x());
         ((TH1F*)HlistOtherHistos->FindObject("LocalPosition_normalized"))->Fill(local_position.x() / module_width);
         double module_thickness = moduleThickness(thedetid);
         int ifirststrip = cluster->firstStrip();
         int theapvpairid = int(float(ifirststrip) / 256.);
         TH1F* histopointer = (TH1F*)HlistAPVPairs->FindObject(Form("ChargeAPVPair_%i_%i", thedetid, theapvpairid));
         if (histopointer) {
           short cCharge = 0;
           for (unsigned int iampl = 0; iampl < ampls.size(); iampl++) {
             cCharge += ampls[iampl];
           }
           double cluster_charge_over_path = ((double)cCharge) * fabs(cos(local_angle)) / (10. * module_thickness);
           histopointer->Fill(cluster_charge_over_path);
         }
       } else {
         if (sistripmatchedhit)
           ((TH1F*)HlistOtherHistos->FindObject("SiStripRecHitType"))->Fill(2.);
       }
     }
   }
 }

◆ algoBeginJob()

void SiStripGainCosmicCalculator::algoBeginJob ( const edm::EventSetup & iSetup )

overrideprivatevirtual

Reimplemented from ConditionDBWriter< SiStripApvGain >.

Definition at line 70 of file SiStripGainCosmicCalculator.cc.

References SiStripDetCabling::addActiveDetectorsRawIds(), gather_cfg::cout, detCablingToken_, detModulesToBeExcluded, TrackerGeometry::dets(), beamvalidation::exit(), spr::find(), edm::EventSetup::getData(), SiStripSubStructure::getTIBDetectors(), SiStripSubStructure::getTOBDetectors(), HlistAPVPairs, HlistOtherHistos, AlCaHLTBitMon_ParallelJobs::p, SelectedDetIds, siStripDetCabling_, thickness_map, tkGeom_, tkGeomToken_, total_nr_of_events, tTopo_, and tTopoToken_.

                                                                           {
   tTopo_ = &iSetup.getData(tTopoToken_);
   siStripDetCabling_ = &iSetup.getData(detCablingToken_);
   tkGeom_ = &iSetup.getData(tkGeomToken_);
 
   std::cout << "SiStripGainCosmicCalculator::algoBeginJob called" << std::endl;
   total_nr_of_events = 0;
   HlistAPVPairs = new TObjArray();
   HlistOtherHistos = new TObjArray();
   //
   HlistOtherHistos->Add(new TH1F(Form("APVPairCorrections"), Form("APVPairCorrections"), 50, -1., 4.));
   HlistOtherHistos->Add(new TH1F(Form("APVPairCorrectionsTIB1mono"), Form("APVPairCorrectionsTIB1mono"), 50, -1., 4.));
   HlistOtherHistos->Add(
       new TH1F(Form("APVPairCorrectionsTIB1stereo"), Form("APVPairCorrectionsTIB1stereo"), 50, -1., 4.));
   HlistOtherHistos->Add(new TH1F(Form("APVPairCorrectionsTIB2"), Form("APVPairCorrectionsTIB2"), 50, -1., 4.));
   HlistOtherHistos->Add(new TH1F(Form("APVPairCorrectionsTOB1"), Form("APVPairCorrectionsTOB1"), 50, -1., 4.));
   HlistOtherHistos->Add(new TH1F(Form("APVPairCorrectionsTOB2"), Form("APVPairCorrectionsTOB2"), 50, -1., 4.));
   HlistOtherHistos->Add(new TH1F(Form("LocalAngle"), Form("LocalAngle"), 70, -0.1, 3.4));
   HlistOtherHistos->Add(new TH1F(Form("LocalAngleAbsoluteCosine"), Form("LocalAngleAbsoluteCosine"), 48, -0.1, 1.1));
   HlistOtherHistos->Add(new TH1F(Form("LocalPosition_cm"), Form("LocalPosition_cm"), 100, -5., 5.));
   HlistOtherHistos->Add(new TH1F(Form("LocalPosition_normalized"), Form("LocalPosition_normalized"), 100, -1.1, 1.1));
   TH1F* local_histo = new TH1F(Form("SiStripRecHitType"), Form("SiStripRecHitType"), 2, 0.5, 2.5);
   HlistOtherHistos->Add(local_histo);
   local_histo->GetXaxis()->SetBinLabel(1, "simple");
   local_histo->GetXaxis()->SetBinLabel(2, "matched");
 
   // get cabling and find out list of active detectors
   std::vector<uint32_t> activeDets;
   activeDets.clear();
   SelectedDetIds.clear();
   siStripDetCabling_->addActiveDetectorsRawIds(activeDets);
   //    SelectedDetIds = activeDets; // all active detector modules
   // use SiStripSubStructure for selecting certain regions
   SiStripSubStructure::getTIBDetectors(
       activeDets, SelectedDetIds, tTopo_, 0, 0, 0, 0);  // this adds rawDetIds to SelectedDetIds
   SiStripSubStructure::getTOBDetectors(
       activeDets, SelectedDetIds, tTopo_, 0, 0, 0);  // this adds rawDetIds to SelectedDetIds
   // get tracker geometry and find nr. of apv pairs for each active detector
   for (auto det : tkGeom_->dets()) {  // loop over detector modules
     if (dynamic_cast<const StripGeomDetUnit*>(det) != nullptr) {
       uint32_t detid = det->geographicalId().rawId();
       // get thickness for all detector modules, not just for active, this is strange
       double module_thickness =
           det->surface()
               .bounds()
               .thickness();  // get thickness of detector from GeomDet (DetContainer == vector<GeomDet*>)
       thickness_map.insert(std::make_pair(detid, module_thickness));
       //
       const bool is_active_detector =
           std::end(SelectedDetIds) != std::find(std::begin(SelectedDetIds), std::end(SelectedDetIds), detid);
       //
       const bool exclude_this_detid =
           std::end(detModulesToBeExcluded) !=
           std::find(std::begin(detModulesToBeExcluded), std::end(detModulesToBeExcluded), detid);
       //
       if (is_active_detector &&
           (!exclude_this_detid)) {  // check whether is active detector and that should not be excluded
         const StripTopology& p = dynamic_cast<const StripGeomDetUnit*>(det)->specificTopology();
         unsigned short NAPVPairs = p.nstrips() / 256;
         if (NAPVPairs < 2 || NAPVPairs > 3) {
           edm::LogError("SiStripGainCosmicCalculator")
               << "Problem with Number of strips in detector: " << p.nstrips() << " Exiting program";
           exit(1);
         }
         for (int iapp = 0; iapp < NAPVPairs; iapp++) {
           TString hid = Form("ChargeAPVPair_%i_%i", detid, iapp);
           HlistAPVPairs->Add(
               new TH1F(hid, hid, 45, 0., 1350.));  // multiply by 3 to take into account division by width
         }
       }
     }
   }
 }

◆ algoEndJob()

void SiStripGainCosmicCalculator::algoEndJob ( )

overrideprivatevirtual

Reimplemented from ConditionDBWriter< SiStripApvGain >.

Definition at line 68 of file SiStripGainCosmicCalculator.cc.

68 {}

◆ getNewObject()

std::unique_ptr< SiStripApvGain > SiStripGainCosmicCalculator::getNewObject ( )

overrideprivatevirtual

Implements ConditionDBWriter< SiStripApvGain >.

Definition at line 287 of file SiStripGainCosmicCalculator.cc.

References FedChannelConnection::ccuAddr(), FedChannelConnection::ccuChan(), gather_cfg::cout, FedChannelConnection::fecCrate(), FedChannelConnection::fecRing(), FedChannelConnection::fecSlot(), HcalObjRepresent::Fill(), SiStripDetCabling::getConnection(), getPeakOfLandau(), HlistAPVPairs, HlistOtherHistos, moduleThickness(), moduleWidth(), getGTfromDQMFile::obj, writedatasetfile::outputfile, outputFileName, outputHistogramsInRootFile, DetId::rawId(), siStripDetCabling_, DetId::subdetId(), StripSubdetector::TIB, TrackerTopology::tibGlued(), TrackerTopology::tibLayer(), TrackerTopology::tibStereo(), StripSubdetector::TOB, TrackerTopology::tobLayer(), total_nr_of_events, and tTopo_.

                                                                         {
   std::cout << "SiStripGainCosmicCalculator::getNewObject called" << std::endl;
 
   std::cout << "total_nr_of_events=" << total_nr_of_events << std::endl;
   // book some more histograms
   TH1F* ChargeOfEachAPVPair = new TH1F("ChargeOfEachAPVPair", "ChargeOfEachAPVPair", 1, 0, 1);
   ChargeOfEachAPVPair->SetCanExtend(TH1::kAllAxes);
   TH1F* EntriesApvPairs = new TH1F("EntriesApvPairs", "EntriesApvPairs", 1, 0, 1);
   EntriesApvPairs->SetCanExtend(TH1::kAllAxes);
   TH1F* NrOfEntries =
       new TH1F("NrOfEntries", "NrOfEntries", 351, -0.5, 350.5);  // NrOfEntries->SetCanExtend(TH1::kAllAxes);
   TH1F* ModuleThickness = new TH1F("ModuleThickness", "ModuleThickness", 2, 0.5, 2.5);
   HlistOtherHistos->Add(ModuleThickness);
   ModuleThickness->GetXaxis()->SetBinLabel(1, "320mu");
   ModuleThickness->GetXaxis()->SetBinLabel(2, "500mu");
   ModuleThickness->SetYTitle("Nr APVPairs");
   TH1F* ModuleWidth = new TH1F("ModuleWidth", "ModuleWidth", 5, 0.5, 5.5);
   HlistOtherHistos->Add(ModuleWidth);
   ModuleWidth->GetXaxis()->SetBinLabel(1, "6.144cm");
   ModuleWidth->GetXaxis()->SetBinLabel(2, "7.14cm");
   ModuleWidth->GetXaxis()->SetBinLabel(3, "9.3696cm");
   ModuleWidth->GetXaxis()->SetBinLabel(4, "10.49cm");
   ModuleWidth->GetXaxis()->SetBinLabel(5, "12.03cm");
   ModuleWidth->SetYTitle("Nr APVPairs");
   // loop over single histograms and extract peak value of charge
   HlistAPVPairs->Sort();  // sort alfabetically
   TIter hiterator(HlistAPVPairs);
   double MeanCharge = 0.;
   double NrOfApvPairs = 0.;
   TH1F* MyHisto = (TH1F*)hiterator();
   while (MyHisto) {
     TString histo_title = MyHisto->GetTitle();
     if (histo_title.Contains("ChargeAPVPair_")) {
       std::pair<double, double> two_values = getPeakOfLandau(MyHisto);
       double local_nrofadcs = two_values.first;
       double local_sigma = two_values.second;
       ChargeOfEachAPVPair->Fill(histo_title, local_nrofadcs);
       int ichbin = ChargeOfEachAPVPair->GetXaxis()->FindBin(histo_title.Data());
       ChargeOfEachAPVPair->SetBinError(ichbin, local_sigma);
       EntriesApvPairs->Fill(histo_title, MyHisto->GetEntries());
       NrOfEntries->Fill(MyHisto->GetEntries());
       if (local_nrofadcs > 0) {  // if nr of adcs is negative, the fitting routine could not extract meaningfull numbers
         MeanCharge += local_nrofadcs;
         NrOfApvPairs += 1.;  // count nr of apv pairs since do not know whether nr of bins of histogram is the same
       }
     }
     MyHisto = (TH1F*)hiterator();
   }
   ChargeOfEachAPVPair->LabelsDeflate("X");
   EntriesApvPairs->LabelsDeflate("X");  // trim nr. of bins to match active labels
   HlistOtherHistos->Add(ChargeOfEachAPVPair);
   HlistOtherHistos->Add(EntriesApvPairs);
   HlistOtherHistos->Add(NrOfEntries);
   MeanCharge = MeanCharge / NrOfApvPairs;
   // calculate correction
   TH1F* CorrectionOfEachAPVPair = (TH1F*)ChargeOfEachAPVPair->Clone("CorrectionOfEachAPVPair");
   TH1F* ChargeOfEachAPVPairControlView =
       new TH1F("ChargeOfEachAPVPairControlView", "ChargeOfEachAPVPairControlView", 1, 0, 1);
   ChargeOfEachAPVPairControlView->SetCanExtend(TH1::kAllAxes);
   TH1F* CorrectionOfEachAPVPairControlView =
       new TH1F("CorrectionOfEachAPVPairControlView", "CorrectionOfEachAPVPairControlView", 1, 0, 1);
   CorrectionOfEachAPVPairControlView->SetCanExtend(TH1::kAllAxes);
   std::ofstream APVPairTextOutput("apvpair_corrections.txt");
   APVPairTextOutput << "# MeanCharge = " << MeanCharge << std::endl;
   APVPairTextOutput << "# Nr. of APVPairs = " << NrOfApvPairs << std::endl;
   for (int ibin = 1; ibin <= ChargeOfEachAPVPair->GetNbinsX(); ibin++) {
     TString local_bin_label = ChargeOfEachAPVPair->GetXaxis()->GetBinLabel(ibin);
     double local_charge_over_path = ChargeOfEachAPVPair->GetBinContent(ibin);
     if (local_bin_label.Contains("ChargeAPVPair_") &&
         local_charge_over_path > 0.0000001) {  // calculate correction only for meaningful numbers
       uint32_t extracted_detid;
       std::istringstream read_label((local_bin_label(14, 9)).Data());
       read_label >> extracted_detid;
       unsigned short extracted_apvpairid;
       std::istringstream read_apvpair((local_bin_label(24, 1)).Data());
       read_apvpair >> extracted_apvpairid;
       double local_error_of_charge = ChargeOfEachAPVPair->GetBinError(ibin);
       double local_correction = -0.5;
       double local_error_correction = 0.;
       local_correction =
           MeanCharge / local_charge_over_path;  // later use ExpectedChargeDeposition instead of MeanCharge
       local_error_correction = local_correction * local_error_of_charge / local_charge_over_path;
       if (local_error_correction > 1.8) {  // understand why error too large sometimes
         std::cout << "too large error " << local_error_correction << " for histogram " << local_bin_label << std::endl;
       }
       double nr_of_entries = EntriesApvPairs->GetBinContent(ibin);
       APVPairTextOutput << local_bin_label << " " << local_correction << " " << local_charge_over_path << " "
                         << nr_of_entries << std::endl;
       CorrectionOfEachAPVPair->SetBinContent(ibin, local_correction);
       CorrectionOfEachAPVPair->SetBinError(ibin, local_error_correction);
       ((TH1F*)HlistOtherHistos->FindObject("APVPairCorrections"))->Fill(local_correction);
       DetId thedetId = DetId(extracted_detid);
       unsigned int generalized_layer = 0;
       // calculate generalized_layer:  31,32 = TIB1, 33 = TIB2, 33 = TIB3, 51 = TOB1, 52 = TOB2, 60 = TEC
       if (thedetId.subdetId() == StripSubdetector::TIB) {
         generalized_layer =
             10 * thedetId.subdetId() + tTopo_->tibLayer(thedetId.rawId()) + tTopo_->tibStereo(thedetId.rawId());
         if (tTopo_->tibLayer(thedetId.rawId()) == 2) {
           generalized_layer++;
           if (tTopo_->tibGlued(thedetId.rawId()))
             edm::LogError("ClusterMTCCFilter") << "WRONGGGG" << std::endl;
         }
       } else {
         generalized_layer = 10 * thedetId.subdetId();
         if (thedetId.subdetId() == StripSubdetector::TOB) {
           generalized_layer += tTopo_->tobLayer(thedetId.rawId());
         }
       }
       if (generalized_layer == 31) {
         ((TH1F*)HlistOtherHistos->FindObject("APVPairCorrectionsTIB1mono"))->Fill(local_correction);
       }
       if (generalized_layer == 32) {
         ((TH1F*)HlistOtherHistos->FindObject("APVPairCorrectionsTIB1stereo"))->Fill(local_correction);
       }
       if (generalized_layer == 33) {
         ((TH1F*)HlistOtherHistos->FindObject("APVPairCorrectionsTIB2"))->Fill(local_correction);
       }
       if (generalized_layer == 51) {
         ((TH1F*)HlistOtherHistos->FindObject("APVPairCorrectionsTOB1"))->Fill(local_correction);
       }
       if (generalized_layer == 52) {
         ((TH1F*)HlistOtherHistos->FindObject("APVPairCorrectionsTOB2"))->Fill(local_correction);
       }
       // control view
       const FedChannelConnection& fedchannelconnection =
           siStripDetCabling_->getConnection(extracted_detid, extracted_apvpairid);
       std::ostringstream local_key;
       // in S. Mersi's analysis the APVPair id seems to be used instead of the lldChannel, hence use the same here
       local_key << "fecCrate" << fedchannelconnection.fecCrate() << "_fecSlot" << fedchannelconnection.fecSlot()
                 << "_fecRing" << fedchannelconnection.fecRing() << "_ccuAddr" << fedchannelconnection.ccuAddr()
                 << "_ccuChan" << fedchannelconnection.ccuChan() << "_apvPair" << extracted_apvpairid;
       TString control_key = local_key.str();
       ChargeOfEachAPVPairControlView->Fill(control_key, local_charge_over_path);
       int ibin1 = ChargeOfEachAPVPairControlView->GetXaxis()->FindBin(control_key);
       ChargeOfEachAPVPairControlView->SetBinError(ibin1, local_error_of_charge);
       CorrectionOfEachAPVPairControlView->Fill(control_key, local_correction);
       int ibin2 = CorrectionOfEachAPVPairControlView->GetXaxis()->FindBin(control_key);
       CorrectionOfEachAPVPairControlView->SetBinError(ibin2, local_error_correction);
       // thickness of each module
       double module_thickness = moduleThickness(extracted_detid);
       if (fabs(module_thickness - 0.032) < 0.001)
         ModuleThickness->Fill(1);
       if (fabs(module_thickness - 0.05) < 0.001)
         ModuleThickness->Fill(2);
       // width of each module
       double module_width = moduleWidth(extracted_detid);
       if (fabs(module_width - 6.144) < 0.01)
         ModuleWidth->Fill(1);
       if (fabs(module_width - 7.14) < 0.01)
         ModuleWidth->Fill(2);
       if (fabs(module_width - 9.3696) < 0.01)
         ModuleWidth->Fill(3);
       if (fabs(module_width - 10.49) < 0.01)
         ModuleWidth->Fill(4);
       if (fabs(module_width - 12.03) < 0.01)
         ModuleWidth->Fill(5);
     }
   }
   HlistOtherHistos->Add(CorrectionOfEachAPVPair);
   ChargeOfEachAPVPairControlView->LabelsDeflate("X");
   CorrectionOfEachAPVPairControlView->LabelsDeflate("X");
   HlistOtherHistos->Add(ChargeOfEachAPVPairControlView);
   HlistOtherHistos->Add(CorrectionOfEachAPVPairControlView);
   // output histograms to file
 
   if (outputHistogramsInRootFile) {
     TFile* outputfile = new TFile(outputFileName, "RECREATE");
     HlistAPVPairs->Write();
     HlistOtherHistos->Write();
     outputfile->Close();
   }
 
   auto obj = std::make_unique<SiStripApvGain>();
 
   //   for(std::map<uint32_t,OptoScanAnalysis*>::const_iterator it = analyses.begin(); it != analyses.end(); it++){
   //     //Generate Gain for det detid
   //     std::vector<float> theSiStripVector;
   //     for(unsigned short j=0; j<it->second; j++){
   //       float gain;
 
   //       //      if(sigmaGain_/meanGain_ < 0.00001) gain = meanGain_;
   //       //      else{
   //       gain = CLHEP::RandGauss::shoot(meanGain_, sigmaGain_);
   //       if(gain<=minimumPosValue_) gain=minimumPosValue_;
   //       //      }
 
   //       if (printdebug_)
   //    edm::LogInfo("SiStripGainCalculator") << "detid " << it->first << " \t"
   //                          << " apv " << j << " \t"
   //                          << gain    << " \t"
   //                          << std::endl;
   //       theSiStripVector.push_back(gain);
   //     }
   //     SiStripApvGain::Range range(theSiStripVector.begin(),theSiStripVector.end());
   //     if ( ! obj->put(it->first,range) )
   //       edm::LogError("SiStripGainCalculator")<<"[SiStripGainCalculator::beginJob] detid already exists"<<std::endl;
   //   }
 
   return obj;
 }

◆ getPeakOfLandau()

std::pair< double, double > SiStripGainCosmicCalculator::getPeakOfLandau ( TH1F * inputHisto )

private

Definition at line 213 of file SiStripGainCosmicCalculator.cc.

References nano_mu_local_reco_cff::chi2, gather_cfg::cout, relativeConstraints::error, MaxChi2OverNDF, and MinNrEntries.

Referenced by getNewObject().

                       {  // automated fitting with finding of the appropriate nr. of ADCs
   // set some default dummy value and return if no entries
   double adcs = -0.5;
   double error = 0.;
   double nr_of_entries = inputHisto->GetEntries();
   if (nr_of_entries < MinNrEntries) {
     return std::make_pair(adcs, error);
   }
   //
   //  // fit with initial setting of  parameter values
   //  double rms_of_histogram = inputHisto->GetRMS();
   //  TF1 *landaufit = new TF1("landaufit","landau",0.,450.);
   //  landaufit->SetParameters(nr_of_entries,mean_of_histogram,rms_of_histogram);
   //  inputHisto->Fit("landaufit","0Q+");
   //  delete landaufit;
   //
   // perform fit with standard landau
   // make our own copy to avoid problems with threads
   std::unique_ptr<TF1> fitfunction(new TF1("landaufit", "landau"));
   inputHisto->Fit(fitfunction.get(), "0Q");
   adcs = fitfunction->GetParameter("MPV");
   error = fitfunction->GetParError(1);  // MPV is parameter 1 (0=constant, 1=MPV, 2=Sigma)
   double chi2 = fitfunction->GetChisquare();
   double ndf = fitfunction->GetNDF();
   double chi2overndf = chi2 / ndf;
   // in case things went wrong, try to refit in smaller range
   if (adcs < 2. || (error / adcs) > 1.8) {
     inputHisto->Fit(fitfunction.get(), "0Q", nullptr, 0., 400.);
     std::cout << "refitting landau for histogram " << inputHisto->GetTitle() << std::endl;
     std::cout << "initial error/adcs =" << error << " / " << adcs << std::endl;
     std::cout << "new     error/adcs =" << fitfunction->GetParError(1) << " / " << fitfunction->GetParameter("MPV")
               << std::endl;
     adcs = fitfunction->GetParameter("MPV");
     error = fitfunction->GetParError(1);  // MPV is parameter 1 (0=constant, 1=MPV, 2=Sigma)
     chi2 = fitfunction->GetChisquare();
     ndf = fitfunction->GetNDF();
     chi2overndf = chi2 / ndf;
   }
   // if still wrong, give up
   if (adcs < 2. || chi2overndf > MaxChi2OverNDF) {
     adcs = -0.5;
     error = 0.;
   }
   return std::make_pair(adcs, error);
 }

◆ moduleThickness()

double SiStripGainCosmicCalculator::moduleThickness ( const uint32_t detid )

private

Definition at line 274 of file SiStripGainCosmicCalculator.cc.

References gather_cfg::cout, TrackerGeometry::idToDetUnit(), ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::it, and tkGeom_.

Referenced by algoAnalyze(), and getNewObject().

 {                                                                          //dk: copied from A. Giammanco and hacked
   double module_thickness = 0.;
   const GeomDetUnit* it = tkGeom_->idToDetUnit(DetId(detid));
   if (dynamic_cast<const StripGeomDetUnit*>(it) == nullptr && dynamic_cast<const PixelGeomDetUnit*>(it) == nullptr) {
     std::cout << "this detID doesn't seem to belong to the Tracker" << std::endl;
   } else {
     module_thickness = it->surface().bounds().thickness();
   }
   return module_thickness;
 }

◆ moduleWidth()

double SiStripGainCosmicCalculator::moduleWidth ( const uint32_t detid )

private

Definition at line 261 of file SiStripGainCosmicCalculator.cc.

References gather_cfg::cout, TrackerGeometry::idToDetUnit(), ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::it, and tkGeom_.

Referenced by algoAnalyze(), and getNewObject().

 {  //dk: copied from A. Giammanco and hacked,  module_width values : 10.49 12.03 6.144 7.14 9.3696
   double module_width = 0.;
   const GeomDetUnit* it = tkGeom_->idToDetUnit(DetId(detid));
   if (dynamic_cast<const StripGeomDetUnit*>(it) == nullptr && dynamic_cast<const PixelGeomDetUnit*>(it) == nullptr) {
     std::cout << "this detID doesn't seem to belong to the Tracker" << std::endl;
   } else {
     module_width = it->surface().bounds().width();
   }
   return module_width;
 }