#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 final

bool	wantsGlobalRuns () const final

bool	wantsInputProcessBlocks () const final

bool	wantsProcessBlocks () const 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

bool	wantsStreamRuns () const

	~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

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

std::vector< ESProxyIndex > 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

bool	registeredToConsumeMany (TypeID const &, 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::ESRecordsToProxyIndices 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 ProductType , BranchType B = InEvent>
void	consumesMany ()

void	consumesMany (const TypeToGet &id)

template<BranchType B>
void	consumesMany (const TypeToGet &id)

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, duplicaterechits_cfi::trackCollection, TrackLabel, 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(), mps_fire::end, 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 hltPixelTracks_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 Surface::bounds(), gather_cfg::cout, TrackerGeometry::idToDetUnit(), GeomDet::surface(), Bounds::thickness(), 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 Surface::bounds(), gather_cfg::cout, TrackerGeometry::idToDetUnit(), GeomDet::surface(), tkGeom_, and Bounds::width().

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;
 }