#include <SiStripGainCosmicCalculator.h>

Inheritance diagram for SiStripGainCosmicCalculator:

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

	~SiStripGainCosmicCalculator ()

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

virtual	~ConditionDBWriter ()

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

	EDAnalyzer ()

ModuleDescription const &	moduleDescription () const

std::string	workerType () const

virtual	~EDAnalyzer ()

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

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::vector< ModuleDescription const * > &modules, 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	updateLookup (BranchType iBranchType, ProductResolverIndexHelper const &, bool iPrefetchMayGet)

virtual	~EDConsumerBase () noexcept(false)

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

void	algoBeginJob (const edm::EventSetup &)

void	algoEndJob ()

SiStripApvGain *	getNewObject ()

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

double	moduleThickness (const uint32_t detid, const edm::EventSetup *iSetup)

double	moduleWidth (const uint32_t detid, const edm::EventSetup *iSetup)

Private Attributes
std::vector< uint32_t >	detModulesToBeExcluded

const edm::EventSetup *	eventSetupCopy_

double	ExpectedChargeDeposition

TObjArray *	HlistAPVPairs

TObjArray *	HlistOtherHistos

double	MaxChi2OverNDF

unsigned int	MinNrEntries

TString	outputFileName

bool	outputHistogramsInRootFile

bool	printdebug_

std::vector< uint32_t >	SelectedDetIds

std::map< uint32_t, double >	thickness_map

uint32_t	total_nr_of_events

std::string	TrackLabel

std::string	TrackProducer

const TrackerTopology *	tTopo

Additional Inherited Members
Public Types inherited from edm::EDAnalyzer
typedef EDAnalyzer	ModuleType

Public Types inherited from edm::EDConsumerBase
typedef ProductLabels	Labels

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

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &)

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)

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 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)

Detailed Description

Definition at line 25 of file SiStripGainCosmicCalculator.h.

Constructor & Destructor Documentation

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

explicit

Definition at line 37 of file SiStripGainCosmicCalculator.cc.

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

                                                                                        : 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);
   tTopo = nullptr;
 }

SiStripGainCosmicCalculator::~SiStripGainCosmicCalculator ( )

Definition at line 63 of file SiStripGainCosmicCalculator.cc.

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

Member Function Documentation

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

privatevirtual

Reimplemented from ConditionDBWriter< SiStripApvGain >.

Definition at line 145 of file SiStripGainCosmicCalculator.cc.

References SiStripRecHit2D::cluster(), funct::cos(), HcalObjRepresent::Fill(), edm::Event::getByLabel(), HlistAPVPairs, HlistOtherHistos, createfilelist::int, TrackingRecHit::localPosition(), moduleThickness(), moduleWidth(), edm::Handle< T >::product(), total_nr_of_events, findElectronsInSiStrips_cfi::trackCollection, TrackLabel, l1t::tracks, 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, &iSetup);
         ((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, &iSetup);
         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.);
       }
     }
   }
 }

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

privatevirtual

Reimplemented from ConditionDBWriter< SiStripApvGain >.

Definition at line 70 of file SiStripGainCosmicCalculator.cc.

