58 #include "TObjString.h"
66 #include <unordered_map>
97 void algoEndJob()
override;
102 std::unique_ptr<SiStripApvGain> getNewObject()
override;
109 double trajChi2OverN);
112 void getPeakOfLandau(TH1* InputHisto,
double* FitResults,
double LowRange = 0,
double HighRange = 5400);
290 return PseudoDetId1 == PseudoDetId2;
295 std::unordered_map<unsigned int, stAPVGain*>
APVsColl;
326 if (strcmp(
AlgoMode.c_str(),
"WriteOnDB") == 0)
364 tmp =
dqmStore_->
book2D(
"APV_PathLength",
"APV_PathLength", 72785, 0, 72784, 100, 0.2, 1.4);
381 tmp =
dqmStore_->
book2D(
"Tracks_Pt_Vs_Eta",
"Tracks_Pt_Vs_Eta", 30, 0, 3, 50, 0, 200);
384 tmp =
dqmStore_->
book2D(
"Charge_Vs_PathTIB",
"Charge_Vs_PathTIB", 100, 0.2, 1.4, 500, 0, 2000);
386 tmp =
dqmStore_->
book2D(
"Charge_Vs_PathTID",
"Charge_Vs_PathTID", 100, 0.2, 1.4, 500, 0, 2000);
388 tmp =
dqmStore_->
book2D(
"Charge_Vs_PathTOB",
"Charge_Vs_PathTOB", 100, 0.2, 1.4, 500, 0, 2000);
390 tmp =
dqmStore_->
book2D(
"Charge_Vs_PathTEC",
"Charge_Vs_PathTEC", 100, 0.2, 1.4, 500, 0, 2000);
392 tmp =
dqmStore_->
book2D(
"Charge_Vs_PathTEC1",
"Charge_Vs_PathTEC1", 100, 0.2, 1.4, 500, 0, 2000);
394 tmp =
dqmStore_->
book2D(
"Charge_Vs_PathTEC2",
"Charge_Vs_PathTEC2", 100, 0.2, 1.4, 500, 0, 2000);
425 tmp =
dqmStore_->
book2D(
"Charge_Vs_PathLength",
"Charge_Vs_PathLength", 100, 0.2, 1.4, 500, 0, 2000);
427 tmp =
dqmStore_->
book2D(
"Charge_Vs_PathLength320",
"Charge_Vs_PathLength", 100, 0.2, 1.4, 500, 0, 2000);
429 tmp =
dqmStore_->
book2D(
"Charge_Vs_PathLength500",
"Charge_Vs_PathLength", 100, 0.2, 1.4, 500, 0, 2000);
432 tmp =
dqmStore_->
book2D(
"Charge_Vs_TransversAngle",
"Charge_Vs_TransversAngle", 220, -20, 200, 500, 0, 2000);
434 tmp =
dqmStore_->
book2D(
"Charge_Vs_Alpha",
"Charge_Vs_Alpha", 220, -20, 200, 500, 0, 2000);
436 tmp =
dqmStore_->
book2D(
"Charge_Vs_Beta",
"Charge_Vs_Beta", 220, -20, 200, 500, 0, 2000);
439 tmp =
dqmStore_->
book2D(
"NStrips_Vs_TransversAngle",
"NStrips_Vs_TransversAngle", 220, -20, 200, 10, 0, 10);
441 tmp =
dqmStore_->
book2D(
"NStrips_Vs_Alpha",
"NStrips_Vs_Alpha", 220, -20, 200, 10, 0, 10);
443 tmp =
dqmStore_->
book2D(
"NStrips_Vs_Beta",
"NStrips_Vs_Beta", 220, -20, 200, 10, 0, 10);
458 if (strcmp(
AlgoMode.c_str(),
"MultiJob") != 0) {
459 tmp =
dqmStore_->
book2D(
"MPV_Vs_EtaTIB",
"MPV_Vs_EtaTIB", 50, -3.0, 3.0, 300, 0, 600);
461 tmp =
dqmStore_->
book2D(
"MPV_Vs_EtaTID",
"MPV_Vs_EtaTID", 50, -3.0, 3.0, 300, 0, 600);
463 tmp =
dqmStore_->
book2D(
"MPV_Vs_EtaTOB",
"MPV_Vs_EtaTOB", 50, -3.0, 3.0, 300, 0, 600);
465 tmp =
dqmStore_->
book2D(
"MPV_Vs_EtaTEC",
"MPV_Vs_EtaTEC", 50, -3.0, 3.0, 300, 0, 600);
467 tmp =
dqmStore_->
book2D(
"MPV_Vs_EtaTEC1",
"MPV_Vs_EtaTEC1", 50, -3.0, 3.0, 300, 0, 600);
469 tmp =
dqmStore_->
book2D(
"MPV_Vs_EtaTEC2",
"MPV_Vs_EtaTEC2", 50, -3.0, 3.0, 300, 0, 600);
472 tmp =
dqmStore_->
book2D(
"MPV_Vs_PhiTIB",
"MPV_Vs_PhiTIB", 50, -3.2, 3.2, 300, 0, 600);
474 tmp =
dqmStore_->
book2D(
"MPV_Vs_PhiTID",
"MPV_Vs_PhiTID", 50, -3.2, 3.2, 300, 0, 600);
476 tmp =
dqmStore_->
book2D(
"MPV_Vs_PhiTOB",
"MPV_Vs_PhiTOB", 50, -3.2, 3.2, 300, 0, 600);
478 tmp =
dqmStore_->
book2D(
"MPV_Vs_PhiTEC",
"MPV_Vs_PhiTEC", 50, -3.2, 3.2, 300, 0, 600);
480 tmp =
dqmStore_->
book2D(
"MPV_Vs_PhiTEC1",
"MPV_Vs_PhiTEC1", 50, -3.2, 3.2, 300, 0, 600);
482 tmp =
dqmStore_->
book2D(
"MPV_Vs_PhiTEC2",
"MPV_Vs_PhiTEC2", 50, -3.2, 3.2, 300, 0, 600);
526 tmp =
dqmStore_->
book1D(
"FWHM_Vs_PathLength320",
"FWHM_Vs_PathLength", 100, 0.2, 1.4);
528 tmp =
dqmStore_->
book1D(
"FWHM_Vs_PathLength500",
"FWHM_Vs_PathLength", 100, 0.2, 1.4);
531 tmp =
dqmStore_->
book1D(
"MPV_Vs_TransversAngle",
"MPV_Vs_TransversAngle", 220, -20, 200);
540 tmp =
dqmStore_->
book2D(
"Error_Vs_Entries",
"Error_Vs_Entries", 500, 0, 10000, 100, 0, 50);
547 tmp =
dqmStore_->
book2D(
"NoMPV_Vs_EtaPhi",
"NoMPV_Vs_EtaPhi", 50, -3.0, 3.0, 50, -3.2, 3.2);
550 tmp =
dqmStore_->
book1D(
"NumberOfEntriesByAPV",
"NumberOfEntriesByAPV", 1000, 0, 10000);
569 auto const&
Det = tkGeom->
dets();
575 printf(
"\n#####################\n\nERROR --> gainHandle is not valid\n\n#####################\n\n");
581 for (
unsigned int i = 0;
i <
Det.size();
i++) {
587 auto DetUnit = dynamic_cast<StripGeomDetUnit const*>(
Det[
i]);
592 unsigned int NAPV = Topo.
nstrips() / 128;
594 double Phi = DetUnit->position().basicVector().phi();
595 double Eta = DetUnit->position().basicVector().eta();
596 double R = DetUnit->position().basicVector().transverse();
597 double Thick = DetUnit->surface().bounds().thickness();
599 for (
unsigned int j = 0;
j < NAPV;
j++) {
607 APV->PreviousGain = 1;
611 APV->Thickness = Thick;
646 if (strcmp(
AlgoMode.c_str(),
"WriteOnDB") == 0 || strcmp(
AlgoMode.c_str(),
"Merge") == 0) {
647 TH1::AddDirectory(kTRUE);
649 TFile*
file =
nullptr;
651 printf(
"Loading New Input File : %s\n",
VInputFiles[
f].c_str());
653 FILE* doesFileExist = fopen(
VInputFiles[
f].c_str(),
"r");
654 if (!doesFileExist) {
655 printf(
"File %s doesn't exist\n",
VInputFiles[
f].c_str());
658 fclose(doesFileExist);
663 printf(
"### Bug With File %s\n### File will be skipped \n",
VInputFiles[
f].c_str());
694 TH1F* JobInfo_tmp = (TH1F*)
file->FindObjectAny(
"JobInfo");
695 NEvent += (
unsigned int)JobInfo_tmp->GetBinContent(JobInfo_tmp->GetXaxis()->FindBin(1));
696 unsigned int tmp_SRun = (
unsigned int)JobInfo_tmp->GetBinContent(JobInfo_tmp->GetXaxis()->FindBin(3));
697 unsigned int tmp_SEvent = (
unsigned int)JobInfo_tmp->GetBinContent(JobInfo_tmp->GetXaxis()->FindBin(4));
698 unsigned int tmp_ERun = (
unsigned int)JobInfo_tmp->GetBinContent(JobInfo_tmp->GetXaxis()->FindBin(6));
699 unsigned int tmp_EEvent = (
unsigned int)JobInfo_tmp->GetBinContent(JobInfo_tmp->GetXaxis()->FindBin(7));
704 if (tmp_SRun <
SRun) {
707 }
else if (tmp_SRun ==
SRun && tmp_SEvent <
SEvent) {
711 if (tmp_ERun >
ERun) {
714 }
else if (tmp_ERun ==
ERun && tmp_EEvent >
EEvent) {
718 printf(
"Deleting Current Input File\n");
730 if (strcmp(
AlgoMode.c_str(),
"MultiJob") != 0) {
731 TH1D* Proj =
nullptr;
732 double FitResults[5];
736 printf(
"Fitting Histograms \t %6.2f%%\n", (100.0 *
I) /
APVsColl.size());
742 Proj = (TH1D*)Proj->Clone();
748 int SecondAPVId =
APV->APVId;
749 if (SecondAPVId % 2 == 0) {
750 SecondAPVId = SecondAPVId + 1;
752 SecondAPVId = SecondAPVId - 1;
757 TH1D* Proj2 =
APV_Charge->ProjectionY(
" ", bin2, bin2,
"e");
758 if (Proj2 !=
nullptr) {
772 TH1D* Proj2 =
APV_Charge->ProjectionY(
" ", bin2, bin2,
"e");
773 if (Proj2 !=
nullptr) {
784 APV->MPV = FitResults[0];
786 if (FitResults[0] != -0.5 && FitResults[1] <
MaxMPVError) {
789 if (
APV->Thickness < 0.04)
791 if (
APV->Thickness > 0.04)
804 if (
APV->Thickness < 0.04)
806 if (
APV->Thickness > 0.04)
818 if (
APV->Thickness < 0.04)
820 if (
APV->Thickness > 0.04)
837 if (
APV->Thickness < 0.04)
839 if (
APV->Thickness > 0.04)
861 if (
APV->Thickness < 0.04)
863 if (
APV->Thickness > 0.04)
871 if (FitResults[0] != -0.5) {
890 unsigned int GOOD = 0;
891 unsigned int BAD = 0;
892 double MPVmean =
MPVs->GetMean();
897 APV->Gain =
APV->MPV / MPVmean;
909 APV->Gain *=
APV->PreviousGain;
916 if (FitResults[0] == -0.5)
930 if (FitResults[0] == -0.5)
944 if (FitResults[0] == -0.5)
1009 if (FitResults[0] == -0.5)
1019 if (FitResults[0] == -0.5)
1029 if (FitResults[0] == -0.5)
1039 if (FitResults[0] == -0.5)
1049 if (FitResults[0] == -0.5)
1059 if (FitResults[0] == -0.5)
1069 if (FitResults[0] == -0.5)
1079 if (FitResults[0] == -0.5)
1089 if (FitResults[0] == -0.5)
1097 fprintf(Gains,
"NEvents = %i\n",
NEvent);
1098 fprintf(Gains,
"Number of APVs = %lu\n", static_cast<unsigned long>(
APVsColl.size()));
1099 fprintf(Gains,
"GoodFits = %i BadFits = %i ratio = %f\n",
GOOD,
BAD, (100.0 *
GOOD) / (
GOOD +
BAD));
1103 "%i | %i | PreviousGain = %7.5f NewGain = %7.5f\n",
1110 std::vector<int> DetIdOfBuggedAPV;
1111 fprintf(Gains,
"----------------------------------------------------------------------\n");
1114 if (
APV->MPV > 0 &&
APV->MPV < 200) {
1115 bool tmpBug =
false;
1116 for (
unsigned int b = 0;
b < DetIdOfBuggedAPV.size() && !tmpBug;
b++) {
1117 if (DetIdOfBuggedAPV[
b] ==
APV->DetId)
1121 fprintf(Gains,
"%i,\n",
APV->DetId);
1122 DetIdOfBuggedAPV.push_back(
APV->DetId);
1142 printf(
"\n#####################\n\nERROR --> gainHandle is not valid\n\n#####################\n\n");
1151 APV->PreviousGain = *(detGainRange.first +
APV->APVId);
1154 if (
APV->PreviousGain < 0)
1155 APV->PreviousGain = 1;
1157 printf(
"GAIN HANDLE IS NOT VALID\n");
1164 if (strcmp(
AlgoMode.c_str(),
"WriteOnDB") == 0)
1198 for (Trajectory::RecHitContainer::const_iterator rechit = transRecHits.begin(); rechit != transRecHits.end();
1200 if ((*rechit)->isValid())
1201 ndof += (*rechit)->dimension();
1209 vector<TrajectoryMeasurement> measurements = traj.
measurements();
1233 for (vector<TrajectoryMeasurement>::const_iterator measurement_it = measurements.begin();
1234 measurement_it != measurements.end();
1241 const SiStripRecHit1D* sistripsimple1dhit = dynamic_cast<const SiStripRecHit1D*>(
hit);
1245 if (sistripsimplehit) {
1248 }
else if (sistripmatchedhit) {
1252 }
else if (sistripsimple1dhit) {
1265 double trajChi2OverN) {
1267 double cosine = trackDirection.
z() / trackDirection.
mag();
1271 int APVId = FirstStrip / 128;
1274 bool Overlaping =
false;
1276 unsigned int NHighStrip = 0;
1281 if (FirstStrip == 0)
1283 if (FirstStrip == 128)
1285 if (FirstStrip == 256)
1287 if (FirstStrip == 384)
1289 if (FirstStrip == 512)
1291 if (FirstStrip == 640)
1294 if (FirstStrip <= 127 && FirstStrip + Ampls.size() > 127)
1296 if (FirstStrip <= 255 && FirstStrip + Ampls.size() > 255)
1298 if (FirstStrip <= 383 && FirstStrip + Ampls.size() > 383)
1300 if (FirstStrip <= 511 && FirstStrip + Ampls.size() > 511)
1302 if (FirstStrip <= 639 && FirstStrip + Ampls.size() > 639)
1305 if (FirstStrip + Ampls.size() == 127)
1307 if (FirstStrip + Ampls.size() == 255)
1309 if (FirstStrip + Ampls.size() == 383)
1311 if (FirstStrip + Ampls.size() == 511)
1313 if (FirstStrip + Ampls.size() == 639)
1315 if (FirstStrip + Ampls.size() == 767)
1329 for (
unsigned int a = 0;
a < Ampls.size();
a++) {
1331 if (Ampls[
a] >= 254)
1336 double path = (10.0 *
APV->Thickness) / fabs(cosine);
1337 double ClusterChargeOverPath = (double)
Charge /
path;
1346 if (
APV->Thickness < 0.04)
1348 if (
APV->Thickness > 0.04)
1358 if (
APV->Thickness < 0.04)
1360 if (
APV->Thickness > 0.04)
1364 double trans = atan2(trackDirection.
y(), trackDirection.
x()) * (180 / 3.14159265);
1366 acos(trackDirection.
x() /
sqrt(
pow(trackDirection.
x(), 2) +
pow(trackDirection.
z(), 2))) * (180 / 3.14159265);
1368 acos(trackDirection.
y() /
sqrt(
pow(trackDirection.
x(), 2) +
pow(trackDirection.
z(), 2))) * (180 / 3.14159265);
1390 ClusterChargeOverPath = ClusterChargeOverPath /
APV->PreviousGain;
1397 return ClusterChargeOverPath;
1401 const uint32_t detid,
1410 if (dynamic_cast<const StripGeomDetUnit*>(it) ==
nullptr && dynamic_cast<const PixelGeomDetUnit*>(it) ==
nullptr) {
1411 std::cout <<
"this detID doesn't seem to belong to the Tracker" << std::endl;
1416 const TrapezoidalPlaneBounds* trapezoidalBounds(dynamic_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
1417 const RectangularPlaneBounds* rectangularBounds(dynamic_cast<const RectangularPlaneBounds*>(&(plane.bounds())));
1419 double DistFromBorder = 1.0;
1423 if (trapezoidalBounds) {
1424 std::array<const float, 4>
const&
parameters = (*trapezoidalBounds).parameters();
1428 }
else if (rectangularBounds) {
1436 if (fabs(HitLocalPos.
y()) + HitLocalError.
yy() >= (
HalfLength - DistFromBorder))
1444 double adcs_err = 0.;
1445 double width = -0.5;
1446 double width_err = 0;
1447 double chi2overndf = -0.5;
1449 double nr_of_entries = InputHisto->GetEntries();
1452 FitResults[0] = adcs;
1453 FitResults[1] = adcs_err;
1454 FitResults[2] =
width;
1455 FitResults[3] = width_err;
1456 FitResults[4] = chi2overndf;
1461 TF1* MyLandau =
new TF1(
"MyLandau",
"landau", LowRange, HighRange);
1462 MyLandau->SetParameter(
"MPV", 300);
1464 InputHisto->Fit(MyLandau,
"QR WW");
1467 adcs = MyLandau->GetParameter(
"MPV");
1468 adcs_err = MyLandau->GetParError(1);
1469 width = MyLandau->GetParameter(2);
1470 width_err = MyLandau->GetParError(2);
1471 chi2overndf = MyLandau->GetChisquare() / MyLandau->GetNDF();
1482 FitResults[0] = adcs;
1483 FitResults[1] = adcs_err;
1484 FitResults[2] =
width;
1485 FitResults[3] = width_err;
1486 FitResults[4] = chi2overndf;
1492 cout <<
"START getNewObject\n";
1495 if (!(strcmp(
AlgoMode.c_str(),
"WriteOnDB") == 0 || strcmp(
AlgoMode.c_str(),
"SingleJob") == 0))
1496 return std::make_unique<SiStripApvGain>();
1498 auto obj = std::make_unique<SiStripApvGain>();
1499 std::vector<float>* theSiStripVector =
nullptr;
1500 int PreviousDetId = -1;
1503 if (
APV ==
nullptr) {
1507 if (
APV->DetId != PreviousDetId) {
1508 if (theSiStripVector !=
nullptr) {
1510 if (!
obj->put(PreviousDetId,
range))
1511 printf(
"Bug to put detId = %i\n", PreviousDetId);
1514 theSiStripVector =
new std::vector<float>;
1515 PreviousDetId =
APV->DetId;
1517 printf(
"%i | %i | PreviousGain = %7.5f NewGain = %7.5f\n",
APV->DetId,
APV->APVId,
APV->PreviousGain,
APV->Gain);
1518 theSiStripVector->push_back(
APV->Gain);
1522 if (theSiStripVector !=
nullptr) {
1524 if (!
obj->put(PreviousDetId,
range))
1525 printf(
"Bug to put detId = %i\n", PreviousDetId);
1528 cout <<
"END getNewObject\n";