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());I++;
651 Proj =
APV_Charge->ProjectionY(
" ",bin,bin,
"e");
652 Proj = (TH1D*)Proj->Clone();
653 if(Proj==
NULL)
continue;
657 int SecondAPVId = APV->
APVId;
658 if(SecondAPVId%2==0){
659 SecondAPVId = SecondAPVId+1;
661 SecondAPVId = SecondAPVId-1;
666 TH1D* Proj2 =
APV_Charge->ProjectionY(
" ",bin2,bin2,
"e");
679 TH1D* Proj2 =
APV_Charge->ProjectionY(
" ",bin2,bin2,
"e");
692 APV->
MPV = FitResults[0];
694 if(FitResults[0]!=-0.5 && FitResults[1]<
MaxMPVError){
745 if(FitResults[0]!=-0.5){
757 if(Proj==
NULL)
continue;
764 unsigned int GOOD = 0;
765 unsigned int BAD = 0;
766 double MPVmean =
MPVs->GetMean();
772 APV->
Gain = APV->
MPV / MPVmean;
942 fprintf(Gains,
"NEvents = %i\n",
NEvent);
943 fprintf(Gains,
"Number of APVs = %lu\n",static_cast<unsigned long>(
APVsColl.size()));
944 fprintf(Gains,
"GoodFits = %i BadFits = %i ratio = %f\n",GOOD,BAD,(100.0*GOOD)/(GOOD+BAD));
950 std::vector<int> DetIdOfBuggedAPV;
951 fprintf(Gains,
"----------------------------------------------------------------------\n");
954 if(APV->
MPV>0 && APV->
MPV<200){
956 for(
unsigned int b=0;
b<DetIdOfBuggedAPV.size()&&!tmpBug;
b++){
if(DetIdOfBuggedAPV[
b]==APV->
DetId)tmpBug=
true;}
957 if(!tmpBug){fprintf(Gains,
"%i,\n",APV->
DetId);DetIdOfBuggedAPV.push_back(APV->
DetId);}
979 if(!gainHandle.
isValid()){printf(
"\n#####################\n\nERROR --> gainHandle is not valid\n\n#####################\n\n");
exit(0);}
992 printf(
"GAIN HANDLE IS NOT VALID\n");
1006 if( strcmp(
AlgoMode.c_str(),
"WriteOnDB")==0 )
return;
1026 const Track track = *it->val;
1039 for(Trajectory::RecHitContainer::const_iterator rechit = transRecHits.begin(); rechit != transRecHits.end(); ++rechit)
1040 if ((*rechit)->isValid()) ndof += (*rechit)->dimension();
1047 vector<TrajectoryMeasurement> measurements = traj.
measurements();
1070 for(vector<TrajectoryMeasurement>::const_iterator measurement_it = measurements.begin(); measurement_it!=measurements.end(); measurement_it++){
1073 if( !trajState.
isValid() )
continue;
1076 const SiStripRecHit1D* sistripsimple1dhit =
dynamic_cast<const SiStripRecHit1D*
>(hit);
1077 const SiStripRecHit2D* sistripsimplehit =
dynamic_cast<const SiStripRecHit2D*
>(hit);
1078 const SiStripMatchedRecHit2D* sistripmatchedhit =
dynamic_cast<const SiStripMatchedRecHit2D*
>(hit);
1080 if(sistripsimplehit){
1082 }
else if(sistripmatchedhit){
1085 }
else if(sistripsimple1dhit){
1099 double cosine = trackDirection.
z()/trackDirection.
mag();
1102 const vector<uint8_t>& Ampls = Cluster->
amplitudes();
1105 int APVId = FirstStrip/128;
1108 bool Overlaping =
false;
1110 unsigned int NHighStrip = 0;
1115 if(FirstStrip==0 )Overlaping=
true;
1116 if(FirstStrip==128 )Overlaping=
true;
1117 if(FirstStrip==256 )Overlaping=
true;
1118 if(FirstStrip==384 )Overlaping=
true;
1119 if(FirstStrip==512 )Overlaping=
true;
1120 if(FirstStrip==640 )Overlaping=
true;
1122 if(FirstStrip<=127 && FirstStrip+Ampls.size()>127)Overlaping=
true;
1123 if(FirstStrip<=255 && FirstStrip+Ampls.size()>255)Overlaping=
true;
1124 if(FirstStrip<=383 && FirstStrip+Ampls.size()>383)Overlaping=
true;
1125 if(FirstStrip<=511 && FirstStrip+Ampls.size()>511)Overlaping=
true;
1126 if(FirstStrip<=639 && FirstStrip+Ampls.size()>639)Overlaping=
true;
1128 if(FirstStrip+Ampls.size()==127 )Overlaping=
true;
1129 if(FirstStrip+Ampls.size()==255 )Overlaping=
true;
1130 if(FirstStrip+Ampls.size()==383 )Overlaping=
true;
1131 if(FirstStrip+Ampls.size()==511 )Overlaping=
true;
1132 if(FirstStrip+Ampls.size()==639 )Overlaping=
true;
1133 if(FirstStrip+Ampls.size()==767 )Overlaping=
true;
1134 if(Overlaping)
return -1;
1146 for(
unsigned int a=0;
a<Ampls.size();
a++){Charge+=Ampls[
a];
if(Ampls[
a]>=254)Saturation++;
if(Ampls[
a]>=20)NHighStrip++;}
1148 double ClusterChargeOverPath = (double)Charge / path ;
1164 double trans = atan2(trackDirection.
y(),trackDirection.
x())*(180/3.14159265);
1165 double alpha = acos(trackDirection.
x() /
sqrt(
pow(trackDirection.
x(),2) +
pow(trackDirection.
z(),2) ) ) * (180/3.14159265);
1166 double beta = acos(trackDirection.
y() /
sqrt(
pow(trackDirection.
x(),2) +
pow(trackDirection.
z(),2) ) ) * (180/3.14159265);
1168 if(path>0.4 && path<0.45){
1194 return ClusterChargeOverPath;
1205 if (dynamic_cast<const StripGeomDetUnit*>(it)==0 && dynamic_cast<const PixelGeomDetUnit*>(it)==0) {
1206 std::cout <<
"this detID doesn't seem to belong to the Tracker" << std::endl;
1211 const TrapezoidalPlaneBounds* trapezoidalBounds( dynamic_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
1212 const RectangularPlaneBounds* rectangularBounds( dynamic_cast<const RectangularPlaneBounds*>(&(plane.bounds())));
1214 double DistFromBorder = 1.0;
1218 if(trapezoidalBounds)
1220 std::array<const float, 4>
const &
parameters = (*trapezoidalBounds).parameters();
1221 HalfLength = parameters[3];
1224 }
else if(rectangularBounds){
1227 }
else{
return false;}
1230 if (fabs(HitLocalPos.
y())+HitLocalError.
yy() >= (HalfLength - DistFromBorder) )
return false;
1238 double adcs_err = 0.;
1239 double width = -0.5;
1240 double width_err = 0;
1241 double chi2overndf = -0.5;
1243 double nr_of_entries = InputHisto->GetEntries();
1246 FitResults[0] = adcs;
1247 FitResults[1] = adcs_err;
1248 FitResults[2] =
width;
1249 FitResults[3] = width_err;
1250 FitResults[4] = chi2overndf;
1255 TF1* MyLandau =
new TF1(
"MyLandau",
"landau",LowRange, HighRange);
1256 MyLandau->SetParameter(
"MPV",300);
1258 InputHisto->Fit(
"MyLandau",
"QR WW");
1259 TF1 * fitfunction = (TF1*) InputHisto->GetListOfFunctions()->First();
1262 adcs = fitfunction->GetParameter(
"MPV");
1263 adcs_err = fitfunction->GetParError(1);
1264 width = fitfunction->GetParameter(2);
1265 width_err = fitfunction->GetParError(2);
1266 chi2overndf = fitfunction->GetChisquare() / fitfunction->GetNDF();
1270 adcs = -0.5; adcs_err = 0.;
1271 width = -0.5; width_err = 0;
1275 FitResults[0] = adcs;
1276 FitResults[1] = adcs_err;
1277 FitResults[2] =
width;
1278 FitResults[3] = width_err;
1279 FitResults[4] = chi2overndf;
1289 cout <<
"START getNewObject\n";
1296 std::vector<float>* theSiStripVector =
NULL;
1297 int PreviousDetId = -1;
1301 if(APV==
NULL){ printf(
"Bug\n");
continue; }
1302 if(APV->
DetId != PreviousDetId){
1303 if(theSiStripVector!=
NULL){
1305 if ( !obj->
put(PreviousDetId,range) ) printf(
"Bug to put detId = %i\n",PreviousDetId);
1308 theSiStripVector =
new std::vector<float>;
1309 PreviousDetId = APV->
DetId;
1312 theSiStripVector->push_back(APV->
Gain);
1316 if(theSiStripVector!=
NULL){
1318 if ( !obj->
put(PreviousDetId,range) ) printf(
"Bug to put detId = %i\n",PreviousDetId);
1321 cout <<
"END getNewObject\n";
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
MonitorElement * book1D(const char *name, const char *title, int nchX, double lowX, double highX)
Book 1D histogram.
friend struct const_iterator
const_iterator end() const
last iterator over the map (read only)
unsigned int MinNrEntries
std::vector< stAPVGain * > APVsCollOrdered
void cd(void)
go to top directory (ie. root)
#define DEFINE_FWK_MODULE(type)
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
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
1: Failed selection (without additional info)
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
tuple path
else: Piece not in the list, fine.
double pt() const
track transverse momentum
std::pair< ContainerIterator, ContainerIterator > Range
const std::complex< double > I
TH1F * MPV_Vs_PathLength320
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
TH2F * Charge_Vs_PathLength
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
ConstRecHitContainer RecHitContainer
TH2F * Charge_Vs_TransversAngle
T const * product() const
TH1F * getTH1F(void) const
T const * product() const
void getPeakOfLandau(TH1 *InputHisto, double *FitResults, double LowRange=0, double HighRange=5400)
std::string TrajToTrackLabel
TH2F * Charge_Vs_PathLength320
GlobalVector globalMomentum() const
void save(const std::string &filename, const std::string &path="", const std::string &pattern="", const std::string &rewrite="", const uint32_t run=0, const uint32_t lumi=0, SaveReferenceTag ref=SaveWithReference, int minStatus=dqm::qstatus::STATUS_OK, const std::string &fileupdate="RECREATE", const bool resetMEsAfterWriting=false)
TH1F * NHighStripInCluster
const_iterator begin() const
first iterator over the map (read only)
TH2F * getTH2F(void) const
MonitorElement * book2D(const char *name, const char *title, int nchX, double lowX, double highX, int nchY, double lowY, double highY)
Book 2D histogram.
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