61 #include "TObjString.h" 69 #include <unordered_map> 88 void algoEndJob()
override ;
100 void getPeakOfLandau(TH1* InputHisto,
double* FitResults,
double LowRange=0,
double HighRange=5400);
280 template <
class T>
bool operator () (
const T& PseudoDetId1,
const T& PseudoDetId2) {
return PseudoDetId1==PseudoDetId2; }
284 std::unordered_map<unsigned int, stAPVGain*>
APVsColl;
316 if( strcmp(
AlgoMode.c_str(),
"WriteOnDB")==0 )
411 if( strcmp(
AlgoMode.c_str(),
"MultiJob")!=0 ){
487 auto const &
Det = tkGeom->
dets();
493 if(!gainHandle.
isValid()){printf(
"\n#####################\n\nERROR --> gainHandle is not valid\n\n#####################\n\n");
exit(0);}
497 for(
unsigned int i=0;
i<
Det.size();
i++){
505 if(!DetUnit)
continue;
508 unsigned int NAPV = Topo.
nstrips()/128;
510 double Phi = DetUnit->position().basicVector().phi();
511 double Eta = DetUnit->position().basicVector().eta();
512 double R = DetUnit->position().basicVector().transverse();
513 double Thick = DetUnit->surface().bounds().thickness();
515 for(
unsigned int j=0;j<NAPV;j++){
567 if( strcmp(
AlgoMode.c_str(),
"WriteOnDB")==0 || strcmp(
AlgoMode.c_str(),
"Merge")==0){
568 TH1::AddDirectory(kTRUE);
570 TFile*
file =
nullptr;
572 printf(
"Loading New Input File : %s\n",
VInputFiles[
f].c_str());
574 FILE* doesFileExist = fopen(
VInputFiles[
f].c_str(),
"r" );
576 printf(
"File %s doesn't exist\n",
VInputFiles[
f].c_str());
579 fclose(doesFileExist);
582 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;}
583 APV_Charge ->Add( (TH1*) file->FindObjectAny(
"APV_Charge") , 1);
584 APV_Momentum ->Add( (TH1*) file->FindObjectAny(
"APV_Momentum") , 1);
585 APV_PathLength ->Add( (TH1*) file->FindObjectAny(
"APV_PathLength") , 1);
587 Tracks_P_Vs_Eta ->Add( (TH1*) file->FindObjectAny(
"Tracks_P_Vs_Eta") , 1);
588 Tracks_Pt_Vs_Eta ->Add( (TH1*) file->FindObjectAny(
"Tracks_Pt_Vs_Eta") , 1);
598 HTrackHits ->Add( (TH1*) file->FindObjectAny(
"TrackHits") , 1);
607 Charge_Vs_Alpha ->Add( (TH1*) file->FindObjectAny(
"Charge_Vs_Alpha") , 1);
608 NStrips_Vs_Alpha ->Add( (TH1*) file->FindObjectAny(
"NStrips_Vs_Alpha") , 1);
609 HFirstStrip ->Add( (TH1*) file->FindObjectAny(
"FirstStrip") , 1);
611 TH1F* JobInfo_tmp = (TH1F*) file->FindObjectAny(
"JobInfo");
612 NEvent += (
unsigned int) JobInfo_tmp->GetBinContent(JobInfo_tmp->GetXaxis()->FindBin(1));
613 unsigned int tmp_SRun = (
unsigned int) JobInfo_tmp->GetBinContent(JobInfo_tmp->GetXaxis()->FindBin(3));
614 unsigned int tmp_SEvent = (
unsigned int) JobInfo_tmp->GetBinContent(JobInfo_tmp->GetXaxis()->FindBin(4));
615 unsigned int tmp_ERun = (
unsigned int) JobInfo_tmp->GetBinContent(JobInfo_tmp->GetXaxis()->FindBin(6));
616 unsigned int tmp_EEvent = (
unsigned int) JobInfo_tmp->GetBinContent(JobInfo_tmp->GetXaxis()->FindBin(7));
626 printf(
"Deleting Current Input File\n");
639 if( strcmp(
AlgoMode.c_str(),
"MultiJob")!=0 ){
640 TH1D* Proj =
nullptr;
641 double* FitResults =
new double[5];
644 if( I%3650==0 ) printf(
"Fitting Histograms \t %6.2f%%\n",(100.0*I)/
APVsColl.size());
649 Proj =
APV_Charge->ProjectionY(
" ",bin,bin,
"e");
650 Proj = (TH1D*)Proj->Clone();
651 if(Proj==
nullptr)
continue;
655 int SecondAPVId = APV->
APVId;
656 if(SecondAPVId%2==0){
657 SecondAPVId = SecondAPVId+1;
659 SecondAPVId = SecondAPVId-1;
664 TH1D* Proj2 =
APV_Charge->ProjectionY(
" ",bin2,bin2,
"e");
677 TH1D* Proj2 =
APV_Charge->ProjectionY(
" ",bin2,bin2,
"e");
690 APV->
MPV = FitResults[0];
692 if(FitResults[0]!=-0.5 && FitResults[1]<
MaxMPVError){
743 if(FitResults[0]!=-0.5){
755 if(Proj==
nullptr)
continue;
762 unsigned int GOOD = 0;
763 unsigned int BAD = 0;
764 double MPVmean =
MPVs->GetMean();
770 APV->
Gain = APV->
MPV / MPVmean;
940 fprintf(Gains,
"NEvents = %i\n",
NEvent);
941 fprintf(Gains,
"Number of APVs = %lu\n",static_cast<unsigned long>(
APVsColl.size()));
942 fprintf(Gains,
"GoodFits = %i BadFits = %i ratio = %f\n",GOOD,BAD,(100.0*GOOD)/(GOOD+BAD));
948 std::vector<int> DetIdOfBuggedAPV;
949 fprintf(Gains,
"----------------------------------------------------------------------\n");
952 if(APV->
MPV>0 && APV->
MPV<200){
954 for(
unsigned int b=0;
b<DetIdOfBuggedAPV.size()&&!tmpBug;
b++){
if(DetIdOfBuggedAPV[
b]==APV->
DetId)tmpBug=
true;}
955 if(!tmpBug){fprintf(Gains,
"%i,\n",APV->
DetId);DetIdOfBuggedAPV.push_back(APV->
DetId);}
977 if(!gainHandle.
isValid()){printf(
"\n#####################\n\nERROR --> gainHandle is not valid\n\n#####################\n\n");
exit(0);}
990 printf(
"GAIN HANDLE IS NOT VALID\n");
1004 if( strcmp(
AlgoMode.c_str(),
"WriteOnDB")==0 )
return;
1037 for(Trajectory::RecHitContainer::const_iterator rechit = transRecHits.begin(); rechit != transRecHits.end(); ++rechit)
1038 if ((*rechit)->isValid()) ndof += (*rechit)->dimension();
1045 vector<TrajectoryMeasurement> measurements = traj.
measurements();
1068 for(vector<TrajectoryMeasurement>::const_iterator measurement_it = measurements.begin(); measurement_it!=measurements.end(); measurement_it++){
1071 if( !trajState.
isValid() )
continue;
1078 if(sistripsimplehit){
1080 }
else if(sistripmatchedhit){
1083 }
else if(sistripsimple1dhit){
1097 double cosine = trackDirection.
z()/trackDirection.
mag();
1101 int APVId = FirstStrip/128;
1104 bool Overlaping =
false;
1106 unsigned int NHighStrip = 0;
1111 if(FirstStrip==0 )Overlaping=
true;
1112 if(FirstStrip==128 )Overlaping=
true;
1113 if(FirstStrip==256 )Overlaping=
true;
1114 if(FirstStrip==384 )Overlaping=
true;
1115 if(FirstStrip==512 )Overlaping=
true;
1116 if(FirstStrip==640 )Overlaping=
true;
1118 if(FirstStrip<=127 && FirstStrip+Ampls.size()>127)Overlaping=
true;
1119 if(FirstStrip<=255 && FirstStrip+Ampls.size()>255)Overlaping=
true;
1120 if(FirstStrip<=383 && FirstStrip+Ampls.size()>383)Overlaping=
true;
1121 if(FirstStrip<=511 && FirstStrip+Ampls.size()>511)Overlaping=
true;
1122 if(FirstStrip<=639 && FirstStrip+Ampls.size()>639)Overlaping=
true;
1124 if(FirstStrip+Ampls.size()==127 )Overlaping=
true;
1125 if(FirstStrip+Ampls.size()==255 )Overlaping=
true;
1126 if(FirstStrip+Ampls.size()==383 )Overlaping=
true;
1127 if(FirstStrip+Ampls.size()==511 )Overlaping=
true;
1128 if(FirstStrip+Ampls.size()==639 )Overlaping=
true;
1129 if(FirstStrip+Ampls.size()==767 )Overlaping=
true;
1130 if(Overlaping)
return -1;
1142 for(
unsigned int a=0;
a<Ampls.size();
a++){Charge+=Ampls[
a];
if(Ampls[
a]>=254)Saturation++;
if(Ampls[
a]>=20)NHighStrip++;}
1144 double ClusterChargeOverPath = (double)Charge / path ;
1160 double trans = atan2(trackDirection.
y(),trackDirection.
x())*(180/3.14159265);
1161 double alpha = acos(trackDirection.
x() /
sqrt(
pow(trackDirection.
x(),2) +
pow(trackDirection.
z(),2) ) ) * (180/3.14159265);
1162 double beta = acos(trackDirection.
y() /
sqrt(
pow(trackDirection.
x(),2) +
pow(trackDirection.
z(),2) ) ) * (180/3.14159265);
1164 if(path>0.4 && path<0.45){
1190 return ClusterChargeOverPath;
1201 if (dynamic_cast<const StripGeomDetUnit*>(it)==
nullptr && dynamic_cast<const PixelGeomDetUnit*>(it)==
nullptr) {
1202 std::cout <<
"this detID doesn't seem to belong to the Tracker" << std::endl;
1207 const TrapezoidalPlaneBounds* trapezoidalBounds( dynamic_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
1208 const RectangularPlaneBounds* rectangularBounds( dynamic_cast<const RectangularPlaneBounds*>(&(plane.bounds())));
1210 double DistFromBorder = 1.0;
1214 if(trapezoidalBounds)
1216 std::array<const float, 4>
const &
parameters = (*trapezoidalBounds).parameters();
1217 HalfLength = parameters[3];
1220 }
else if(rectangularBounds){
1223 }
else{
return false;}
1226 if (fabs(HitLocalPos.
y())+HitLocalError.
yy() >= (HalfLength - DistFromBorder) )
return false;
1234 double adcs_err = 0.;
1235 double width = -0.5;
1236 double width_err = 0;
1237 double chi2overndf = -0.5;
1239 double nr_of_entries = InputHisto->GetEntries();
1242 FitResults[0] = adcs;
1243 FitResults[1] = adcs_err;
1244 FitResults[2] =
width;
1245 FitResults[3] = width_err;
1246 FitResults[4] = chi2overndf;
1251 TF1* MyLandau =
new TF1(
"MyLandau",
"landau",LowRange, HighRange);
1252 MyLandau->SetParameter(
"MPV",300);
1254 InputHisto->Fit(MyLandau,
"QR WW");
1257 adcs = MyLandau->GetParameter(
"MPV");
1258 adcs_err = MyLandau->GetParError(1);
1259 width = MyLandau->GetParameter(2);
1260 width_err = MyLandau->GetParError(2);
1261 chi2overndf = MyLandau->GetChisquare() / MyLandau->GetNDF();
1265 adcs = -0.5; adcs_err = 0.;
1266 width = -0.5; width_err = 0;
1270 FitResults[0] = adcs;
1271 FitResults[1] = adcs_err;
1272 FitResults[2] =
width;
1273 FitResults[3] = width_err;
1274 FitResults[4] = chi2overndf;
1284 cout <<
"START getNewObject\n";
1291 std::vector<float>* theSiStripVector =
nullptr;
1292 int PreviousDetId = -1;
1296 if(APV==
nullptr){ printf(
"Bug\n");
continue; }
1297 if(APV->
DetId != PreviousDetId){
1298 if(theSiStripVector!=
nullptr){
1300 if ( !obj->
put(PreviousDetId,range) ) printf(
"Bug to put detId = %i\n",PreviousDetId);
1303 theSiStripVector =
new std::vector<float>;
1304 PreviousDetId = APV->
DetId;
1307 theSiStripVector->push_back(APV->
Gain);
1311 if(theSiStripVector!=
nullptr){
1313 if ( !obj->
put(PreviousDetId,range) ) printf(
"Bug to put detId = %i\n",PreviousDetId);
1316 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
#define DEFINE_FWK_MODULE(type)
SiStripCluster const & monoCluster() const
LocalVector localDirection() const
LocalPoint localPosition() const
TH1F * NSatStripInCluster
bool IsFarFromBorder(TrajectoryStateOnSurface trajState, const uint32_t detid, const edm::EventSetup *iSetup)
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
bool put(const uint32_t &detID, Range input)
uint32_t rawId() const
get the raw id
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
~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
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
int subdetId() const
get the contents of the subdetector field (not cast into any detector's numbering enum) ...
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
TH1F * getTH1F(void) const
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)
TH2F * getTH2F(void) const
std::string TrajToTrackProducer
T const * product() const
SiStripApvGain * getNewObject() override
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_
DQMStore * dqmStore_infile
unsigned int tecSide(const DetId &id) const
TH1F * MPV_Vs_PathLength500
const SiStripApvGain::Range getRange(uint32_t detID) const