72 static constexpr float defaultGainTick = 690. / 640.;
80 edm::LogError(
"SiStripGainPCLHarvester") <<
"gainHandle is not valid\n";
93 APV->isMasked = SiStripQuality_->
IsApvBad(APV->DetId, APV->APVId);
96 edm::LogError(
"SiStripGainPCLHarvester") <<
"NUMBER OF GAIN TAG IS EXPECTED TO BE 2\n";
100 float newPreviousGain = gainHandle->
getApvGain(APV->APVId, gainHandle->
getRange(APV->DetId, 1), 1);
101 if (APV->PreviousGain != 1 and newPreviousGain != APV->PreviousGain)
102 edm::LogWarning(
"SiStripGainPCLHarvester") <<
"WARNING: ParticleGain in the global tag changed\n";
103 APV->PreviousGain = newPreviousGain;
105 float newPreviousGainTick =
106 APV->isMasked ? defaultGainTick : gainHandle->
getApvGain(APV->APVId, gainHandle->
getRange(APV->DetId, 0), 0);
107 if (APV->PreviousGainTick != 1 and newPreviousGainTick != APV->PreviousGainTick) {
109 <<
"WARNING: TickMarkGain in the global tag changed\n" 111 <<
" APV->SubDet: " << APV->SubDet <<
" APV->APVId:" << APV->APVId << std::endl
112 <<
" APV->PreviousGainTick: " << APV->PreviousGainTick <<
" newPreviousGainTick: " << newPreviousGainTick
115 APV->PreviousGainTick = newPreviousGainTick;
121 edm::LogInfo(
"SiStripGainsPCLHarvester") <<
"Starting harvesting statistics" << std::endl;
126 if (!stag.empty() && stag[0] !=
'_')
127 stag.insert(0, 1,
'_');
133 if (Charge_Vs_Index ==
nullptr) {
135 <<
"Harvesting: could not retrieve " << cvi.c_str() <<
", statistics will not be summed!" << std::endl;
138 <<
"Harvesting " << (Charge_Vs_Index)->getTH2S()->GetEntries() <<
" more clusters" << std::endl;
142 std::unique_ptr<SiStripApvGain> theAPVGains = this->
getNewObject(Charge_Vs_Index);
151 throw std::runtime_error(
"PoolDBService required.");
153 edm::LogInfo(
"SiStripGainsPCLHarvester") <<
"Will not produce payload!" << std::endl;
165 std::vector<APVGain::APVmon> new_charge_histos;
166 std::vector<std::pair<std::string, std::string>> cnames =
168 for (
unsigned int i = 0;
i < cnames.size();
i++) {
173 new_charge_histos.push_back(
APVGain::APVmon(
id, side, plane, monitor));
181 ibooker_.
book2DD(
"MPVvsEtaTIB",
"MPV vs Eta TIB", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
183 ibooker_.
book2DD(
"MPVvsEtaTID",
"MPV vs Eta TID", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
185 ibooker_.
book2DD(
"MPVvsEtaTOB",
"MPV vs Eta TOB", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
187 ibooker_.
book2DD(
"MPVvsEtaTEC",
"MPV vs Eta TEC", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
189 ibooker_.
book2DD(
"MPVvsEtaTEC1",
"MPV vs Eta TEC-thin", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
191 ibooker_.
book2DD(
"MPVvsEtaTEC2",
"MPV vs Eta TEC-thick", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
194 ibooker_.
book2DD(
"MPVvsPhiTIB",
"MPV vs Phi TIB", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
196 ibooker_.
book2DD(
"MPVvsPhiTID",
"MPV vs Phi TID", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
198 ibooker_.
book2DD(
"MPVvsPhiTOB",
"MPV vs Phi TOB", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
200 ibooker_.
book2DD(
"MPVvsPhiTEC",
"MPV vs Phi TEC", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
202 ibooker_.
book2DD(
"MPVvsPhiTEC1",
"MPV vs Phi TEC-thin ", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
204 ibooker_.
book2DD(
"MPVvsPhiTEC2",
"MPV vs Phi TEC-thick", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
206 MonitorElement* NoMPVfit = ibooker_.
book2DD(
"NoMPVfit",
"Modules with bad Landau Fit", 350, -350, 350, 240, 0, 120);
207 MonitorElement* NoMPVmasked = ibooker_.
book2DD(
"NoMPVmasked",
"Masked Modules", 350, -350, 350, 240, 0, 120);
222 MonitorElement* MPVsTECthick = ibooker_.
book1DD(
"MPV_TEC2",
"MPV TEC thick", MPVbin, MPVmin, MPVmax);
229 MonitorElement* MPVErrorVsMPV = ibooker_.
book2DD(
"MPVErrorVsMPV",
"MPV Error vs MPV", 300, 0, 600, 150, 0, 150);
230 MonitorElement* MPVErrorVsEta = ibooker_.
book2DD(
"MPVErrorVsEta",
"MPV Error vs Eta", 50, -3.0, 3.0, 150, 0, 150);
231 MonitorElement* MPVErrorVsPhi = ibooker_.
book2DD(
"MPVErrorVsPhi",
"MPV Error vs Phi", 50, -3.4, 3.4, 150, 0, 150);
232 MonitorElement* MPVErrorVsN = ibooker_.
book2DD(
"MPVErrorVsN",
"MPV Error vs N", 500, 0, 1000, 150, 0, 150);
234 MonitorElement* DiffWRTPrevGainTIB = ibooker_.
book1DD(
"DiffWRTPrevGainTIB",
"Diff w.r.t. PrevGain TIB", 250, 0, 2);
235 MonitorElement* DiffWRTPrevGainTID = ibooker_.
book1DD(
"DiffWRTPrevGainTID",
"Diff w.r.t. PrevGain TID", 250, 0, 2);
236 MonitorElement* DiffWRTPrevGainTOB = ibooker_.
book1DD(
"DiffWRTPrevGainTOB",
"Diff w.r.t. PrevGain TOB", 250, 0, 2);
237 MonitorElement* DiffWRTPrevGainTEC = ibooker_.
book1DD(
"DiffWRTPrevGainTEC",
"Diff w.r.t. PrevGain TEC", 250, 0, 2);
240 ibooker_.
book2DD(
"GainVsPrevGainTIB",
"Gain vs PrevGain TIB", 100, 0, 2, 100, 0, 2);
242 ibooker_.
book2DD(
"GainVsPrevGainTID",
"Gain vs PrevGain TID", 100, 0, 2, 100, 0, 2);
244 ibooker_.
book2DD(
"GainVsPrevGainTOB",
"Gain vs PrevGain TOB", 100, 0, 2, 100, 0, 2);
246 ibooker_.
book2DD(
"GainVsPrevGainTEC",
"Gain vs PrevGain TEC", 100, 0, 2, 100, 0, 2);
253 unsigned int Index = APV->Index;
254 unsigned int SubDet = APV->SubDet;
257 float Eta = APV->Eta;
259 float Phi = APV->Phi;
260 float Thickness = APV->Thickness;
261 double FitMPV = APV->FitMPV;
262 double FitMPVErr = APV->FitMPVErr;
263 double Gain = APV->Gain;
264 double NEntries = APV->NEntries;
265 double PreviousGain = APV->PreviousGain;
272 TH2S* chvsidx = (Charge_Vs_Index)->getTH2S();
273 int bin = chvsidx->GetXaxis()->FindBin(Index);
274 TH1D* Proj = chvsidx->ProjectionY(
"proj", bin, bin);
275 for (
int binId = 0; binId < Proj->GetXaxis()->GetNbins(); binId++) {
276 double new_charge = Proj->GetXaxis()->GetBinCenter(binId) /
Gain;
277 if (Proj->GetBinContent(binId) != 0.) {
278 for (
unsigned int h = 0;
h < charge_histos.size();
h++) {
279 TH1D* chisto = (charge_histos[
h])->getTH1D();
280 for (
int e = 0;
e < Proj->GetBinContent(binId);
e++)
281 chisto->Fill(new_charge);
289 NoMPVmasked->
Fill(z, R);
291 NoMPVfit->
Fill(z, R);
298 if (Thickness < 0.04)
299 MPVs320->
Fill(FitMPV);
300 if (Thickness > 0.04)
301 MPVs500->
Fill(FitMPV);
303 MPVError->
Fill(FitMPVErr);
304 MPVErrorVsMPV->
Fill(FitMPV, FitMPVErr);
305 MPVErrorVsEta->
Fill(Eta, FitMPVErr);
306 MPVErrorVsPhi->
Fill(Phi, FitMPVErr);
307 MPVErrorVsN->
Fill(NEntries, FitMPVErr);
310 MPV_Vs_EtaTIB->
Fill(Eta, FitMPV);
311 MPV_Vs_PhiTIB->
Fill(Phi, FitMPV);
312 MPVsTIB->
Fill(FitMPV);
314 }
else if (SubDet == 4) {
315 MPV_Vs_EtaTID->
Fill(Eta, FitMPV);
316 MPV_Vs_PhiTID->
Fill(Phi, FitMPV);
317 MPVsTID->
Fill(FitMPV);
319 MPVsTIDM->
Fill(FitMPV);
321 MPVsTIDP->
Fill(FitMPV);
323 }
else if (SubDet == 5) {
324 MPV_Vs_EtaTOB->
Fill(Eta, FitMPV);
325 MPV_Vs_PhiTOB->
Fill(Phi, FitMPV);
326 MPVsTOB->
Fill(FitMPV);
328 }
else if (SubDet == 6) {
329 MPV_Vs_EtaTEC->
Fill(Eta, FitMPV);
330 MPV_Vs_PhiTEC->
Fill(Phi, FitMPV);
331 MPVsTEC->
Fill(FitMPV);
333 MPVsTECM->
Fill(FitMPV);
335 MPVsTECP->
Fill(FitMPV);
336 if (Thickness < 0.04) {
337 MPV_Vs_EtaTECthin->
Fill(Eta, FitMPV);
338 MPV_Vs_PhiTECthin->
Fill(Phi, FitMPV);
339 MPVsTECthin->
Fill(FitMPV);
341 MPVsTECP1->
Fill(FitMPV);
343 MPVsTECM1->
Fill(FitMPV);
345 if (Thickness > 0.04) {
346 MPV_Vs_EtaTECthick->
Fill(Eta, FitMPV);
347 MPV_Vs_PhiTECthick->
Fill(Phi, FitMPV);
348 MPVsTECthick->
Fill(FitMPV);
350 MPVsTECP2->
Fill(FitMPV);
352 MPVsTECM2->
Fill(FitMPV);
357 if (SubDet == 3 && PreviousGain != 0.)
358 DiffWRTPrevGainTIB->
Fill(Gain / PreviousGain);
359 else if (SubDet == 4 && PreviousGain != 0.)
360 DiffWRTPrevGainTID->
Fill(Gain / PreviousGain);
361 else if (SubDet == 5 && PreviousGain != 0.)
362 DiffWRTPrevGainTOB->
Fill(Gain / PreviousGain);
363 else if (SubDet == 6 && PreviousGain != 0.)
364 DiffWRTPrevGainTEC->
Fill(Gain / PreviousGain);
367 GainVsPrevGainTIB->
Fill(PreviousGain, Gain);
368 else if (SubDet == 4)
369 GainVsPrevGainTID->
Fill(PreviousGain, Gain);
370 else if (SubDet == 5)
371 GainVsPrevGainTOB->
Fill(PreviousGain, Gain);
372 else if (SubDet == 6)
373 GainVsPrevGainTEC->
Fill(PreviousGain, Gain);
380 TH1F* Proj =
nullptr;
381 double FitResults[6];
382 double MPVmean = 300;
384 if (Charge_Vs_Index ==
nullptr) {
386 <<
"Harvesting: could not execute algoComputeMPVandGain method because " <<
m_calibrationMode 387 <<
" statistics cannot be retrieved.\n" 388 <<
"Please check if input contains " <<
m_calibrationMode <<
" data." << std::endl;
392 TH2S* chvsidx = (Charge_Vs_Index)->getTH2S();
394 printf(
"Progressing Bar :0%% 20%% 40%% 60%% 80%% 100%%\n");
395 printf(
"Fitting Charge Distribution :");
396 int TreeStep =
APVsColl.size() / 50;
399 if (I % TreeStep == 0) {
403 std::shared_ptr<stAPVGain>
APV = it->second;
405 APV->Bin = chvsidx->GetXaxis()->FindBin(APV->Index);
408 APV->Gain = APV->PreviousGain;
413 Proj = (TH1F*)(chvsidx->ProjectionY(
414 "", chvsidx->GetXaxis()->FindBin(APV->Index), chvsidx->GetXaxis()->FindBin(APV->Index),
"e"));
420 int SecondAPVId = APV->APVId;
421 if (SecondAPVId % 2 == 0) {
422 SecondAPVId = SecondAPVId + 1;
424 SecondAPVId = SecondAPVId - 1;
426 std::shared_ptr<stAPVGain> APV2 =
APVsColl[(APV->DetId << 4) | SecondAPVId];
428 APV2->Bin = chvsidx->GetXaxis()->FindBin(APV2->Index);
429 TH1F* Proj2 = (TH1F*)(chvsidx->ProjectionY(
"", APV2->Bin, APV2->Bin,
"e"));
435 for (
unsigned int i = 0;
i < 16;
i++) {
436 auto tmpit =
APVsColl.find((APV->DetId << 4) |
i);
439 std::shared_ptr<stAPVGain> APV2 = tmpit->second;
440 if (APV2->DetId != APV->DetId || APV2->APVId == APV->APVId)
443 APV2->Bin = chvsidx->GetXaxis()->FindBin(APV2->Index);
444 TH1F* Proj2 = (TH1F*)(chvsidx->ProjectionY(
"", APV2->Bin, APV2->Bin,
"e"));
456 APV->FitMPV = FitResults[0];
457 APV->FitMPVErr = FitResults[1];
458 APV->FitWidth = FitResults[2];
459 APV->FitWidthErr = FitResults[3];
460 APV->FitChi2 = FitResults[4];
461 APV->FitNorm = FitResults[5];
462 APV->NEntries = Proj->GetEntries();
465 APV->Gain = APV->FitMPV / MPVmean;
469 APV->Gain = APV->PreviousGain;
483 FitResults[0] = -0.5;
485 FitResults[2] = -0.5;
487 FitResults[4] = -0.5;
494 TF1 MyLandau(
"MyLandau",
"landau", LowRange, HighRange);
495 MyLandau.SetParameter(1, 300);
496 InputHisto->Fit(&MyLandau,
"0QR WW");
499 FitResults[0] = MyLandau.GetParameter(1);
500 FitResults[1] = MyLandau.GetParError(1);
501 FitResults[2] = MyLandau.GetParameter(2);
502 FitResults[3] = MyLandau.GetParError(2);
503 FitResults[4] = MyLandau.GetChisquare() / MyLandau.GetNDF();
504 FitResults[5] = MyLandau.GetParameter(0);
509 if (FitResults[0] <= 0)
525 auto const&
Det = newBareTkGeomPtr->
dets();
527 unsigned int Index = 0;
529 for (
unsigned int i = 0;
i <
Det.size();
i++) {
540 unsigned int NAPV = Topo.
nstrips() / 128;
542 for (
unsigned int j = 0;
j < NAPV;
j++) {
543 auto APV = std::make_shared<stAPVGain>();
552 APV->FitWidthErr = -1;
556 APV->PreviousGain = 1;
557 APV->PreviousGainTick = 1;
558 APV->x = DetUnit->position().basicVector().x();
559 APV->y = DetUnit->position().basicVector().y();
560 APV->z = DetUnit->position().basicVector().z();
561 APV->Eta = DetUnit->position().basicVector().eta();
562 APV->Phi = DetUnit->position().basicVector().phi();
563 APV->R = DetUnit->position().basicVector().transverse();
564 APV->Thickness = DetUnit->surface().bounds().thickness();
566 APV->isMasked =
false;
576 for (
unsigned int i = 0;
i <
Det.size();
586 unsigned int NROCRow = Topo.
nrows() / (80.);
587 unsigned int NROCCol = Topo.
ncolumns() / (52.);
589 for (
unsigned int j = 0;
j < NROCRow;
j++) {
590 for (
unsigned int i = 0;
i < NROCCol;
i++) {
591 auto APV = std::make_shared<stAPVGain>();
595 APV->APVId = (
j << 3 |
i);
600 APV->FitWidthErr = -1;
603 APV->PreviousGain = 1;
604 APV->PreviousGainTick = 1;
605 APV->x = DetUnit->position().basicVector().x();
606 APV->y = DetUnit->position().basicVector().y();
607 APV->z = DetUnit->position().basicVector().z();
608 APV->Eta = DetUnit->position().basicVector().eta();
609 APV->Phi = DetUnit->position().basicVector().phi();
610 APV->R = DetUnit->position().basicVector().transverse();
611 APV->Thickness = DetUnit->surface().bounds().thickness();
612 APV->isMasked =
false;
641 if (Charge_Vs_Index ==
nullptr) {
643 <<
"produceTagFilter -> Return false: could not retrieve the " <<
m_calibrationMode <<
" statistics.\n" 644 <<
"Please check if input contains " <<
m_calibrationMode <<
" data." << std::endl;
648 float integral = (Charge_Vs_Index)->getTH2S()->Integral();
652 <<
"produceTagFilter -> Return false: Statistics is too low : " << integral << std::endl;
658 <<
"produceTagFilter -> Return false: ratio of GOOD/TOTAL is too low: " << (1.0 *
GOOD) / (
GOOD +
BAD)
667 std::unique_ptr<SiStripApvGain>
obj = std::unique_ptr<SiStripApvGain>(
new SiStripApvGain());
671 <<
"getNewObject -> will not produce a paylaod because produceTagFilter returned false " << std::endl;
677 std::vector<float> theSiStripVector;
678 unsigned int PreviousDetId = 0;
681 if (APV ==
nullptr) {
685 if (APV->SubDet <= 2)
687 if (APV->DetId != PreviousDetId) {
688 if (!theSiStripVector.empty()) {
690 if (!obj->put(PreviousDetId,
range))
691 printf(
"Bug to put detId = %i\n", PreviousDetId);
693 theSiStripVector.clear();
694 PreviousDetId = APV->DetId;
696 theSiStripVector.push_back(APV->Gain);
698 LogDebug(
"SiStripGainsPCLHarvester") <<
" DetId: " << APV->DetId <<
" APV: " << APV->APVId
699 <<
" Gain: " << APV->Gain << std::endl;
701 if (!theSiStripVector.empty()) {
703 if (!obj->put(PreviousDetId,
range))
704 printf(
"Bug to put detId = %i\n", PreviousDetId);
std::unique_ptr< SiStripApvGain > getNewObject(const MonitorElement *Charge_Vs_Index)
bool IsApvBad(const uint32_t &detid, const short &apvNb) const
static constexpr auto TEC
T getUntrackedParameter(std::string const &, T const &) const
FWCore Framework interface EventSetupRecordImplementation h
Helper function to determine trigger accepts.
virtual int nrows() const =0
edm::ESHandle< TrackerGeometry > tkGeom_
void dqmEndJob(DQMStore::IBooker &ibooker_, DQMStore::IGetter &igetter_) override
size_t getNumberOfTags() const
void setCurrentFolder(std::string const &fullpath)
MonitorElement * book1DD(TString const &name, TString const &title, int nchX, double lowX, double highX)
std::vector< MonitorElement * > FetchMonitor(std::vector< APVmon >, uint32_t, const TrackerTopology *topo=0)
constexpr uint32_t rawId() const
get the raw id
const TrackerTopology * tTopo_
int subdetectorPlane(uint32_t, const TrackerTopology *)
std::vector< std::string > dqm_tag_
void beginRun(edm::Run const &run, edm::EventSetup const &isetup) override
bool produceTagFilter(const MonitorElement *Charge_Vs_Index)
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
void getPeakOfLandau(TH1 *InputHisto, double *FitResults, double LowRange=50, double HighRange=5400)
void algoComputeMPVandGain(const MonitorElement *Charge_Vs_Index)
static float getApvGain(const uint16_t &apv, const SiStripApvGain::Range &range)
U second(std::pair< T, U > const &p)
SiStripGainsPCLHarvester(const edm::ParameterSet &ps)
const DetContainer & dets() const override
Returm a vector of all GeomDet (including all GeomDetUnits)
Data available (reporting)
void addDefault(ParameterSetDescription const &psetDescription)
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
double tagCondition_GoodFrac
bool doChargeMonitorPerPlane
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector's numbering enum) ...
std::pair< ContainerIterator, ContainerIterator > Range
const std::complex< double > I
void writeOne(T *payload, Time_t time, const std::string &recordName, bool withlogging=false)
static constexpr auto TOB
std::vector< std::string > VChargeHisto
Integral< F, X >::type integral(const F &f)
virtual void checkAndRetrieveTopology(const edm::EventSetup &setup)
std::string m_calibrationMode
bool IsGoodLandauFit(double *FitResults)
static constexpr auto TIB
void gainQualityMonitor(DQMStore::IBooker &ibooker_, const MonitorElement *Charge_Vs_Index) const
virtual int nstrips() const =0
int subdetectorSide(uint32_t, const TrackerTopology *)
const TrackerGeometry * bareTkGeomPtr_
void endRun(edm::Run const &run, edm::EventSetup const &isetup) override
double tagCondition_NClusters
virtual void checkBookAPVColls(const edm::EventSetup &setup)
std::vector< std::shared_ptr< stAPVGain > > APVsCollOrdered
virtual int ncolumns() const =0
MonitorElement * book2DD(TString const &name, TString const &title, int nchX, double lowX, double highX, int nchY, double lowY, double highY)
cond::Time_t currentTime() const
std::vector< std::pair< std::string, std::string > > monHnames(std::vector< std::string >, bool, const char *tag)
MonitorElement * get(std::string const &path)
static constexpr auto TID
T const * product() const
int subdetectorId(uint32_t)
constexpr DetId()
Create an empty or null id (also for persistence)
const SiStripApvGain::Range getRange(uint32_t detID) const
std::unordered_map< unsigned int, std::shared_ptr< stAPVGain > > APVsColl