423 int irun = iEvent.
id().
run();
444 int tmpeta1 = (
i>0) ?
i -1 : -
i +14;
445 if (tmpeta1 <0 || tmpeta1 >netamx)
continue;
488 if ((*ho).size()>0) {
491 m_coder = (*conditions_).getHcalCoder(
id);
493 int tmpeta=
id.ieta();
494 int tmpphi=
id.iphi();
496 int tmpeta1 = (tmpeta>0) ? tmpeta -1 : -tmpeta +14;
497 for (
int i=0;
i<(*j).size() &&
i<
nchnmx;
i++) {
499 allhotime->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) +
i, tmpdata[
i]);
500 Nallhotime->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + i, 1.);
504 if ((*hbhe).size()>0) {
507 m_coder = (*conditions_).getHcalCoder(
id);
509 int tmpeta=
id.ieta();
510 int tmpphi=
id.iphi();
511 int tmpdepth =
id.depth();
512 int tmpeta1 = (tmpeta>0) ? tmpeta -15 : -tmpeta + 1;
513 if (tmpdepth==1) tmpeta1 = (tmpeta>0) ? tmpeta -1 : -tmpeta +29;
514 for (
int i=0;
i<(*j).size() &&
i<
nchnmx;
i++) {
517 allhb1->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) +
i, signal);
518 Nallhb1->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) +
i, 1);
519 hb1pedpr->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) +
i, signal);}
521 allhb2->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) +
i, signal);
522 Nallhb2->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) +
i, 1);}
524 allhb3->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) +
i, signal);
525 Nallhb3->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) +
i, 1);}
531 double pival = acos(-1.);
536 if (cosmicmuon->size()>0) {
541 int ntrgpas_gm[
ntrgp_gm]={0,0,0,0,0,0,0,0,0,0};
582 int Noccu_old =
Noccu;
584 for(reco::TrackCollection::const_iterator ncosm = cosmicmuon->begin();
585 ncosm != cosmicmuon->end(); ++ncosm) {
587 if ((*ncosm).ndof() < 15)
continue;
588 if ((*ncosm).normalizedChi2() >30.0)
continue;
592 tmpHOCalib.
trig1 = l1trg;
593 tmpHOCalib.
trig2 = hlttr;
595 int charge = ncosm->charge();
597 double innerr = (*ncosm).innerPosition().Perp2();
598 double outerr = (*ncosm).outerPosition().Perp2();
599 int iiner = (innerr <outerr) ? 1 : 0;
610 double posx, posy, posz;
611 double momx, momy, momz;
614 posx = (*ncosm).innerPosition().X();
615 posy = (*ncosm).innerPosition().Y();
616 posz = (*ncosm).innerPosition().Z();
618 momx = (*ncosm).innerMomentum().X();
619 momy = (*ncosm).innerMomentum().Y();
620 momz = (*ncosm).innerMomentum().Z();
623 posx = (*ncosm).outerPosition().X();
624 posy = (*ncosm).outerPosition().Y();
625 posz = (*ncosm).outerPosition().Z();
627 momx = (*ncosm).outerMomentum().X();
628 momy = (*ncosm).outerMomentum().Y();
629 momz = (*ncosm).outerMomentum().Z();
635 CLHEP::Hep3Vector tmpmuon3v(posx, posy, posz);
636 CLHEP::Hep3Vector tmpmuondir(momx, momy, momz);
638 bool samedir = (tmpmuon3v.dot(tmpmuondir) >0) ?
true :
false;
639 for (
int i=0;
i<3;
i++) {tmpHOCalib.
caloen[
i] = 0.0;}
641 for(reco::TrackCollection::const_iterator ncosmcor = cosmicmuon->begin();
642 ncosmcor != cosmicmuon->end(); ++ncosmcor) {
643 if (ncosmcor==ncosm)
continue;
645 CLHEP::Hep3Vector tmpmuon3vcor;
646 CLHEP::Hep3Vector tmpmom3v;
648 tmpmuon3vcor = CLHEP::Hep3Vector((*ncosmcor).innerPosition().X(),(*ncosmcor).innerPosition().Y(),(*ncosmcor).innerPosition().Z());
649 tmpmom3v = CLHEP::Hep3Vector((*ncosmcor).innerMomentum().X(),(*ncosmcor).innerMomentum().Y(),(*ncosmcor).innerMomentum().Z());
651 tmpmuon3vcor = CLHEP::Hep3Vector((*ncosmcor).outerPosition().X(),(*ncosmcor).outerPosition().Y(),(*ncosmcor).outerPosition().Z());
652 tmpmom3v = CLHEP::Hep3Vector((*ncosmcor).outerMomentum().X(),(*ncosmcor).outerMomentum().Y(),(*ncosmcor).outerMomentum().Z());
655 if (tmpmom3v.mag()<0.2 || (*ncosmcor).ndof()<5)
continue;
657 double angle = tmpmuon3v.angle(tmpmuon3vcor);
658 if (angle < 7.5*pival/180.) {inearbymuon=1;}
661 if (angle <7.5*pival/180.) { tmpHOCalib.
caloen[0] +=1.;}
662 if (angle <15.0*pival/180.) { tmpHOCalib.
caloen[1] +=1.;}
663 if (angle <35.0*pival/180.) { tmpHOCalib.
caloen[2] +=1.;}
675 calt !=calotower->end(); calt++) {
677 double ith = (*calt).momentum().theta();
678 double iph = (*calt).momentum().phi();
680 CLHEP::Hep3Vector calo3v(
sin(ith)*
cos(iph),
sin(ith)*
sin(iph),
cos(ith));
682 double angle = tmpmuon3v.angle(calo3v);
684 if (angle < 7.5*pival/180.) {tmpHOCalib.
caloen[0] += calt->emEnergy()+calt->hadEnergy();}
685 if (angle < 15*pival/180.) {tmpHOCalib.
caloen[1] += calt->emEnergy()+calt->hadEnergy();}
686 if (angle < 35*pival/180.) {tmpHOCalib.
caloen[2] += calt->emEnergy()+calt->hadEnergy();}
691 if (tmpHOCalib.
caloen[0] >10.0)
continue;
695 double mom =
sqrt(momx*momx + momy*momy +momz*momz);
703 tmpHOCalib.
trkdr = (*ncosm).d0();
704 tmpHOCalib.
trkdz = (*ncosm).dz();
706 tmpHOCalib.
nmuon = cosmicmuon->size();
707 tmpHOCalib.
trkvx = glbpt.x();
708 tmpHOCalib.
trkvy = glbpt.y();
709 tmpHOCalib.
trkvz = glbpt.z();
711 tmpHOCalib.
trkth = trkdir.theta();
712 tmpHOCalib.
trkph = trkdir.phi();
714 tmpHOCalib.
ndof = (inearbymuon ==0) ? (
int)(*ncosm).ndof() : -(int)(*ncosm).ndof();
715 tmpHOCalib.
chisq = (*ncosm).normalizedChi2();
716 tmpHOCalib.
therr = 0.;
717 tmpHOCalib.
pherr = 0.;
721 tmpHOCalib.
therr = innercov(1,1);
722 tmpHOCalib.
pherr = innercov(2,2);
725 tmpHOCalib.
therr = outercov(1,1);
726 tmpHOCalib.
pherr = outercov(2,2);
733 myHelix.setMaterialMode(
false);
734 myHelix.applyRadX0Correction(
true);
736 double phiho = trkpos.phi();
737 if (phiho<0) phiho +=2*pival;
739 int iphisect_dt=int(6*(phiho+pival/18.)/pival);
740 if (iphisect_dt>=12) iphisect_dt=0;
745 for (
int kl = 0; kl<=2; kl++) {
747 int iphisecttmp = (kl<2) ? iphisect_dt + kl : iphisect_dt - 1;
748 if (iphisecttmp <0) iphisecttmp = 11;
749 if (iphisecttmp >=12) iphisecttmp = 0;
751 double phipos = iphisecttmp*pival/6.;
752 double phirot = phipos;
765 for (
int ik=1; ik>=0; ik--) {
767 double radial = 407.0;
768 if (ik==0) radial = 382.0;
778 if (steppingHelixstateinfo_.
isValid()) {
781 CLHEP::Hep3Vector hotrkdir2(steppingHelixstateinfo_.
momentum().
x(), steppingHelixstateinfo_.
momentum().
y(),steppingHelixstateinfo_.
momentum().
z());
785 double xx = lclvt0.
x();
786 double yy = lclvt0.y();
789 if ((
std::abs(yy) < 130 && xx >-64.7 && xx <138.2)
792 iphisect = iphisecttmp;
796 if (iphisect != iphisecttmp)
continue;
808 tmpHOCalib.
hoang = CLHEP::Hep3Vector(zLocal.
x(),zLocal.
y(),zLocal.
z()).
dot(hotrkdir2.unit());
823 for (
int i=0;
i<9;
i++) {tmpHOCalib.
hosig[
i]=-100.0;}
824 for (
int i=0;
i<18;
i++) {tmpHOCalib.
hocorsig[
i]=-100.0;}
825 for (
int i=0;
i<9;
i++) {tmpHOCalib.
hbhesig[
i]=-100.0;}
826 tmpHOCalib.
hocro = -100;
827 tmpHOCalib.
htime = -1000;
837 if (iphiho<0) isect -=2000000;
838 tmpHOCalib.
isect = isect;
850 if (
iring==1) {etamn=5; etamx = 10;}
851 if (
iring==2) {etamn=11; etamx = 16;}
852 if (
iring==-1){etamn=-10; etamx = -5;}
853 if (
iring==-2){etamn=-16; etamx = -11;}
858 phimx =2*int((iphiho+1)/2.);
861 phimn = 3*int((iphiho+1)/3.) - 1;
865 if (phimn <1) phimn +=
nphimx;
866 if (phimx >72) phimx -=
nphimx;
877 for (
int i=0;
i<9;
i++) {tmpHOCalib.
hbhesig[
i]=-100.0;}
880 if ((*hbhe).size() >0) {
885 m_coder = (*conditions_).getHcalCoder(
id);
887 int tmpeta=
id.ieta();
888 int tmpphi=
id.iphi();
891 int deta = tmpeta-ietaho;
892 if (tmpeta==-1 && ietaho== 1) deta = -1;
893 if (tmpeta== 1 && ietaho==-1) deta = 1;
894 int dphi = tmpphi-iphiho;
896 if (dphi==71) dphi=-1;
897 if (dphi==-71) dphi=1;
902 if (ipass2 ==0 )
continue;
905 for (
int i=0;
i<(*j).size() &&
i<
nchnmx;
i++) {
910 for (
int i=1;
i<(*j).size() &&
i<=8;
i++) {
915 if (3*(deta+1)+dphi+1<9) tmpHOCalib.
hbhesig[3*(deta+1)+dphi+1] = signal;
926 if ((*hbheht).size()>0) {
927 if(!(*hbheht).size())
throw (
int)(*hbheht).size();
933 int tmpeta=
id.
ieta();
934 int tmpphi=
id.iphi();
936 int deta = tmpeta-ietaho;
937 if (tmpeta==-1 && ietaho== 1) deta = -1;
938 if (tmpeta== 1 && ietaho==-1) deta = 1;
939 int dphi = tmpphi-iphiho;
941 if (dphi==71) dphi=-1;
942 if (dphi==-71) dphi=1;
946 if ( ipass2 ==0 )
continue;
948 float signal = (*j).energy();
950 if (3*(deta+1)+dphi+1<9) tmpHOCalib.
hbhesig[3*(deta+1)+dphi+1] = signal;
959 if ((*ho).size()>0) {
960 int isFilled[netamx*
nphimx];
961 for (
int j=0;
j<netamx*
nphimx;
j++) {isFilled[
j]=0;}
971 m_coder = (*conditions_).getHcalCoder(
id);
973 int tmpeta=
id.ieta();
974 int tmpphi=
id.iphi();
977 if (tmpeta >=etamn && tmpeta <=etamx) {
979 ipass1 = (tmpphi >=phimn && tmpphi <=phimx ) ? 1 : 0;
981 ipass1 = (tmpphi==71 || tmpphi ==72 || tmpphi==1) ? 1 : 0;
985 int deta = tmpeta-ietaho;
986 if (tmpeta==-1 && ietaho== 1) deta = -1;
987 if (tmpeta== 1 && ietaho==-1) deta = 1;
989 int dphi = tmpphi -iphiho;
991 if (dphi==71) dphi=-1;
992 if (dphi==-71) dphi=1;
997 int tmpeta1 = (tmpeta>0) ? tmpeta -1 : -tmpeta +14;
999 float tmpdata[
nchnmx]={0,0,0,0,0,0,0,0,0,0};
1000 float sigvall[
nsigpk]={0,0,0,0,0,0,0};
1002 for (
int i=0;
i<(*j).size() &&
i<
nchnmx;
i++) {
1004 if (deta==0 && dphi==0) {
1005 double tmpE = tmpdata[
i] -
pedestal[tmpeta1][tmpphi-1][(*j).sample(
i).capid()];
1013 for (
int ncap=0; ncap<
nsigpk; ncap++) {
1015 sigvall[ncap] +=tmpdata[
i];
1019 if (
i==(*j).size()-1) {
1022 for (
int ij=0; ij<
nsigpk; ij++) {
1023 if (sigvall[ij] > mxled) {mxled = sigvall[ij]; imxled=ij;}
1026 for (
int ij=0; ij<4; ij++) {
1027 pedx +=pedestal[tmpeta1][tmpphi-1][ij];
1029 if (mxled-pedx >2 && mxled-pedx <20 ) {
1031 for (
int jk=0; jk<
ntrgp_gm; jk++) {
1032 if (ntrgpas_gm[jk]>0) {
1033 hopeak[jk]->Fill(nphimx*tmpeta1 + tmpphi-1, imxled+
nstrbn);
1036 if (tmpdata[5]+tmpdata[6] >1) {
1037 horatio->Fill(nphimx*tmpeta1 + tmpphi-1, (tmpdata[5]-tmpdata[6])/(tmpdata[5]+tmpdata[6]));
1039 for (
int ij=0; ij<(*j).size() && ij<
nchnmx; ij++) {
1040 hotime[
ntrgp_gm]->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + ij, tmpdata[ij]);
1041 Nhotime[
ntrgp_gm]->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + ij, 1.);
1042 for (
int jk=0; jk<
ntrgp_gm; jk++) {
1043 if (ntrgpas_gm[jk]>0) {
1044 hotime[jk]->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + ij, tmpdata[ij]);
1045 Nhotime[jk]->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + ij, 1.);
1055 if (
std::abs(tmpeta) <=15 && deta==0 && dphi ==0) {
1058 for (
int i =0;
i<nchnmx &&
i< (*j).size();
i++) {
1059 if (
i >=sigstr &&
i<=sigend)
continue;
1060 signal += tmpdata[
i] - pedestal[tmpeta1][tmpphi-1][(*j).sample(
i).capid()];
1061 if (++icnt >=4)
break;
1063 tmpHOCalib.
hocro = signal;
1067 if (ipass1 ==0 && ipass2 ==0 && cosmicmuon->size()<=2) {
1069 if ((iphiho >=1 && iphiho<=nphimx/2 && tmpphi >=1 && tmpphi <=nphimx/2) ||
1070 (iphiho >nphimx/2 && iphiho<=nphimx && tmpphi >nphimx/2 && tmpphi <=
nphimx)) {
1071 if (isFilled[nphimx*tmpeta1+tmpphi-1]==0) {
1072 isFilled[nphimx*tmpeta1+tmpphi-1]=1;
1073 for (
int i=0;
i<(*j).size() &&
i<
nchnmx;
i++) {
1074 hopedtime->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) +
i, tmpdata[
i]);
1075 Nhopedtime->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + i, 1.);
1076 hopedpr->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + i, tmpdata[i]);
1084 if (ipass1 ==0 && ipass2 ==0 )
continue;
1087 for (
int i=sigstr;
i<(*j).size() &&
i<=sigend;
i++) {
1088 signal += tmpdata[
i] - pedestal[tmpeta1][tmpphi-1][(*j).sample(
i).capid()];
1090 if (signal <-100 || signal >100000) signal = -100;
1092 if (signal >-100 &&
Noccu == Noccu_old) {
1093 for (
int i=0;
i<5;
i++) {
1101 if (ipass1 ==0 && ipass2 ==0 )
continue;
1104 int tmpdph = tmpphi-phimn;
1105 if (tmpdph<0) tmpdph = 2;
1107 int ilog = 2*(tmpeta-etamn)+tmpdph;
1110 ilog = 3*(tmpeta-etamn)+tmpdph;
1112 ilog = 3*(etamx-tmpeta)+tmpdph;
1115 if (ilog>-1 && ilog<18) {
1116 tmpHOCalib.
hocorsig[ilog] = signal;
1121 if (3*(deta+1)+dphi+1<9) tmpHOCalib.
hosig[3*(deta+1)+dphi+1] = signal;
1155 tmpHOCalib.
htime = sumEt/
max(sumE,1.
e-6);
1162 if ((*hoht).size()>0) {
1167 int tmpeta=
id.
ieta();
1168 int tmpphi=
id.iphi();
1171 if (tmpeta >=etamn && tmpeta <=etamx) {
1172 if (phimn < phimx) {
1173 ipass1 = (tmpphi >=phimn && tmpphi <=phimx ) ? 1 : 0;
1175 ipass1 = (tmpphi==71 || tmpphi ==72 || tmpphi==1) ? 1 : 0;
1179 int deta = tmpeta-ietaho;
1180 if (tmpeta==-1 && ietaho== 1) deta = -1;
1181 if (tmpeta== 1 && ietaho==-1) deta = 1;
1183 int dphi = tmpphi -iphiho;
1185 if (dphi==71) dphi=-1;
1186 if (dphi==-71) dphi=1;
1189 float signal = (*j).energy();
1191 int tmpeta1 = (tmpeta>0) ? tmpeta -1 : -tmpeta +14;
1192 if (signal >-100 &&
Noccu == Noccu_old) {
1193 for (
int i=0;
i<5;
i++) {
1203 if (ipass1 ==0 && ipass2 ==0 )
continue;
1206 int tmpdph = tmpphi-phimn;
1207 if (tmpdph<0) tmpdph = 2;
1209 int ilog = 2*(tmpeta-etamn)+tmpdph;
1212 ilog = 3*(tmpeta-etamn)+tmpdph;
1214 ilog = 3*(etamx-tmpeta)+tmpdph;
1217 if (ilog>-1 && ilog<18) {
1218 tmpHOCalib.
hocorsig[ilog] = signal;
1224 if (3*(deta+1)+dphi+1<9) {
1225 tmpHOCalib.
hosig[3*(deta+1)+dphi+1] = signal;
1229 if (deta==0 && dphi ==0) {
1230 tmpHOCalib.
htime = (*j).time();
1231 int crphi = tmpphi + 6;
1232 if (crphi >72) crphi -=72;
1237 int etacr= idcr.
ieta();
1238 int phicr= idcr.
iphi();
1239 if (tmpeta==etacr && crphi ==phicr) {
1252 hostore->push_back(tmpHOCalib);
1260 iEvent.
put(hostore,
"HOCalibVariableCollection");
Basic3DVector< float > DirectionType
TProfile * hopeak[ntrgp_gm+1]
edm::EDGetTokenT< CaloTowerCollection > tok_tower_
FreeTrajectoryState getFreeTrajectoryState(const reco::Track &tk, const MagneticField *field, int itag, bool dir)
edm::EDGetTokenT< reco::TrackCollection > tok_muons_
bool getByToken(EDGetToken token, Handle< PROD > &result) const
HcalCalibrations calibped
Sin< T >::type sin(const T &t)
ReturnType plane(const PositionType &pos, const RotationType &rot) const
std::vector< HODataFrame >::const_iterator const_iterator
edm::ESHandle< HcalDbService > conditions_
double pedestal(int fCapId) const
get pedestal for capid=0..3
const HcalQIEShape * m_shape
GlobalVector momentum() const
void findHOEtaPhi(int iphsect, int &ietaho, int &iphiho)
const HcalQIECoder * m_coder
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
GlobalPoint position() const
Vector3DBase< typename PreciseFloatType< T, U >::Type, FrameTag > cross(const Vector3DBase< U, FrameTag > &v) const
Cos< T >::type cos(const T &t)
int ieta() const
get the cell ieta
Abs< T >::type abs(const T &t)
Basic3DVector< float > PositionType
float pedestal[netamx][nphimx][ncidmx]
int iphi() const
get the cell iphi
TH1F * Nhotime[ntrgp_gm+1]
std::vector< HOCalibVariables > HOCalibVariableCollection
collection of HOcalibration variabale
edm::EDGetTokenT< HORecHitCollection > tok_ho_
T dot(const Basic3DVector &v) const
Scalar product, or "dot" product, with a vector of same type.
TH1F * hotime[ntrgp_gm+1]
HcalDetId id() const
get the id
edm::EDGetTokenT< HBHERecHitCollection > tok_hbhe_
math::Error< dimension >::type CovarianceMatrix
5 parameter covariance matrix
T angle(T x1, T y1, T z1, T x2, T y2, T z2)
float charge(const HcalQIEShape &fShape, unsigned fAdc, unsigned fCapId) const
ADC [0..127] + capid [0..3] -> fC conversion.
Public Types inherited from edm::EDProducer
1.8.5
