427 int irun = iEvent.
id().
run();
448 int tmpeta1 = (
i>0) ?
i -1 : -
i +14;
449 if (tmpeta1 <0 || tmpeta1 >netamx)
continue;
492 if ((*ho).size()>0) {
495 m_coder = (*conditions_).getHcalCoder(
id);
497 int tmpeta=
id.ieta();
498 int tmpphi=
id.iphi();
500 int tmpeta1 = (tmpeta>0) ? tmpeta -1 : -tmpeta +14;
501 for (
int i=0;
i<(*j).size() &&
i<
nchnmx;
i++) {
503 allhotime->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) +
i, tmpdata[
i]);
504 Nallhotime->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + i, 1.);
508 if ((*hbhe).size()>0) {
511 m_coder = (*conditions_).getHcalCoder(
id);
513 int tmpeta=
id.ieta();
514 int tmpphi=
id.iphi();
515 int tmpdepth =
id.depth();
516 int tmpeta1 = (tmpeta>0) ? tmpeta -15 : -tmpeta + 1;
517 if (tmpdepth==1) tmpeta1 = (tmpeta>0) ? tmpeta -1 : -tmpeta +29;
518 for (
int i=0;
i<(*j).size() &&
i<
nchnmx;
i++) {
521 allhb1->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) +
i, signal);
522 Nallhb1->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) +
i, 1);
523 hb1pedpr->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) +
i, signal);}
525 allhb2->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) +
i, signal);
526 Nallhb2->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) +
i, 1);}
528 allhb3->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) +
i, signal);
529 Nallhb3->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) +
i, 1);}
535 double pival = acos(-1.);
540 if (cosmicmuon->size()>0) {
545 int ntrgpas_gm[
ntrgp_gm]={0,0,0,0,0,0,0,0,0,0};
586 int Noccu_old =
Noccu;
588 for(reco::TrackCollection::const_iterator ncosm = cosmicmuon->begin();
589 ncosm != cosmicmuon->end(); ++ncosm) {
591 if ((*ncosm).ndof() < 15)
continue;
592 if ((*ncosm).normalizedChi2() >30.0)
continue;
596 tmpHOCalib.
trig1 = l1trg;
597 tmpHOCalib.
trig2 = hlttr;
599 int charge = ncosm->charge();
601 double innerr = (*ncosm).innerPosition().Perp2();
602 double outerr = (*ncosm).outerPosition().Perp2();
603 int iiner = (innerr <outerr) ? 1 : 0;
614 double posx, posy, posz;
615 double momx, momy, momz;
618 posx = (*ncosm).innerPosition().X();
619 posy = (*ncosm).innerPosition().Y();
620 posz = (*ncosm).innerPosition().Z();
622 momx = (*ncosm).innerMomentum().X();
623 momy = (*ncosm).innerMomentum().Y();
624 momz = (*ncosm).innerMomentum().Z();
627 posx = (*ncosm).outerPosition().X();
628 posy = (*ncosm).outerPosition().Y();
629 posz = (*ncosm).outerPosition().Z();
631 momx = (*ncosm).outerMomentum().X();
632 momy = (*ncosm).outerMomentum().Y();
633 momz = (*ncosm).outerMomentum().Z();
639 CLHEP::Hep3Vector tmpmuon3v(posx, posy, posz);
640 CLHEP::Hep3Vector tmpmuondir(momx, momy, momz);
642 bool samedir = (tmpmuon3v.dot(tmpmuondir) >0) ?
true :
false;
643 for (
int i=0;
i<3;
i++) {tmpHOCalib.
caloen[
i] = 0.0;}
645 for(reco::TrackCollection::const_iterator ncosmcor = cosmicmuon->begin();
646 ncosmcor != cosmicmuon->end(); ++ncosmcor) {
647 if (ncosmcor==ncosm)
continue;
649 CLHEP::Hep3Vector tmpmuon3vcor;
650 CLHEP::Hep3Vector tmpmom3v;
652 tmpmuon3vcor = CLHEP::Hep3Vector((*ncosmcor).innerPosition().X(),(*ncosmcor).innerPosition().Y(),(*ncosmcor).innerPosition().Z());
653 tmpmom3v = CLHEP::Hep3Vector((*ncosmcor).innerMomentum().X(),(*ncosmcor).innerMomentum().Y(),(*ncosmcor).innerMomentum().Z());
655 tmpmuon3vcor = CLHEP::Hep3Vector((*ncosmcor).outerPosition().X(),(*ncosmcor).outerPosition().Y(),(*ncosmcor).outerPosition().Z());
656 tmpmom3v = CLHEP::Hep3Vector((*ncosmcor).outerMomentum().X(),(*ncosmcor).outerMomentum().Y(),(*ncosmcor).outerMomentum().Z());
659 if (tmpmom3v.mag()<0.2 || (*ncosmcor).ndof()<5)
continue;
661 double angle = tmpmuon3v.angle(tmpmuon3vcor);
662 if (angle < 7.5*pival/180.) {inearbymuon=1;}
665 if (angle <7.5*pival/180.) { tmpHOCalib.
caloen[0] +=1.;}
666 if (angle <15.0*pival/180.) { tmpHOCalib.
caloen[1] +=1.;}
667 if (angle <35.0*pival/180.) { tmpHOCalib.
caloen[2] +=1.;}
679 calt !=calotower->end(); calt++) {
681 double ith = (*calt).momentum().theta();
682 double iph = (*calt).momentum().phi();
684 CLHEP::Hep3Vector calo3v(
sin(ith)*
cos(iph),
sin(ith)*
sin(iph),
cos(ith));
686 double angle = tmpmuon3v.angle(calo3v);
688 if (angle < 7.5*pival/180.) {tmpHOCalib.
caloen[0] += calt->emEnergy()+calt->hadEnergy();}
689 if (angle < 15*pival/180.) {tmpHOCalib.
caloen[1] += calt->emEnergy()+calt->hadEnergy();}
690 if (angle < 35*pival/180.) {tmpHOCalib.
caloen[2] += calt->emEnergy()+calt->hadEnergy();}
695 if (tmpHOCalib.
caloen[0] >10.0)
continue;
699 double mom =
sqrt(momx*momx + momy*momy +momz*momz);
707 tmpHOCalib.
trkdr = (*ncosm).d0();
708 tmpHOCalib.
trkdz = (*ncosm).dz();
710 tmpHOCalib.
nmuon = cosmicmuon->size();
711 tmpHOCalib.
trkvx = glbpt.x();
712 tmpHOCalib.
trkvy = glbpt.y();
713 tmpHOCalib.
trkvz = glbpt.z();
715 tmpHOCalib.
trkth = trkdir.theta();
716 tmpHOCalib.
trkph = trkdir.phi();
718 tmpHOCalib.
ndof = (inearbymuon ==0) ? (
int)(*ncosm).ndof() : -(int)(*ncosm).ndof();
719 tmpHOCalib.
chisq = (*ncosm).normalizedChi2();
720 tmpHOCalib.
therr = 0.;
721 tmpHOCalib.
pherr = 0.;
725 tmpHOCalib.
therr = innercov(1,1);
726 tmpHOCalib.
pherr = innercov(2,2);
729 tmpHOCalib.
therr = outercov(1,1);
730 tmpHOCalib.
pherr = outercov(2,2);
737 myHelix.setMaterialMode(
false);
738 myHelix.applyRadX0Correction(
true);
740 double phiho = trkpos.phi();
741 if (phiho<0) phiho +=2*pival;
743 int iphisect_dt=int(6*(phiho+pival/18.)/pival);
744 if (iphisect_dt>=12) iphisect_dt=0;
749 for (
int kl = 0; kl<=2; kl++) {
751 int iphisecttmp = (kl<2) ? iphisect_dt + kl : iphisect_dt - 1;
752 if (iphisecttmp <0) iphisecttmp = 11;
753 if (iphisecttmp >=12) iphisecttmp = 0;
755 double phipos = iphisecttmp*pival/6.;
756 double phirot = phipos;
769 for (
int ik=1; ik>=0; ik--) {
771 double radial = 407.0;
772 if (ik==0) radial = 382.0;
781 if (steppingHelixstateinfo_.
isValid()) {
784 CLHEP::Hep3Vector hotrkdir2(steppingHelixstateinfo_.
momentum().
x(), steppingHelixstateinfo_.
momentum().
y(),steppingHelixstateinfo_.
momentum().
z());
788 double xx = lclvt0.
x();
789 double yy = lclvt0.y();
792 if ((
std::abs(yy) < 130 && xx >-64.7 && xx <138.2)
795 iphisect = iphisecttmp;
799 if (iphisect != iphisecttmp)
continue;
811 tmpHOCalib.
hoang = CLHEP::Hep3Vector(zLocal.
x(),zLocal.
y(),zLocal.
z()).
dot(hotrkdir2.unit());
826 for (
int i=0;
i<9;
i++) {tmpHOCalib.
hosig[
i]=-100.0;}
827 for (
int i=0;
i<18;
i++) {tmpHOCalib.
hocorsig[
i]=-100.0;}
828 for (
int i=0;
i<9;
i++) {tmpHOCalib.
hbhesig[
i]=-100.0;}
829 tmpHOCalib.
hocro = -100;
830 tmpHOCalib.
htime = -1000;
840 if (iphiho<0) isect -=2000000;
841 tmpHOCalib.
isect = isect;
853 if (
iring==1) {etamn=5; etamx = 10;}
854 if (
iring==2) {etamn=11; etamx = 16;}
855 if (
iring==-1){etamn=-10; etamx = -5;}
856 if (
iring==-2){etamn=-16; etamx = -11;}
861 phimx =2*int((iphiho+1)/2.);
864 phimn = 3*int((iphiho+1)/3.) - 1;
868 if (phimn <1) phimn +=
nphimx;
869 if (phimx >72) phimx -=
nphimx;
880 for (
int i=0;
i<9;
i++) {tmpHOCalib.
hbhesig[
i]=-100.0;}
883 if ((*hbhe).size() >0) {
888 m_coder = (*conditions_).getHcalCoder(
id);
890 int tmpeta=
id.ieta();
891 int tmpphi=
id.iphi();
894 int deta = tmpeta-ietaho;
895 if (tmpeta==-1 && ietaho== 1) deta = -1;
896 if (tmpeta== 1 && ietaho==-1) deta = 1;
897 int dphi = tmpphi-iphiho;
899 if (dphi==71) dphi=-1;
900 if (dphi==-71) dphi=1;
905 if (ipass2 ==0 )
continue;
908 for (
int i=0;
i<(*j).size() &&
i<
nchnmx;
i++) {
913 for (
int i=1;
i<(*j).size() &&
i<=8;
i++) {
918 if (3*(deta+1)+dphi+1<9) tmpHOCalib.
hbhesig[3*(deta+1)+dphi+1] = signal;
929 if ((*hbheht).size()>0) {
930 if(!(*hbheht).size())
throw (
int)(*hbheht).size();
936 int tmpeta=
id.
ieta();
937 int tmpphi=
id.iphi();
939 int deta = tmpeta-ietaho;
940 if (tmpeta==-1 && ietaho== 1) deta = -1;
941 if (tmpeta== 1 && ietaho==-1) deta = 1;
942 int dphi = tmpphi-iphiho;
944 if (dphi==71) dphi=-1;
945 if (dphi==-71) dphi=1;
949 if ( ipass2 ==0 )
continue;
951 float signal = (*j).energy();
953 if (3*(deta+1)+dphi+1<9) tmpHOCalib.
hbhesig[3*(deta+1)+dphi+1] = signal;
962 if ((*ho).size()>0) {
963 int isFilled[netamx*
nphimx];
964 for (
int j=0;
j<netamx*
nphimx;
j++) {isFilled[
j]=0;}
974 m_coder = (*conditions_).getHcalCoder(
id);
976 int tmpeta=
id.ieta();
977 int tmpphi=
id.iphi();
980 if (tmpeta >=etamn && tmpeta <=etamx) {
982 ipass1 = (tmpphi >=phimn && tmpphi <=phimx ) ? 1 : 0;
984 ipass1 = (tmpphi==71 || tmpphi ==72 || tmpphi==1) ? 1 : 0;
988 int deta = tmpeta-ietaho;
989 if (tmpeta==-1 && ietaho== 1) deta = -1;
990 if (tmpeta== 1 && ietaho==-1) deta = 1;
992 int dphi = tmpphi -iphiho;
994 if (dphi==71) dphi=-1;
995 if (dphi==-71) dphi=1;
1000 int tmpeta1 = (tmpeta>0) ? tmpeta -1 : -tmpeta +14;
1002 float tmpdata[
nchnmx]={0,0,0,0,0,0,0,0,0,0};
1003 float sigvall[
nsigpk]={0,0,0,0,0,0,0};
1005 for (
int i=0;
i<(*j).size() &&
i<
nchnmx;
i++) {
1007 if (deta==0 && dphi==0) {
1008 double tmpE = tmpdata[
i] -
pedestal[tmpeta1][tmpphi-1][(*j).sample(
i).capid()];
1016 for (
int ncap=0; ncap<
nsigpk; ncap++) {
1018 sigvall[ncap] +=tmpdata[
i];
1022 if (
i==(*j).size()-1) {
1025 for (
int ij=0; ij<
nsigpk; ij++) {
1026 if (sigvall[ij] > mxled) {mxled = sigvall[ij]; imxled=ij;}
1029 for (
int ij=0; ij<4; ij++) {
1030 pedx +=pedestal[tmpeta1][tmpphi-1][ij];
1032 if (mxled-pedx >2 && mxled-pedx <20 ) {
1034 for (
int jk=0; jk<
ntrgp_gm; jk++) {
1035 if (ntrgpas_gm[jk]>0) {
1036 hopeak[jk]->Fill(nphimx*tmpeta1 + tmpphi-1, imxled+
nstrbn);
1039 if (tmpdata[5]+tmpdata[6] >1) {
1040 horatio->Fill(nphimx*tmpeta1 + tmpphi-1, (tmpdata[5]-tmpdata[6])/(tmpdata[5]+tmpdata[6]));
1042 for (
int ij=0; ij<(*j).size() && ij<
nchnmx; ij++) {
1043 hotime[
ntrgp_gm]->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + ij, tmpdata[ij]);
1044 Nhotime[
ntrgp_gm]->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + ij, 1.);
1045 for (
int jk=0; jk<
ntrgp_gm; jk++) {
1046 if (ntrgpas_gm[jk]>0) {
1047 hotime[jk]->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + ij, tmpdata[ij]);
1048 Nhotime[jk]->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + ij, 1.);
1058 if (
std::abs(tmpeta) <=15 && deta==0 && dphi ==0) {
1061 for (
int i =0;
i<nchnmx &&
i< (*j).size();
i++) {
1062 if (
i >=sigstr &&
i<=sigend)
continue;
1063 signal += tmpdata[
i] - pedestal[tmpeta1][tmpphi-1][(*j).sample(
i).capid()];
1064 if (++icnt >=4)
break;
1066 tmpHOCalib.
hocro = signal;
1070 if (ipass1 ==0 && ipass2 ==0 && cosmicmuon->size()<=2) {
1072 if ((iphiho >=1 && iphiho<=nphimx/2 && tmpphi >=1 && tmpphi <=nphimx/2) ||
1073 (iphiho >nphimx/2 && iphiho<=nphimx && tmpphi >nphimx/2 && tmpphi <=
nphimx)) {
1074 if (isFilled[nphimx*tmpeta1+tmpphi-1]==0) {
1075 isFilled[nphimx*tmpeta1+tmpphi-1]=1;
1076 for (
int i=0;
i<(*j).size() &&
i<
nchnmx;
i++) {
1077 hopedtime->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) +
i, tmpdata[
i]);
1078 Nhopedtime->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + i, 1.);
1079 hopedpr->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + i, tmpdata[i]);
1087 if (ipass1 ==0 && ipass2 ==0 )
continue;
1090 for (
int i=sigstr;
i<(*j).size() &&
i<=sigend;
i++) {
1091 signal += tmpdata[
i] - pedestal[tmpeta1][tmpphi-1][(*j).sample(
i).capid()];
1093 if (signal <-100 || signal >100000) signal = -100;
1095 if (signal >-100 &&
Noccu == Noccu_old) {
1096 for (
int i=0;
i<5;
i++) {
1104 if (ipass1 ==0 && ipass2 ==0 )
continue;
1107 int tmpdph = tmpphi-phimn;
1108 if (tmpdph<0) tmpdph = 2;
1110 int ilog = 2*(tmpeta-etamn)+tmpdph;
1113 ilog = 3*(tmpeta-etamn)+tmpdph;
1115 ilog = 3*(etamx-tmpeta)+tmpdph;
1118 if (ilog>-1 && ilog<18) {
1119 tmpHOCalib.
hocorsig[ilog] = signal;
1124 if (3*(deta+1)+dphi+1<9) tmpHOCalib.
hosig[3*(deta+1)+dphi+1] = signal;
1158 tmpHOCalib.
htime = sumEt/
max(sumE,1.
e-6);
1165 if ((*hoht).size()>0) {
1170 int tmpeta=
id.
ieta();
1171 int tmpphi=
id.iphi();
1174 if (tmpeta >=etamn && tmpeta <=etamx) {
1175 if (phimn < phimx) {
1176 ipass1 = (tmpphi >=phimn && tmpphi <=phimx ) ? 1 : 0;
1178 ipass1 = (tmpphi==71 || tmpphi ==72 || tmpphi==1) ? 1 : 0;
1182 int deta = tmpeta-ietaho;
1183 if (tmpeta==-1 && ietaho== 1) deta = -1;
1184 if (tmpeta== 1 && ietaho==-1) deta = 1;
1186 int dphi = tmpphi -iphiho;
1188 if (dphi==71) dphi=-1;
1189 if (dphi==-71) dphi=1;
1192 float signal = (*j).energy();
1194 int tmpeta1 = (tmpeta>0) ? tmpeta -1 : -tmpeta +14;
1195 if (signal >-100 &&
Noccu == Noccu_old) {
1196 for (
int i=0;
i<5;
i++) {
1206 if (ipass1 ==0 && ipass2 ==0 )
continue;
1209 int tmpdph = tmpphi-phimn;
1210 if (tmpdph<0) tmpdph = 2;
1212 int ilog = 2*(tmpeta-etamn)+tmpdph;
1215 ilog = 3*(tmpeta-etamn)+tmpdph;
1217 ilog = 3*(etamx-tmpeta)+tmpdph;
1220 if (ilog>-1 && ilog<18) {
1221 tmpHOCalib.
hocorsig[ilog] = signal;
1227 if (3*(deta+1)+dphi+1<9) {
1228 tmpHOCalib.
hosig[3*(deta+1)+dphi+1] = signal;
1232 if (deta==0 && dphi ==0) {
1233 tmpHOCalib.
htime = (*j).time();
1234 int crphi = tmpphi + 6;
1235 if (crphi >72) crphi -=72;
1240 int etacr= idcr.
ieta();
1241 int phicr= idcr.
iphi();
1242 if (tmpeta==etacr && crphi ==phicr) {
1255 hostore->push_back(tmpHOCalib);
1263 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
const T & max(const T &a, const T &b)
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
