CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
EcalGeomPhiSymHelper.cc
Go to the documentation of this file.
2 
4 
5 
6 // Geometry
12 
13 //Channel status
14 
17 #include <fstream>
18 
20  const EcalChannelStatus* chStatus,
21  int statusThresold){
22 
23 
24  for (int ieta=0; ieta<kBarlRings; ieta++) nBads_barl[ieta] = 0;
25  for (int ring=0; ring<kEndcEtaRings; ring++) nBads_endc[ring] = 0;
26 
27  for (int ix=0; ix<kEndcWedgesX; ix++) {
28  for (int iy=0; iy<kEndcWedgesY; iy++) {
29 
30  cellPhi_[ix][iy]=0.;
31  cellArea_[ix][iy]=0.;
32  endcapRing_[ix][iy]=-1;
33  }
34  }
35 
36 
37  // loop over all barrel crystals
38  const std::vector<DetId>& barrelCells = geometry->getValidDetIds(DetId::Ecal, EcalBarrel);
39  std::vector<DetId>::const_iterator barrelIt;
40 
41  for (barrelIt=barrelCells.begin(); barrelIt!=barrelCells.end(); barrelIt++) {
42  EBDetId eb(*barrelIt);
43 
44  int sign = eb.zside()>0 ? 1 : 0;
45 
46  int chs= (*chStatus)[*barrelIt].getStatusCode() & 0x001F;
47  if( chs <= statusThresold)
48  goodCell_barl[abs(eb.ieta())-1][eb.iphi()-1][sign] = true;
49 
50  if( !goodCell_barl[abs(eb.ieta())-1][eb.iphi()-1][sign] )
51  nBads_barl[abs(eb.ieta())-1]++;
52 
53  }
54 
55 
56 
57  const CaloSubdetectorGeometry *endcapGeometry =
59 
60  for (int ix=0; ix<kEndcWedgesX; ix++) {
61  for (int iy=0; iy<kEndcWedgesY; iy++) {
62  cellPos_[ix][iy] = GlobalPoint(0.,0.,0.);
63  cellPhi_[ix][iy]=0.;
64  cellArea_[ix][iy]=0.;
65  endcapRing_[ix][iy]=-1;
66  }
67  }
68 
69  const std::vector<DetId>& endcapCells = geometry->getValidDetIds(DetId::Ecal, EcalEndcap);
70  std::vector<DetId>::const_iterator endcapIt;
71  for (endcapIt=endcapCells.begin(); endcapIt!=endcapCells.end(); endcapIt++) {
72 
73  const CaloCellGeometry *cellGeometry = endcapGeometry->getGeometry(*endcapIt);
74  EEDetId ee(*endcapIt);
75  int ix=ee.ix()-1;
76  int iy=ee.iy()-1;
77 
78  int sign = ee.zside()>0 ? 1 : 0;
79 
80  // store all crystal positions
81  cellPos_[ix][iy] = cellGeometry->getPosition();
82  cellPhi_[ix][iy] = cellGeometry->getPosition().phi();
83 
84  // calculate and store eta-phi area for each crystal front face Shoelace formuls
85  const CaloCellGeometry::CornersVec& cellCorners (cellGeometry->getCorners()) ;
86  cellArea_[ix][iy]=0.;
87 
88  for (int i=0; i<4; i++) {
89  int iplus1 = i==3 ? 0 : i+1;
90  cellArea_[ix][iy] +=
91  cellCorners[i].eta()*float(cellCorners[iplus1].phi()) -
92  cellCorners[iplus1].eta()*float(cellCorners[i].phi());
93 
94  }
95 
96 
97  cellArea_[ix][iy] = fabs(cellArea_[ix][iy])/2.;
98 /*
99  const double deltaPhi =
100  (dynamic_cast<const EcalEndcapGeometry*>(endcapGeometry))->deltaPhi(ee);
101 
102  const double deltaEta =
103  (dynamic_cast<const EcalEndcapGeometry*>(endcapGeometry))->deltaEta(ee) ;
104 
105  cellArea_[ix][iy] = deltaEta*deltaPhi;
106 */
107  int chs= (*chStatus)[*endcapIt].getStatusCode() & 0x001F;
108  if( chs <= statusThresold)
109  goodCell_endc[ix][iy][sign] = true;
110 
111 
112  }
113 
114  // get eta boundaries for each endcap ring
115  etaBoundary_[0]=1.479;
116  etaBoundary_[39]=3.; //It was 4. !!!
117  for (int ring=1; ring<kEndcEtaRings; ring++) {
118  double eta_ring_minus1= cellPos_[ring-1][50].eta() ;
119  double eta_ring = cellPos_[ring][50].eta() ;
120  etaBoundary_[ring]=(eta_ring+eta_ring_minus1)/2.;
121  std::cout << "Eta ring " << ring << " : " << eta_ring << std::endl;
122  }
123 
124 
125 
126 
127 
128  // determine to which ring each endcap crystal belongs,
129  // the number of crystals in each ring,
130  // and the mean eta-phi area of the crystals in each ring
131 
132  for (int ring=0; ring<kEndcEtaRings; ring++) {
133  nRing_[ring]=0;
134  meanCellArea_[ring]=0.;
135  for (int ix=0; ix<kEndcWedgesX; ix++) {
136  for (int iy=0; iy<kEndcWedgesY; iy++) {
137  if (fabs(cellPos_[ix][iy].eta())>etaBoundary_[ring] &&
138  fabs(cellPos_[ix][iy].eta())<etaBoundary_[ring+1]) {
139  meanCellArea_[ring]+=cellArea_[ix][iy];
140  endcapRing_[ix][iy]=ring;
141  nRing_[ring]++;
142 
143 
144 
145  for(int sign=0; sign<kSides; sign++){
146  if( !goodCell_endc[ix][iy][sign] )
147  nBads_endc[ring]++;
148  } //sign
149 
150  } //if
151  } //ix
152  } //iy
153 
155 
156 
157  } //ring
158 
159  // fill phi_endc[ip][ring] vector
160  for (int ring=0; ring<kEndcEtaRings; ring++) {
161 
162  for (int i=0; i<kMaxEndciPhi; i++)
163  phi_endc_[i][ring]=0.;
164 
165  float philast=-999.;
166  for (int ip=0; ip<nRing_[ring]; ip++) {
167  float phimin=999.;
168  for (int ix=0; ix<kEndcWedgesX; ix++) {
169  for (int iy=0; iy<kEndcWedgesY; iy++) {
170  if (endcapRing_[ix][iy]==ring) {
171  if (cellPhi_[ix][iy]<phimin && cellPhi_[ix][iy]>philast) {
172  phimin=cellPhi_[ix][iy];
173  } //if edges
174  } //if ring
175  } //iy
176  } //ix
177  phi_endc_[ip][ring]=phimin;
178  philast=phimin;
179  } //ip
180 
181  } //ring
182 
183  // Print out detid->ring association
184  std::fstream eeringsf("endcaprings.dat",std::ios::out);
185  for (endcapIt=endcapCells.begin(); endcapIt!=endcapCells.end();endcapIt++){
186  EEDetId eedet(*endcapIt);
187  eeringsf<< eedet.hashedIndex()<< " "
188  << endcapRing_[eedet.ix()-1][eedet.iy()-1] << " "
189  << cellPhi_ [eedet.ix()-1][eedet.iy()-1] << " "
190  << cellArea_[eedet.ix()-1][eedet.iy()-1]/
191  meanCellArea_[endcapRing_[eedet.ix()-1][eedet.iy()-1]] << std::endl;
192 
193  }
194 }
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
Definition: CaloGeometry.cc:43
int i
Definition: DBlmapReader.cc:9
void setup(const CaloGeometry *geometry, const EcalChannelStatus *chstatus, int statusThreshold)
int ix() const
Definition: EEDetId.h:71
int nRing_[kEndcEtaRings]
int nBads_endc[kEndcEtaRings]
static const int kBarlRings
bool goodCell_barl[kBarlRings][kBarlWedges][kSides]
Geom::Phi< T > phi() const
Definition: PV3DBase.h:68
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
GlobalPoint cellPos_[kEndcWedgesX][kEndcWedgesY]
#define abs(x)
Definition: mlp_lapack.h:159
static const int kMaxEndciPhi
static const int kSides
T eta() const
const GlobalPoint & getPosition() const
Returns the position of reference for this cell.
double phi_endc_[kMaxEndciPhi][kEndcEtaRings]
virtual const CaloCellGeometry * getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
int iphi() const
get the crystal iphi
Definition: EBDetId.h:46
static const int kEndcWedgesX
int endcapRing_[kEndcWedgesX][kEndcWedgesY]
static const int kEndcEtaRings
int zside() const
Definition: EEDetId.h:65
int iy() const
Definition: EEDetId.h:77
int ieta() const
get the crystal ieta
Definition: EBDetId.h:44
tuple out
Definition: dbtoconf.py:99
int hashedIndex() const
Definition: EEDetId.h:177
std::vector< DetId > getValidDetIds() const
Get the list of all valid detector ids.
Definition: CaloGeometry.cc:90
double cellPhi_[kEndcWedgesX][kEndcWedgesY]
T eta() const
Definition: PV3DBase.h:75
ESHandle< TrackerGeometry > geometry
double meanCellArea_[kEndcEtaRings]
double cellArea_[kEndcWedgesX][kEndcWedgesY]
tuple cout
Definition: gather_cfg.py:121
int nBads_barl[kBarlRings]
virtual const CornersVec & getCorners() const =0
Returns the corner points of this cell&#39;s volume.
int zside() const
get the z-side of the crystal (1/-1)
Definition: EBDetId.h:40
bool goodCell_endc[kEndcWedgesX][kEndcWedgesX][kSides]
Definition: DDAxes.h:10
double etaBoundary_[kEndcEtaRings+1]
static const int kEndcWedgesY