61 #include "TObjString.h"
69 #include <ext/hash_map>
76 using namespace __gnu_cxx;
77 using __gnu_cxx::hash_map;
81 struct stAPVGain{
unsigned int Index;
int DetId;
int APVId;
int SubDet;
float Eta;
float R;
float Phi;
float Thickness;
double MPV;
double Gain;
double PreviousGain;
char Side;};
91 virtual void algoEndJob()
override ;
103 void getPeakOfLandau(TH1* InputHisto,
double* FitResults,
double LowRange=0,
double HighRange=5400);
283 template <
class T>
bool operator () (
const T& PseudoDetId1,
const T& PseudoDetId2) {
return PseudoDetId1==PseudoDetId2; }
287 __gnu_cxx::hash_map<unsigned int, stAPVGain*, __gnu_cxx::hash<unsigned int>,
isEqual >
APVsColl;
319 if( strcmp(
AlgoMode.c_str(),
"WriteOnDB")==0 )
414 if( strcmp(
AlgoMode.c_str(),
"MultiJob")!=0 ){
490 auto const &
Det = tkGeom->dets();
496 if(!gainHandle.
isValid()){printf(
"\n#####################\n\nERROR --> gainHandle is not valid\n\n#####################\n\n");
exit(0);}
500 for(
unsigned int i=0;
i<
Det.size();
i++){
508 if(!DetUnit)
continue;
511 unsigned int NAPV = Topo.
nstrips()/128;
513 double Phi = DetUnit->position().basicVector().phi();
514 double Eta = DetUnit->position().basicVector().eta();
515 double R = DetUnit->position().basicVector().transverse();
516 double Thick = DetUnit->surface().bounds().thickness();
518 for(
unsigned int j=0;
j<NAPV;
j++){
570 if( strcmp(
AlgoMode.c_str(),
"WriteOnDB")==0 || strcmp(
AlgoMode.c_str(),
"Merge")==0){
571 TH1::AddDirectory(kTRUE);
575 printf(
"Loading New Input File : %s\n",
VInputFiles[
f].c_str());
577 FILE* doesFileExist = fopen(
VInputFiles[
f].c_str(),
"r" );
579 printf(
"File %s doesn't exist\n",
VInputFiles[
f].c_str());
582 fclose(doesFileExist);
585 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;}
586 APV_Charge ->Add( (TH1*) file->FindObjectAny(
"APV_Charge") , 1);
587 APV_Momentum ->Add( (TH1*) file->FindObjectAny(
"APV_Momentum") , 1);
588 APV_PathLength ->Add( (TH1*) file->FindObjectAny(
"APV_PathLength") , 1);
590 Tracks_P_Vs_Eta ->Add( (TH1*) file->FindObjectAny(
"Tracks_P_Vs_Eta") , 1);
591 Tracks_Pt_Vs_Eta ->Add( (TH1*) file->FindObjectAny(
"Tracks_Pt_Vs_Eta") , 1);
601 HTrackHits ->Add( (TH1*) file->FindObjectAny(
"TrackHits") , 1);
610 Charge_Vs_Alpha ->Add( (TH1*) file->FindObjectAny(
"Charge_Vs_Alpha") , 1);
611 NStrips_Vs_Alpha ->Add( (TH1*) file->FindObjectAny(
"NStrips_Vs_Alpha") , 1);
612 HFirstStrip ->Add( (TH1*) file->FindObjectAny(
"FirstStrip") , 1);
614 TH1F* JobInfo_tmp = (TH1F*) file->FindObjectAny(
"JobInfo");
615 NEvent += (
unsigned int) JobInfo_tmp->GetBinContent(JobInfo_tmp->GetXaxis()->FindBin(1));
616 unsigned int tmp_SRun = (
unsigned int) JobInfo_tmp->GetBinContent(JobInfo_tmp->GetXaxis()->FindBin(3));
617 unsigned int tmp_SEvent = (
unsigned int) JobInfo_tmp->GetBinContent(JobInfo_tmp->GetXaxis()->FindBin(4));
618 unsigned int tmp_ERun = (
unsigned int) JobInfo_tmp->GetBinContent(JobInfo_tmp->GetXaxis()->FindBin(6));
619 unsigned int tmp_EEvent = (
unsigned int) JobInfo_tmp->GetBinContent(JobInfo_tmp->GetXaxis()->FindBin(7));
629 printf(
"Deleting Current Input File\n");
642 if( strcmp(
AlgoMode.c_str(),
"MultiJob")!=0 ){
644 double* FitResults =
new double[5];
647 if( I%3650==0 ) printf(
"Fitting Histograms \t %6.2f%%\n",(100.0*I)/
APVsColl.size());
652 Proj =
APV_Charge->ProjectionY(
" ",bin,bin,
"e");
653 Proj = (TH1D*)Proj->Clone();
654 if(Proj==
NULL)
continue;
658 int SecondAPVId = APV->
APVId;
659 if(SecondAPVId%2==0){
660 SecondAPVId = SecondAPVId+1;
662 SecondAPVId = SecondAPVId-1;
667 TH1D* Proj2 =
APV_Charge->ProjectionY(
" ",bin2,bin2,
"e");
680 TH1D* Proj2 =
APV_Charge->ProjectionY(
" ",bin2,bin2,
"e");
693 APV->
MPV = FitResults[0];
695 if(FitResults[0]!=-0.5 && FitResults[1]<
MaxMPVError){
746 if(FitResults[0]!=-0.5){
758 if(Proj==
NULL)
continue;
765 unsigned int GOOD = 0;
766 unsigned int BAD = 0;
767 double MPVmean =
MPVs->GetMean();
773 APV->
Gain = APV->
MPV / MPVmean;
943 fprintf(Gains,
"NEvents = %i\n",
NEvent);
944 fprintf(Gains,
"Number of APVs = %lu\n",static_cast<unsigned long>(
APVsColl.size()));
945 fprintf(Gains,
"GoodFits = %i BadFits = %i ratio = %f\n",GOOD,BAD,(100.0*GOOD)/(GOOD+BAD));
951 std::vector<int> DetIdOfBuggedAPV;
952 fprintf(Gains,
"----------------------------------------------------------------------\n");
955 if(APV->
MPV>0 && APV->
MPV<200){
957 for(
unsigned int b=0;
b<DetIdOfBuggedAPV.size()&&!tmpBug;
b++){
if(DetIdOfBuggedAPV[
b]==APV->
DetId)tmpBug=
true;}
958 if(!tmpBug){fprintf(Gains,
"%i,\n",APV->
DetId);DetIdOfBuggedAPV.push_back(APV->
DetId);}
980 if(!gainHandle.
isValid()){printf(
"\n#####################\n\nERROR --> gainHandle is not valid\n\n#####################\n\n");
exit(0);}
993 printf(
"GAIN HANDLE IS NOT VALID\n");
1007 if( strcmp(
AlgoMode.c_str(),
"WriteOnDB")==0 )
return;
1027 const Track track = *it->val;
1040 for(Trajectory::RecHitContainer::const_iterator rechit = transRecHits.begin(); rechit != transRecHits.end(); ++rechit)
1041 if ((*rechit)->isValid()) ndof += (*rechit)->dimension();
1048 vector<TrajectoryMeasurement> measurements = traj.
measurements();
1071 for(vector<TrajectoryMeasurement>::const_iterator measurement_it = measurements.begin(); measurement_it!=measurements.end(); measurement_it++){
1074 if( !trajState.
isValid() )
continue;
1081 if(sistripsimplehit){
1083 }
else if(sistripmatchedhit){
1086 }
else if(sistripsimple1dhit){
1100 double cosine = trackDirection.
z()/trackDirection.
mag();
1104 int APVId = FirstStrip/128;
1107 bool Overlaping =
false;
1109 unsigned int NHighStrip = 0;
1114 if(FirstStrip==0 )Overlaping=
true;
1115 if(FirstStrip==128 )Overlaping=
true;
1116 if(FirstStrip==256 )Overlaping=
true;
1117 if(FirstStrip==384 )Overlaping=
true;
1118 if(FirstStrip==512 )Overlaping=
true;
1119 if(FirstStrip==640 )Overlaping=
true;
1121 if(FirstStrip<=127 && FirstStrip+Ampls.size()>127)Overlaping=
true;
1122 if(FirstStrip<=255 && FirstStrip+Ampls.size()>255)Overlaping=
true;
1123 if(FirstStrip<=383 && FirstStrip+Ampls.size()>383)Overlaping=
true;
1124 if(FirstStrip<=511 && FirstStrip+Ampls.size()>511)Overlaping=
true;
1125 if(FirstStrip<=639 && FirstStrip+Ampls.size()>639)Overlaping=
true;
1127 if(FirstStrip+Ampls.size()==127 )Overlaping=
true;
1128 if(FirstStrip+Ampls.size()==255 )Overlaping=
true;
1129 if(FirstStrip+Ampls.size()==383 )Overlaping=
true;
1130 if(FirstStrip+Ampls.size()==511 )Overlaping=
true;
1131 if(FirstStrip+Ampls.size()==639 )Overlaping=
true;
1132 if(FirstStrip+Ampls.size()==767 )Overlaping=
true;
1133 if(Overlaping)
return -1;
1145 for(
unsigned int a=0;
a<Ampls.size();
a++){Charge+=Ampls[
a];
if(Ampls[
a]>=254)Saturation++;
if(Ampls[
a]>=20)NHighStrip++;}
1147 double ClusterChargeOverPath = (double)Charge / path ;
1163 double trans = atan2(trackDirection.
y(),trackDirection.
x())*(180/3.14159265);
1164 double alpha = acos(trackDirection.
x() /
sqrt(
pow(trackDirection.
x(),2) +
pow(trackDirection.
z(),2) ) ) * (180/3.14159265);
1165 double beta = acos(trackDirection.
y() /
sqrt(
pow(trackDirection.
x(),2) +
pow(trackDirection.
z(),2) ) ) * (180/3.14159265);
1167 if(path>0.4 && path<0.45){
1193 return ClusterChargeOverPath;
1204 if (dynamic_cast<const StripGeomDetUnit*>(it)==0 && dynamic_cast<const PixelGeomDetUnit*>(it)==0) {
1205 std::cout <<
"this detID doesn't seem to belong to the Tracker" << std::endl;
1210 const TrapezoidalPlaneBounds* trapezoidalBounds( dynamic_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
1211 const RectangularPlaneBounds* rectangularBounds( dynamic_cast<const RectangularPlaneBounds*>(&(plane.bounds())));
1213 double DistFromBorder = 1.0;
1217 if(trapezoidalBounds)
1219 std::array<const float, 4>
const &
parameters = (*trapezoidalBounds).parameters();
1220 HalfLength = parameters[3];
1223 }
else if(rectangularBounds){
1226 }
else{
return false;}
1229 if (fabs(HitLocalPos.
y())+HitLocalError.
yy() >= (HalfLength - DistFromBorder) )
return false;
1237 double adcs_err = 0.;
1238 double width = -0.5;
1239 double width_err = 0;
1240 double chi2overndf = -0.5;
1242 double nr_of_entries = InputHisto->GetEntries();
1245 FitResults[0] = adcs;
1246 FitResults[1] = adcs_err;
1247 FitResults[2] =
width;
1248 FitResults[3] = width_err;
1249 FitResults[4] = chi2overndf;
1254 TF1* MyLandau =
new TF1(
"MyLandau",
"landau",LowRange, HighRange);
1255 MyLandau->SetParameter(
"MPV",300);
1257 InputHisto->Fit(MyLandau,
"QR WW");
1260 adcs = MyLandau->GetParameter(
"MPV");
1261 adcs_err = MyLandau->GetParError(1);
1262 width = MyLandau->GetParameter(2);
1263 width_err = MyLandau->GetParError(2);
1264 chi2overndf = MyLandau->GetChisquare() / MyLandau->GetNDF();
1268 adcs = -0.5; adcs_err = 0.;
1269 width = -0.5; width_err = 0;
1273 FitResults[0] = adcs;
1274 FitResults[1] = adcs_err;
1275 FitResults[2] =
width;
1276 FitResults[3] = width_err;
1277 FitResults[4] = chi2overndf;
1287 cout <<
"START getNewObject\n";
1294 std::vector<float>* theSiStripVector =
NULL;
1295 int PreviousDetId = -1;
1299 if(APV==
NULL){ printf(
"Bug\n");
continue; }
1300 if(APV->
DetId != PreviousDetId){
1301 if(theSiStripVector!=
NULL){
1303 if ( !obj->
put(PreviousDetId,range) ) printf(
"Bug to put detId = %i\n",PreviousDetId);
1306 theSiStripVector =
new std::vector<float>;
1307 PreviousDetId = APV->
DetId;
1310 theSiStripVector->push_back(APV->
Gain);
1314 if(theSiStripVector!=
NULL){
1316 if ( !obj->
put(PreviousDetId,range) ) printf(
"Bug to put detId = %i\n",PreviousDetId);
1319 cout <<
"END getNewObject\n";
ClusterRef cluster() const
double p() const
momentum vector magnitude
T getParameter(std::string const &) const
EventNumber_t event() const
virtual int nstrips() const =0
T getUntrackedParameter(std::string const &, T const &) const
virtual void algoEndJob() override
virtual float length() const =0
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)
virtual 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
DataContainer const & measurements() const
vector< string > VInputFiles
double eta() const
pseudorapidity of momentum vector
TH1F * FWHM_Vs_PathLength
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
virtual 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
__gnu_cxx::hash_map< unsigned int, stAPVGain *, __gnu_cxx::hash< unsigned int >, isEqual > APVsColl
unsigned int MinTrackHits
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)
virtual void algoAnalyze(const edm::Event &, const edm::EventSetup &) override
TH2F * Charge_Vs_TransversAngle
T const * product() const
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
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