References SiStripDetCabling::addActiveDetectorsRawIds(), gather_cfg::cout, detModulesToBeExcluded, TrackerGeometry::dets(), eventSetupCopy_, cmsRelvalreport::exit, edm::EventSetup::get(), SiStripSubStructure::getTIBDetectors(), SiStripSubStructure::getTOBDetectors(), HlistAPVPairs, HlistOtherHistos, StripTopology::nstrips(), AlCaHLTBitMon_ParallelJobs::p, edm::ESHandle< T >::product(), SelectedDetIds, thickness_map, total_nr_of_events, and tTopo.

 {
    //Retrieve tracker topology from geometry
    edm::ESHandle<TrackerTopology> tTopoHandle;
    iSetup.get<TrackerTopologyRcd>().get(tTopoHandle);
    tTopo = tTopoHandle.product();
  
    eventSetupCopy_ = &iSetup;
    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
    edm::ESHandle<SiStripDetCabling> siStripDetCabling; iSetup.get<SiStripDetCablingRcd>().get(siStripDetCabling);
    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 substructure;
    substructure.getTIBDetectors(activeDets, SelectedDetIds, 0, 0, 0, 0); // this adds rawDetIds to SelectedDetIds
    substructure.getTOBDetectors(activeDets, SelectedDetIds, 0, 0, 0);    // this adds rawDetIds to SelectedDetIds
    // get tracker geometry and find nr. of apv pairs for each active detector 
    edm::ESHandle<TrackerGeometry> tkGeom; iSetup.get<TrackerDigiGeometryRecord>().get( tkGeom );     
    for(TrackerGeometry::DetContainer::const_iterator it = tkGeom->dets().begin(); it != tkGeom->dets().end(); it++){ // loop over detector modules
      if( dynamic_cast<const StripGeomDetUnit*>((*it))!=0){
        uint32_t detid= ((*it)->geographicalId()).rawId();
        // get thickness for all detector modules, not just for active, this is strange 
        double module_thickness = (*it)->surface().bounds().thickness(); // get thickness of detector from GeomDet (DetContainer == vector<GeomDet*>)
        thickness_map.insert(std::make_pair(detid,module_thickness));
        //
        bool is_active_detector = false;
        for(std::vector<uint32_t>::iterator iactive = SelectedDetIds.begin(); iactive != SelectedDetIds.end(); iactive++){
          if( *iactive == detid ){
            is_active_detector = true;
            break; // leave for loop if found matching detid
          }
        }
        //
        bool exclude_this_detid = false;
        for( std::vector<uint32_t>::const_iterator imod = detModulesToBeExcluded.begin(); imod != detModulesToBeExcluded.end(); imod++ ){
            if(*imod == detid) exclude_this_detid = true; // found in exclusion list
            break;
        }
        //
        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*>((*it))->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
          }
        }
      }
    }
 }

void SiStripGainCosmicCalculator::algoEndJob ( )

privatevirtual

Reimplemented from ConditionDBWriter< SiStripApvGain >.

Definition at line 67 of file SiStripGainCosmicCalculator.cc.

67 {

68 }

SiStripApvGain * SiStripGainCosmicCalculator::getNewObject ( )

privatevirtual

Implements ConditionDBWriter< SiStripApvGain >.

Definition at line 280 of file SiStripGainCosmicCalculator.cc.

References FedChannelConnection::ccuAddr(), FedChannelConnection::ccuChan(), gather_cfg::cout, eventSetupCopy_, FedChannelConnection::fecCrate(), FedChannelConnection::fecRing(), FedChannelConnection::fecSlot(), HcalObjRepresent::Fill(), edm::EventSetup::get(), SiStripDetCabling::getConnection(), getPeakOfLandau(), HlistAPVPairs, HlistOtherHistos, moduleThickness(), moduleWidth(), MuonAssociatorByHits_cfi::obj, estimatePileup_makeJSON::outputfile, outputFileName, outputHistogramsInRootFile, DetId::rawId(), 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
        edm::ESHandle<SiStripDetCabling> siStripDetCabling; eventSetupCopy_->get<SiStripDetCablingRcd>().get(siStripDetCabling);
        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, eventSetupCopy_);
        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, eventSetupCopy_);
        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();
   }
 
   SiStripApvGain * obj = new 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;
 }

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

private

Definition at line 209 of file SiStripGainCosmicCalculator.cc.

References HiEvtPlane_cfi::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",0,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);
 }

double SiStripGainCosmicCalculator::moduleThickness	(	const uint32_t	detid,
		const edm::EventSetup *	iSetup
	)

private

Definition at line 266 of file SiStripGainCosmicCalculator.cc.

References Surface::bounds(), gather_cfg::cout, edm::EventSetup::get(), TrackerGeometry::idToDetUnit(), GeomDet::surface(), and Bounds::thickness().

Referenced by algoAnalyze(), and getNewObject().

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

double SiStripGainCosmicCalculator::moduleWidth	(	const uint32_t	detid,
		const edm::EventSetup *	iSetup
	)

private

Definition at line 252 of file SiStripGainCosmicCalculator.cc.

References Surface::bounds(), gather_cfg::cout, edm::EventSetup::get(), TrackerGeometry::idToDetUnit(), GeomDet::surface(), 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
   edm::ESHandle<TrackerGeometry> tkGeom; iSetup->get<TrackerDigiGeometryRecord>().get( tkGeom );     
   double module_width=0.;
   const GeomDetUnit* it = tkGeom->idToDetUnit(DetId(detid));
   if (dynamic_cast<const StripGeomDetUnit*>(it)==0 && dynamic_cast<const PixelGeomDetUnit*>(it)==0) {
     std::cout << "this detID doesn't seem to belong to the Tracker" << std::endl;
   }else{
     module_width = it->surface().bounds().width();
   }
   return module_width;
 }