6 #include <unordered_map> 52 #include "TObjString.h" 73 void algoEndJob()
override;
78 std::unique_ptr<SiStripApvGain> getNewObject()
override;
85 double trajChi2OverN);
88 void getPeakOfLandau(TH1* InputHisto,
double* FitResults,
double LowRange = 0,
double HighRange = 5400);
266 return PseudoDetId1 == PseudoDetId2;
271 std::unordered_map<unsigned int, stAPVGain*>
APVsColl;
306 if (strcmp(
AlgoMode.c_str(),
"WriteOnDB") == 0)
311 tTopoToken_ = esConsumes<edm::Transition::BeginRun>();
313 gainToken_ = esConsumes<edm::Transition::BeginRun>();
343 tmp =
dqmStore_->
book2D(
"APV_PathLength",
"APV_PathLength", 72785, 0, 72784, 100, 0.2, 1.4);
360 tmp =
dqmStore_->
book2D(
"Tracks_Pt_Vs_Eta",
"Tracks_Pt_Vs_Eta", 30, 0, 3, 50, 0, 200);
363 tmp =
dqmStore_->
book2D(
"Charge_Vs_PathTIB",
"Charge_Vs_PathTIB", 100, 0.2, 1.4, 500, 0, 2000);
365 tmp =
dqmStore_->
book2D(
"Charge_Vs_PathTID",
"Charge_Vs_PathTID", 100, 0.2, 1.4, 500, 0, 2000);
367 tmp =
dqmStore_->
book2D(
"Charge_Vs_PathTOB",
"Charge_Vs_PathTOB", 100, 0.2, 1.4, 500, 0, 2000);
369 tmp =
dqmStore_->
book2D(
"Charge_Vs_PathTEC",
"Charge_Vs_PathTEC", 100, 0.2, 1.4, 500, 0, 2000);
371 tmp =
dqmStore_->
book2D(
"Charge_Vs_PathTEC1",
"Charge_Vs_PathTEC1", 100, 0.2, 1.4, 500, 0, 2000);
373 tmp =
dqmStore_->
book2D(
"Charge_Vs_PathTEC2",
"Charge_Vs_PathTEC2", 100, 0.2, 1.4, 500, 0, 2000);
404 tmp =
dqmStore_->
book2D(
"Charge_Vs_PathLength",
"Charge_Vs_PathLength", 100, 0.2, 1.4, 500, 0, 2000);
406 tmp =
dqmStore_->
book2D(
"Charge_Vs_PathLength320",
"Charge_Vs_PathLength", 100, 0.2, 1.4, 500, 0, 2000);
408 tmp =
dqmStore_->
book2D(
"Charge_Vs_PathLength500",
"Charge_Vs_PathLength", 100, 0.2, 1.4, 500, 0, 2000);
411 tmp =
dqmStore_->
book2D(
"Charge_Vs_TransversAngle",
"Charge_Vs_TransversAngle", 220, -20, 200, 500, 0, 2000);
413 tmp =
dqmStore_->
book2D(
"Charge_Vs_Alpha",
"Charge_Vs_Alpha", 220, -20, 200, 500, 0, 2000);
415 tmp =
dqmStore_->
book2D(
"Charge_Vs_Beta",
"Charge_Vs_Beta", 220, -20, 200, 500, 0, 2000);
418 tmp =
dqmStore_->
book2D(
"NStrips_Vs_TransversAngle",
"NStrips_Vs_TransversAngle", 220, -20, 200, 10, 0, 10);
420 tmp =
dqmStore_->
book2D(
"NStrips_Vs_Alpha",
"NStrips_Vs_Alpha", 220, -20, 200, 10, 0, 10);
422 tmp =
dqmStore_->
book2D(
"NStrips_Vs_Beta",
"NStrips_Vs_Beta", 220, -20, 200, 10, 0, 10);
437 if (strcmp(
AlgoMode.c_str(),
"MultiJob") != 0) {
438 tmp =
dqmStore_->
book2D(
"MPV_Vs_EtaTIB",
"MPV_Vs_EtaTIB", 50, -3.0, 3.0, 300, 0, 600);
440 tmp =
dqmStore_->
book2D(
"MPV_Vs_EtaTID",
"MPV_Vs_EtaTID", 50, -3.0, 3.0, 300, 0, 600);
442 tmp =
dqmStore_->
book2D(
"MPV_Vs_EtaTOB",
"MPV_Vs_EtaTOB", 50, -3.0, 3.0, 300, 0, 600);
444 tmp =
dqmStore_->
book2D(
"MPV_Vs_EtaTEC",
"MPV_Vs_EtaTEC", 50, -3.0, 3.0, 300, 0, 600);
446 tmp =
dqmStore_->
book2D(
"MPV_Vs_EtaTEC1",
"MPV_Vs_EtaTEC1", 50, -3.0, 3.0, 300, 0, 600);
448 tmp =
dqmStore_->
book2D(
"MPV_Vs_EtaTEC2",
"MPV_Vs_EtaTEC2", 50, -3.0, 3.0, 300, 0, 600);
451 tmp =
dqmStore_->
book2D(
"MPV_Vs_PhiTIB",
"MPV_Vs_PhiTIB", 50, -3.2, 3.2, 300, 0, 600);
453 tmp =
dqmStore_->
book2D(
"MPV_Vs_PhiTID",
"MPV_Vs_PhiTID", 50, -3.2, 3.2, 300, 0, 600);
455 tmp =
dqmStore_->
book2D(
"MPV_Vs_PhiTOB",
"MPV_Vs_PhiTOB", 50, -3.2, 3.2, 300, 0, 600);
457 tmp =
dqmStore_->
book2D(
"MPV_Vs_PhiTEC",
"MPV_Vs_PhiTEC", 50, -3.2, 3.2, 300, 0, 600);
459 tmp =
dqmStore_->
book2D(
"MPV_Vs_PhiTEC1",
"MPV_Vs_PhiTEC1", 50, -3.2, 3.2, 300, 0, 600);
461 tmp =
dqmStore_->
book2D(
"MPV_Vs_PhiTEC2",
"MPV_Vs_PhiTEC2", 50, -3.2, 3.2, 300, 0, 600);
505 tmp =
dqmStore_->
book1D(
"FWHM_Vs_PathLength320",
"FWHM_Vs_PathLength", 100, 0.2, 1.4);
507 tmp =
dqmStore_->
book1D(
"FWHM_Vs_PathLength500",
"FWHM_Vs_PathLength", 100, 0.2, 1.4);
510 tmp =
dqmStore_->
book1D(
"MPV_Vs_TransversAngle",
"MPV_Vs_TransversAngle", 220, -20, 200);
519 tmp =
dqmStore_->
book2D(
"Error_Vs_Entries",
"Error_Vs_Entries", 500, 0, 10000, 100, 0, 50);
526 tmp =
dqmStore_->
book2D(
"NoMPV_Vs_EtaPhi",
"NoMPV_Vs_EtaPhi", 50, -3.0, 3.0, 50, -3.2, 3.2);
529 tmp =
dqmStore_->
book1D(
"NumberOfEntriesByAPV",
"NumberOfEntriesByAPV", 1000, 0, 10000);
550 printf(
"\n#####################\n\nERROR --> gainHandle is not valid\n\n#####################\n\n");
556 for (
unsigned int i = 0;
i <
Det.size();
i++) {
567 unsigned int NAPV = Topo.
nstrips() / 128;
569 double Phi = DetUnit->position().basicVector().phi();
570 double Eta = DetUnit->position().basicVector().eta();
571 double R = DetUnit->position().basicVector().transverse();
572 double Thick = DetUnit->surface().bounds().thickness();
574 for (
unsigned int j = 0;
j < NAPV;
j++) {
582 APV->PreviousGain = 1;
586 APV->Thickness = Thick;
621 if (strcmp(
AlgoMode.c_str(),
"WriteOnDB") == 0 || strcmp(
AlgoMode.c_str(),
"Merge") == 0) {
622 TH1::AddDirectory(kTRUE);
624 TFile*
file =
nullptr;
626 printf(
"Loading New Input File : %s\n",
VInputFiles[
f].c_str());
628 FILE* doesFileExist = fopen(
VInputFiles[
f].c_str(),
"r");
629 if (!doesFileExist) {
630 printf(
"File %s doesn't exist\n",
VInputFiles[
f].c_str());
633 fclose(doesFileExist);
638 printf(
"### Bug With File %s\n### File will be skipped \n",
VInputFiles[
f].c_str());
669 TH1F* JobInfo_tmp = (TH1F*)
file->FindObjectAny(
"JobInfo");
670 NEvent += (
unsigned int)JobInfo_tmp->GetBinContent(JobInfo_tmp->GetXaxis()->FindBin(1));
671 unsigned int tmp_SRun = (
unsigned int)JobInfo_tmp->GetBinContent(JobInfo_tmp->GetXaxis()->FindBin(3));
672 unsigned int tmp_SEvent = (
unsigned int)JobInfo_tmp->GetBinContent(JobInfo_tmp->GetXaxis()->FindBin(4));
673 unsigned int tmp_ERun = (
unsigned int)JobInfo_tmp->GetBinContent(JobInfo_tmp->GetXaxis()->FindBin(6));
674 unsigned int tmp_EEvent = (
unsigned int)JobInfo_tmp->GetBinContent(JobInfo_tmp->GetXaxis()->FindBin(7));
679 if (tmp_SRun <
SRun) {
682 }
else if (tmp_SRun ==
SRun && tmp_SEvent <
SEvent) {
686 if (tmp_ERun >
ERun) {
689 }
else if (tmp_ERun ==
ERun && tmp_EEvent >
EEvent) {
693 printf(
"Deleting Current Input File\n");
705 if (strcmp(
AlgoMode.c_str(),
"MultiJob") != 0) {
706 TH1D* Proj =
nullptr;
707 double FitResults[5];
711 printf(
"Fitting Histograms \t %6.2f%%\n", (100.0 *
I) /
APVsColl.size());
717 Proj = (TH1D*)Proj->Clone();
723 int SecondAPVId =
APV->APVId;
724 if (SecondAPVId % 2 == 0) {
725 SecondAPVId = SecondAPVId + 1;
727 SecondAPVId = SecondAPVId - 1;
732 TH1D* Proj2 =
APV_Charge->ProjectionY(
" ", bin2, bin2,
"e");
733 if (Proj2 !=
nullptr) {
747 TH1D* Proj2 =
APV_Charge->ProjectionY(
" ", bin2, bin2,
"e");
748 if (Proj2 !=
nullptr) {
759 APV->FitMPV = FitResults[0];
761 if (FitResults[0] != -0.5 && FitResults[1] <
MaxMPVError) {
764 if (
APV->Thickness < 0.04)
766 if (
APV->Thickness > 0.04)
779 if (
APV->Thickness < 0.04)
781 if (
APV->Thickness > 0.04)
793 if (
APV->Thickness < 0.04)
795 if (
APV->Thickness > 0.04)
812 if (
APV->Thickness < 0.04)
814 if (
APV->Thickness > 0.04)
836 if (
APV->Thickness < 0.04)
838 if (
APV->Thickness > 0.04)
846 if (FitResults[0] != -0.5) {
865 unsigned int GOOD = 0;
866 unsigned int BAD = 0;
867 double MPVmean = 300.;
870 if (
APV->FitMPV > 0) {
871 APV->Gain =
APV->FitMPV / MPVmean;
883 APV->Gain *=
APV->PreviousGain;
890 if (FitResults[0] == -0.5)
904 if (FitResults[0] == -0.5)
918 if (FitResults[0] == -0.5)
983 if (FitResults[0] == -0.5)
993 if (FitResults[0] == -0.5)
1003 if (FitResults[0] == -0.5)
1013 if (FitResults[0] == -0.5)
1023 if (FitResults[0] == -0.5)
1033 if (FitResults[0] == -0.5)
1043 if (FitResults[0] == -0.5)
1053 if (FitResults[0] == -0.5)
1063 if (FitResults[0] == -0.5)
1071 fprintf(Gains,
"NEvents = %i\n",
NEvent);
1072 fprintf(Gains,
"Number of APVs = %lu\n", static_cast<unsigned long>(
APVsColl.size()));
1073 fprintf(Gains,
"GoodFits = %i BadFits = %i ratio = %f\n",
GOOD,
BAD, (100.0 *
GOOD) / (
GOOD +
BAD));
1077 "%i | %i | PreviousGain = %7.5f NewGain = %7.5f\n",
1084 std::vector<unsigned int> DetIdOfBuggedAPV;
1085 fprintf(Gains,
"----------------------------------------------------------------------\n");
1088 if (
APV->FitMPV > 0 &&
APV->FitMPV < 200) {
1089 bool tmpBug =
false;
1090 for (
unsigned int b = 0;
b < DetIdOfBuggedAPV.size() && !tmpBug;
b++) {
1091 if (DetIdOfBuggedAPV[
b] ==
APV->DetId)
1095 fprintf(Gains,
"%i,\n",
APV->DetId);
1096 DetIdOfBuggedAPV.push_back(
APV->DetId);
1114 if (!gainHandle.isValid()) {
1115 printf(
"\n#####################\n\nERROR --> gainHandle is not valid\n\n#####################\n\n");
1122 if (gainHandle.isValid()) {
1124 APV->PreviousGain = *(detGainRange.first +
APV->APVId);
1127 if (
APV->PreviousGain < 0)
1128 APV->PreviousGain = 1;
1130 printf(
"GAIN HANDLE IS NOT VALID\n");
1137 if (strcmp(
AlgoMode.c_str(),
"WriteOnDB") == 0)
1171 for (Trajectory::RecHitContainer::const_iterator rechit = transRecHits.begin(); rechit != transRecHits.end();
1173 if ((*rechit)->isValid())
1174 ndof += (*rechit)->dimension();
1182 vector<TrajectoryMeasurement> measurements = traj.
measurements();
1206 for (vector<TrajectoryMeasurement>::const_iterator measurement_it = measurements.begin();
1207 measurement_it != measurements.end();
1218 if (sistripsimplehit) {
1221 }
else if (sistripmatchedhit) {
1225 }
else if (sistripsimple1dhit) {
1238 double trajChi2OverN) {
1240 double cosine = trackDirection.
z() / trackDirection.
mag();
1244 int APVId = FirstStrip / 128;
1247 bool Overlaping =
false;
1249 unsigned int NHighStrip = 0;
1254 if (FirstStrip == 0)
1256 if (FirstStrip == 128)
1258 if (FirstStrip == 256)
1260 if (FirstStrip == 384)
1262 if (FirstStrip == 512)
1264 if (FirstStrip == 640)
1267 if (FirstStrip <= 127 && FirstStrip + Ampls.size() > 127)
1269 if (FirstStrip <= 255 && FirstStrip + Ampls.size() > 255)
1271 if (FirstStrip <= 383 && FirstStrip + Ampls.size() > 383)
1273 if (FirstStrip <= 511 && FirstStrip + Ampls.size() > 511)
1275 if (FirstStrip <= 639 && FirstStrip + Ampls.size() > 639)
1278 if (FirstStrip + Ampls.size() == 127)
1280 if (FirstStrip + Ampls.size() == 255)
1282 if (FirstStrip + Ampls.size() == 383)
1284 if (FirstStrip + Ampls.size() == 511)
1286 if (FirstStrip + Ampls.size() == 639)
1288 if (FirstStrip + Ampls.size() == 767)
1302 for (
unsigned int a = 0;
a < Ampls.size();
a++) {
1304 if (Ampls[
a] >= 254)
1309 double path = (10.0 *
APV->Thickness) / fabs(cosine);
1310 double ClusterChargeOverPath = (double)
Charge /
path;
1319 if (
APV->Thickness < 0.04)
1321 if (
APV->Thickness > 0.04)
1331 if (
APV->Thickness < 0.04)
1333 if (
APV->Thickness > 0.04)
1337 double trans = atan2(trackDirection.
y(), trackDirection.
x()) * (180 / 3.14159265);
1339 acos(trackDirection.
x() /
sqrt(
pow(trackDirection.
x(), 2) +
pow(trackDirection.
z(), 2))) * (180 / 3.14159265);
1341 acos(trackDirection.
y() /
sqrt(
pow(trackDirection.
x(), 2) +
pow(trackDirection.
z(), 2))) * (180 / 3.14159265);
1363 ClusterChargeOverPath = ClusterChargeOverPath /
APV->PreviousGain;
1370 return ClusterChargeOverPath;
1374 const uint32_t detid,
1382 if (dynamic_cast<const StripGeomDetUnit*>(it) ==
nullptr && dynamic_cast<const PixelGeomDetUnit*>(it) ==
nullptr) {
1383 std::cout <<
"this detID doesn't seem to belong to the Tracker" << std::endl;
1388 const TrapezoidalPlaneBounds* trapezoidalBounds(dynamic_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
1389 const RectangularPlaneBounds* rectangularBounds(dynamic_cast<const RectangularPlaneBounds*>(&(plane.bounds())));
1391 double DistFromBorder = 1.0;
1395 if (trapezoidalBounds) {
1396 std::array<const float, 4>
const&
parameters = (*trapezoidalBounds).parameters();
1400 }
else if (rectangularBounds) {
1408 if (fabs(HitLocalPos.
y()) + HitLocalError.
yy() >= (
HalfLength - DistFromBorder))
1416 double adcs_err = 0.;
1417 double width = -0.5;
1418 double width_err = 0;
1419 double chi2overndf = -0.5;
1421 double nr_of_entries = InputHisto->GetEntries();
1424 FitResults[0] = adcs;
1425 FitResults[1] = adcs_err;
1426 FitResults[2] =
width;
1427 FitResults[3] = width_err;
1428 FitResults[4] = chi2overndf;
1433 TF1* MyLandau =
new TF1(
"MyLandau",
"landau", LowRange, HighRange);
1434 MyLandau->SetParameter(
"MPV", 300);
1436 InputHisto->Fit(MyLandau,
"QR WW");
1439 adcs = MyLandau->GetParameter(
"MPV");
1440 adcs_err = MyLandau->GetParError(1);
1441 width = MyLandau->GetParameter(2);
1442 width_err = MyLandau->GetParError(2);
1443 chi2overndf = MyLandau->GetChisquare() / MyLandau->GetNDF();
1454 FitResults[0] = adcs;
1455 FitResults[1] = adcs_err;
1456 FitResults[2] =
width;
1457 FitResults[3] = width_err;
1458 FitResults[4] = chi2overndf;
1464 cout <<
"START getNewObject\n";
1467 if (!(strcmp(
AlgoMode.c_str(),
"WriteOnDB") == 0 || strcmp(
AlgoMode.c_str(),
"SingleJob") == 0))
1468 return std::make_unique<SiStripApvGain>();
1470 auto obj = std::make_unique<SiStripApvGain>();
1471 std::vector<float>* theSiStripVector =
nullptr;
1472 unsigned int PreviousDetId = -1;
1475 if (
APV ==
nullptr) {
1479 if (
APV->DetId != PreviousDetId) {
1480 if (theSiStripVector !=
nullptr) {
1482 if (!
obj->put(PreviousDetId,
range))
1483 printf(
"Bug to put detId = %i\n", PreviousDetId);
1486 theSiStripVector =
new std::vector<float>;
1487 PreviousDetId =
APV->DetId;
1489 printf(
"%i | %i | PreviousGain = %7.5f NewGain = %7.5f\n",
APV->DetId,
APV->APVId,
APV->PreviousGain,
APV->Gain);
1490 if (theSiStripVector !=
nullptr) {
1491 theSiStripVector->push_back(
APV->Gain);
1496 if (theSiStripVector !=
nullptr) {
1498 if (!
obj->put(PreviousDetId,
range))
1499 printf(
"Bug to put detId = %i\n", PreviousDetId);
1502 cout <<
"END getNewObject\n";
edm::ESGetToken< SiStripGain, SiStripGainRcd > gainToken_
static constexpr auto TEC
virtual int nstrips() const =0
T getParameter(std::string const &) const
void algoEndJob() override
uint16_t firstStrip() const
virtual float length() const =0
SiStripGainFromData(const edm::ParameterSet &)
T const & getData(const ESGetToken< T, R > &iToken) const noexcept(false)
TH2F * HitLocalPositionBefCut
edm::ESGetToken< TrackerGeometry, TrackerDigiGeometryRecord > tkGeomToken_
const LocalTrajectoryError & localError() const
LocalPoint localPosition() const
friend struct const_iterator
edm::ESGetToken< TrackerTopology, TrackerTopologyRcd > tTopoToken_
unsigned int MinNrEntries
std::vector< stAPVGain * > APVsCollOrdered
unsigned int tidSide(const DetId &id) const
TH1F * NSatStripInCluster
bool IsFarFromBorder(TrajectoryStateOnSurface trajState, const uint32_t detid, const edm::EventSetup *iSetup)
T const * product() const
dqm::legacy::MonitorElement MonitorElement
void algoBeginRun(const edm::Run &, const edm::EventSetup &) override
bool operator()(const T &PseudoDetId1, const T &PseudoDetId2)
LocalError positionError() const
TH1F * FWHM_Vs_PathLength320
const edm::EventSetup * iSetup_
ClusterRef cluster() const
dqm::legacy::DQMStore DQMStore
TH2F * NStrips_Vs_TransversAngle
TH2F * Charge_Vs_PathTEC1
TH1F * FWHM_Vs_PathLength500
DataContainer const & measurements() const
const_iterator end() const
last iterator over the map (read only)
T getUntrackedParameter(std::string const &, T const &) const
TH1F * MPV_Vs_TransversAngle
virtual TH2F * getTH2F() const
vector< string > VInputFiles
TH1F * FWHM_Vs_PathLength
std::unordered_map< unsigned int, stAPVGain * > APVsColl
SiStripCluster const & monoCluster() const
unsigned int tecSide(const DetId &id) const
SiStripCluster const & amplitudes() const
ConstRecHitContainer recHits() const
LocalVector localDirection() const
std::pair< ContainerIterator, ContainerIterator > Range
const std::complex< double > I
TH1F * MPV_Vs_PathLength320
#define DEFINE_FWK_MODULE(type)
TH2F * Charge_Vs_PathTEC2
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector's numbering enum) ...
static constexpr auto TOB
double MaxTrackChiOverNdf
void algoBeginJob(const edm::EventSetup &) override
TH1F * NumberOfEntriesByAPV
unsigned long long TimeValue_t
ESHandle< T > getHandle(const ESGetToken< T, R > &iToken) const
TH2F * Charge_Vs_PathLength500
unsigned int MinTrackHits
ConstRecHitContainer RecHitContainer
static constexpr auto TIB
const Plane & surface() const
The nominal surface of the GeomDet.
TH2F * Charge_Vs_PathLength
double ComputeChargeOverPath(const SiStripCluster *Cluster, TrajectoryStateOnSurface trajState, const edm::EventSetup *iSetup, const Track *track, double trajChi2OverN)
void algoAnalyze(const edm::Event &, const edm::EventSetup &) override
constexpr uint32_t rawId() const
get the raw id
TH2F * Charge_Vs_TransversAngle
virtual TH1F * getTH1F() const
MonitorElement * book2D(TString const &name, TString const &title, int nchX, double lowX, double highX, int nchY, double lowY, double highY, FUNC onbooking=NOOP())
std::unique_ptr< SiStripApvGain > getNewObject() override
GlobalVector globalMomentum() const
SiStripCluster const & stereoCluster() const
void getPeakOfLandau(TH1 *InputHisto, double *FitResults, double LowRange=0, double HighRange=5400)
std::string TrajToTrackLabel
DQM_DEPRECATED void save(std::string const &filename, std::string const &path="")
const_iterator begin() const
first iterator over the map (read only)
TH2F * Charge_Vs_PathLength320
TH1F * NHighStripInCluster
MonitorElement * book1D(TString const &name, TString const &title, int const nchX, double const lowX, double const highX, FUNC onbooking=NOOP())
std::string TrajToTrackProducer
ClusterRef cluster() const
static constexpr auto TID
const edm::Event * iEvent_
TH1F * MPV_Vs_PathLength500
const Bounds & bounds() const