60 #include "TObjString.h" 68 #include <unordered_map> 88 void algoEndJob()
override ;
93 std::unique_ptr<SiStripApvGain> getNewObject()
override;
99 void getPeakOfLandau(TH1* InputHisto,
double* FitResults,
double LowRange=0,
double HighRange=5400);
279 template <
class T>
bool operator () (
const T& PseudoDetId1,
const T& PseudoDetId2) {
return PseudoDetId1==PseudoDetId2; }
283 std::unordered_map<unsigned int, stAPVGain*>
APVsColl;
315 if( strcmp(
AlgoMode.c_str(),
"WriteOnDB")==0 )
409 if( strcmp(
AlgoMode.c_str(),
"MultiJob")!=0 ){
485 auto const &
Det = tkGeom->
dets();
491 if(!gainHandle.
isValid()){printf(
"\n#####################\n\nERROR --> gainHandle is not valid\n\n#####################\n\n");
exit(0);}
495 for(
unsigned int i=0;
i<
Det.size();
i++){
503 if(!DetUnit)
continue;
506 unsigned int NAPV = Topo.
nstrips()/128;
508 double Phi = DetUnit->position().basicVector().phi();
509 double Eta = DetUnit->position().basicVector().eta();
510 double R = DetUnit->position().basicVector().transverse();
511 double Thick = DetUnit->surface().bounds().thickness();
513 for(
unsigned int j=0;j<NAPV;j++){
565 if( strcmp(
AlgoMode.c_str(),
"WriteOnDB")==0 || strcmp(
AlgoMode.c_str(),
"Merge")==0){
566 TH1::AddDirectory(kTRUE);
568 TFile*
file =
nullptr;
570 printf(
"Loading New Input File : %s\n",
VInputFiles[
f].c_str());
572 FILE* doesFileExist = fopen(
VInputFiles[
f].c_str(),
"r" );
574 printf(
"File %s doesn't exist\n",
VInputFiles[
f].c_str());
577 fclose(doesFileExist);
580 file =
new TFile(
VInputFiles[
f].c_str() );
if(!file || file->IsZombie() ){printf(
"### Bug With File %s\n### File will be skipped \n",
VInputFiles[
f].c_str());
continue;}
581 APV_Charge ->Add( (TH1*) file->FindObjectAny(
"APV_Charge") , 1);
582 APV_Momentum ->Add( (TH1*) file->FindObjectAny(
"APV_Momentum") , 1);
583 APV_PathLength ->Add( (TH1*) file->FindObjectAny(
"APV_PathLength") , 1);
585 Tracks_P_Vs_Eta ->Add( (TH1*) file->FindObjectAny(
"Tracks_P_Vs_Eta") , 1);
586 Tracks_Pt_Vs_Eta ->Add( (TH1*) file->FindObjectAny(
"Tracks_Pt_Vs_Eta") , 1);
596 HTrackHits ->Add( (TH1*) file->FindObjectAny(
"TrackHits") , 1);
605 Charge_Vs_Alpha ->Add( (TH1*) file->FindObjectAny(
"Charge_Vs_Alpha") , 1);
606 NStrips_Vs_Alpha ->Add( (TH1*) file->FindObjectAny(
"NStrips_Vs_Alpha") , 1);
607 HFirstStrip ->Add( (TH1*) file->FindObjectAny(
"FirstStrip") , 1);
609 TH1F* JobInfo_tmp = (TH1F*) file->FindObjectAny(
"JobInfo");
610 NEvent += (
unsigned int) JobInfo_tmp->GetBinContent(JobInfo_tmp->GetXaxis()->FindBin(1));
611 unsigned int tmp_SRun = (
unsigned int) JobInfo_tmp->GetBinContent(JobInfo_tmp->GetXaxis()->FindBin(3));
612 unsigned int tmp_SEvent = (
unsigned int) JobInfo_tmp->GetBinContent(JobInfo_tmp->GetXaxis()->FindBin(4));
613 unsigned int tmp_ERun = (
unsigned int) JobInfo_tmp->GetBinContent(JobInfo_tmp->GetXaxis()->FindBin(6));
614 unsigned int tmp_EEvent = (
unsigned int) JobInfo_tmp->GetBinContent(JobInfo_tmp->GetXaxis()->FindBin(7));
624 printf(
"Deleting Current Input File\n");
637 if( strcmp(
AlgoMode.c_str(),
"MultiJob")!=0 ){
638 TH1D* Proj =
nullptr;
639 double FitResults[5];
642 if( I%3650==0 ) printf(
"Fitting Histograms \t %6.2f%%\n",(100.0*I)/
APVsColl.size());
647 Proj =
APV_Charge->ProjectionY(
" ",bin,bin,
"e");
648 Proj = (TH1D*)Proj->Clone();
649 if(Proj==
nullptr)
continue;
653 int SecondAPVId = APV->
APVId;
654 if(SecondAPVId%2==0){
655 SecondAPVId = SecondAPVId+1;
657 SecondAPVId = SecondAPVId-1;
662 TH1D* Proj2 =
APV_Charge->ProjectionY(
" ",bin2,bin2,
"e");
675 TH1D* Proj2 =
APV_Charge->ProjectionY(
" ",bin2,bin2,
"e");
688 APV->
MPV = FitResults[0];
690 if(FitResults[0]!=-0.5 && FitResults[1]<
MaxMPVError){
741 if(FitResults[0]!=-0.5){
753 if(Proj==
nullptr)
continue;
760 unsigned int GOOD = 0;
761 unsigned int BAD = 0;
762 double MPVmean =
MPVs->GetMean();
768 APV->
Gain = APV->
MPV / MPVmean;
938 fprintf(Gains,
"NEvents = %i\n",
NEvent);
939 fprintf(Gains,
"Number of APVs = %lu\n",static_cast<unsigned long>(
APVsColl.size()));
940 fprintf(Gains,
"GoodFits = %i BadFits = %i ratio = %f\n",GOOD,BAD,(100.0*GOOD)/(GOOD+BAD));
946 std::vector<int> DetIdOfBuggedAPV;
947 fprintf(Gains,
"----------------------------------------------------------------------\n");
950 if(APV->
MPV>0 && APV->
MPV<200){
952 for(
unsigned int b=0;
b<DetIdOfBuggedAPV.size()&&!tmpBug;
b++){
if(DetIdOfBuggedAPV[
b]==APV->
DetId)tmpBug=
true;}
953 if(!tmpBug){fprintf(Gains,
"%i,\n",APV->
DetId);DetIdOfBuggedAPV.push_back(APV->
DetId);}
975 if(!gainHandle.
isValid()){printf(
"\n#####################\n\nERROR --> gainHandle is not valid\n\n#####################\n\n");
exit(0);}
988 printf(
"GAIN HANDLE IS NOT VALID\n");
1002 if( strcmp(
AlgoMode.c_str(),
"WriteOnDB")==0 )
return;
1035 for(Trajectory::RecHitContainer::const_iterator rechit = transRecHits.begin(); rechit != transRecHits.end(); ++rechit)
1036 if ((*rechit)->isValid()) ndof += (*rechit)->dimension();
1043 vector<TrajectoryMeasurement> measurements = traj.
measurements();
1066 for(vector<TrajectoryMeasurement>::const_iterator measurement_it = measurements.begin(); measurement_it!=measurements.end(); measurement_it++){
1069 if( !trajState.
isValid() )
continue;
1076 if(sistripsimplehit){
1078 }
else if(sistripmatchedhit){
1081 }
else if(sistripsimple1dhit){
1095 double cosine = trackDirection.
z()/trackDirection.
mag();
1099 int APVId = FirstStrip/128;
1102 bool Overlaping =
false;
1104 unsigned int NHighStrip = 0;
1109 if(FirstStrip==0 )Overlaping=
true;
1110 if(FirstStrip==128 )Overlaping=
true;
1111 if(FirstStrip==256 )Overlaping=
true;
1112 if(FirstStrip==384 )Overlaping=
true;
1113 if(FirstStrip==512 )Overlaping=
true;
1114 if(FirstStrip==640 )Overlaping=
true;
1116 if(FirstStrip<=127 && FirstStrip+Ampls.size()>127)Overlaping=
true;
1117 if(FirstStrip<=255 && FirstStrip+Ampls.size()>255)Overlaping=
true;
1118 if(FirstStrip<=383 && FirstStrip+Ampls.size()>383)Overlaping=
true;
1119 if(FirstStrip<=511 && FirstStrip+Ampls.size()>511)Overlaping=
true;
1120 if(FirstStrip<=639 && FirstStrip+Ampls.size()>639)Overlaping=
true;
1122 if(FirstStrip+Ampls.size()==127 )Overlaping=
true;
1123 if(FirstStrip+Ampls.size()==255 )Overlaping=
true;
1124 if(FirstStrip+Ampls.size()==383 )Overlaping=
true;
1125 if(FirstStrip+Ampls.size()==511 )Overlaping=
true;
1126 if(FirstStrip+Ampls.size()==639 )Overlaping=
true;
1127 if(FirstStrip+Ampls.size()==767 )Overlaping=
true;
1128 if(Overlaping)
return -1;
1140 for(
unsigned int a=0;
a<Ampls.size();
a++){Charge+=Ampls[
a];
if(Ampls[
a]>=254)Saturation++;
if(Ampls[
a]>=20)NHighStrip++;}
1142 double ClusterChargeOverPath = (double)Charge / path ;
1158 double trans = atan2(trackDirection.
y(),trackDirection.
x())*(180/3.14159265);
1159 double alpha = acos(trackDirection.
x() /
sqrt(
pow(trackDirection.
x(),2) +
pow(trackDirection.
z(),2) ) ) * (180/3.14159265);
1160 double beta = acos(trackDirection.
y() /
sqrt(
pow(trackDirection.
x(),2) +
pow(trackDirection.
z(),2) ) ) * (180/3.14159265);
1162 if(path>0.4 && path<0.45){
1188 return ClusterChargeOverPath;
1199 if (dynamic_cast<const StripGeomDetUnit*>(it)==
nullptr && dynamic_cast<const PixelGeomDetUnit*>(it)==
nullptr) {
1200 std::cout <<
"this detID doesn't seem to belong to the Tracker" << std::endl;
1205 const TrapezoidalPlaneBounds* trapezoidalBounds( dynamic_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
1206 const RectangularPlaneBounds* rectangularBounds( dynamic_cast<const RectangularPlaneBounds*>(&(plane.bounds())));
1208 double DistFromBorder = 1.0;
1212 if(trapezoidalBounds)
1214 std::array<const float, 4>
const &
parameters = (*trapezoidalBounds).parameters();
1215 HalfLength = parameters[3];
1218 }
else if(rectangularBounds){
1221 }
else{
return false;}
1224 if (fabs(HitLocalPos.
y())+HitLocalError.
yy() >= (HalfLength - DistFromBorder) )
return false;
1232 double adcs_err = 0.;
1233 double width = -0.5;
1234 double width_err = 0;
1235 double chi2overndf = -0.5;
1237 double nr_of_entries = InputHisto->GetEntries();
1240 FitResults[0] = adcs;
1241 FitResults[1] = adcs_err;
1242 FitResults[2] =
width;
1243 FitResults[3] = width_err;
1244 FitResults[4] = chi2overndf;
1249 TF1* MyLandau =
new TF1(
"MyLandau",
"landau",LowRange, HighRange);
1250 MyLandau->SetParameter(
"MPV",300);
1252 InputHisto->Fit(MyLandau,
"QR WW");
1255 adcs = MyLandau->GetParameter(
"MPV");
1256 adcs_err = MyLandau->GetParError(1);
1257 width = MyLandau->GetParameter(2);
1258 width_err = MyLandau->GetParError(2);
1259 chi2overndf = MyLandau->GetChisquare() / MyLandau->GetNDF();
1263 adcs = -0.5; adcs_err = 0.;
1264 width = -0.5; width_err = 0;
1268 FitResults[0] = adcs;
1269 FitResults[1] = adcs_err;
1270 FitResults[2] =
width;
1271 FitResults[3] = width_err;
1272 FitResults[4] = chi2overndf;
1282 cout <<
"START getNewObject\n";
1285 if( !(strcmp(
AlgoMode.c_str(),
"WriteOnDB")==0 || strcmp(
AlgoMode.c_str(),
"SingleJob")==0) )
1286 return std::make_unique<SiStripApvGain>();
1289 auto obj = std::make_unique<SiStripApvGain>();
1290 std::vector<float>* theSiStripVector =
nullptr;
1291 int PreviousDetId = -1;
1295 if(APV==
nullptr){ printf(
"Bug\n");
continue; }
1296 if(APV->
DetId != PreviousDetId){
1297 if(theSiStripVector!=
nullptr){
1299 if ( !
obj->put(PreviousDetId,range) ) printf(
"Bug to put detId = %i\n",PreviousDetId);
1302 theSiStripVector =
new std::vector<float>;
1303 PreviousDetId = APV->
DetId;
1306 theSiStripVector->push_back(APV->
Gain);
1310 if(theSiStripVector!=
nullptr){
1312 if ( !
obj->put(PreviousDetId,range) ) printf(
"Bug to put detId = %i\n",PreviousDetId);
1315 cout <<
"END getNewObject\n";
ClusterRef cluster() const
double p() const
momentum vector magnitude
T getParameter(std::string const &) const
EventNumber_t event() const
T getUntrackedParameter(std::string const &, T const &) const
virtual float length() const =0
void algoEndJob() override
SiStripGainFromData(const edm::ParameterSet &)
TH2F * HitLocalPositionBefCut
friend struct const_iterator
const_iterator end() const
last iterator over the map (read only)
unsigned int MinNrEntries
std::vector< stAPVGain * > APVsCollOrdered
SiStripCluster const & monoCluster() const
LocalVector localDirection() const
LocalPoint localPosition() const
TH1F * NSatStripInCluster
bool IsFarFromBorder(TrajectoryStateOnSurface trajState, const uint32_t detid, const edm::EventSetup *iSetup)
constexpr uint32_t rawId() const
get the raw id
void algoBeginRun(const edm::Run &, const edm::EventSetup &) override
ConstRecHitContainer recHits() const
const Bounds & bounds() const
uint16_t firstStrip() const
TH1F * FWHM_Vs_PathLength320
const Plane & surface() const
The nominal surface of the GeomDet.
const edm::EventSetup * iSetup_
LocalError positionError() const
TH2F * NStrips_Vs_TransversAngle
TH2F * Charge_Vs_PathTEC1
TH1F * FWHM_Vs_PathLength500
TH1F * MPV_Vs_TransversAngle
const DetContainer & dets() const override
Returm a vector of all GeomDet (including all GeomDetUnits)
const TrackerGeomDet * idToDetUnit(DetId) const override
Return the pointer to the GeomDetUnit corresponding to a given DetId.
DataContainer const & measurements() const
vector< string > VInputFiles
#define DEFINE_FWK_MODULE(type)
~SiStripGainFromData() override
double eta() const
pseudorapidity of momentum vector
TH1F * FWHM_Vs_PathLength
std::unordered_map< unsigned int, stAPVGain * > APVsColl
double pt() const
track transverse momentum
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
TH1F * MPV_Vs_PathLength320
ClusterRef cluster() const
TH2F * Charge_Vs_PathTEC2
const LocalTrajectoryError & localError() const
double MaxTrackChiOverNdf
void algoBeginJob(const edm::EventSetup &) override
TH1F * NumberOfEntriesByAPV
unsigned long long TimeValue_t
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
TH2F * Charge_Vs_PathLength500
unsigned int MinTrackHits
bin
set the eta bin as selection string.
ConstRecHitContainer RecHitContainer
TH2F * Charge_Vs_PathLength
T const * product() const
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
TH2F * Charge_Vs_TransversAngle
virtual int nstrips() const =0
std::unique_ptr< SiStripApvGain > getNewObject() override
void getPeakOfLandau(TH1 *InputHisto, double *FitResults, double LowRange=0, double HighRange=5400)
std::string TrajToTrackLabel
SiStripCluster const & stereoCluster() const
TH2F * Charge_Vs_PathLength320
GlobalVector globalMomentum() const
TH1F * NHighStripInCluster
const_iterator begin() const
first iterator over the map (read only)
std::string TrajToTrackProducer
T const * product() const
TimeValue_t value() const
const std::vector< uint8_t > & amplitudes() const
edm::Timestamp time() const
Power< A, B >::type pow(const A &a, const B &b)
const edm::Event * iEvent_
unsigned int tecSide(const DetId &id) const
TH1F * MPV_Vs_PathLength500
const SiStripApvGain::Range getRange(uint32_t detID) const