129 {
131 Gaudi::svcLocator() -> service(
"MessageSvc",
msgSvc);
132 MsgStream log(
msgSvc,
"MdcCalib");
133
134 IDataProviderSvc* eventSvc =
NULL;
135 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
136 SmartDataPtr<MdcDigiCol> mdcDigiCol(eventSvc,"/Event/Digi/MdcDigiCol");
137 if (!mdcDigiCol) {
138 log << MSG::FATAL << "Could not find event" << endreq;
139 }
140
141 int lay;
142 int cel;
143 bool hitCel[43][288];
144 for(lay=0; lay<43; lay++){
145 for(cel=0; cel<288; cel++){
146 hitCel[lay][cel] = false;
147 }
148 }
149
150 MdcDigiCol::iterator
iter = mdcDigiCol->begin();
151 unsigned fgOverFlow;
152 int digiId = 0;
154 for(;
iter != mdcDigiCol->end();
iter++, digiId++) {
156 id = (aDigi)->identify();
159 fgOverFlow = (aDigi) -> getOverflow();
160
161 bool goodTQ = true;
162 if ( ((fgOverFlow & 3) !=0 ) || ((fgOverFlow & 12) != 0) ||
165
166 bool goodT = true;
167 if ( ((fgOverFlow & 1) !=0 ) || (aDigi->
getTimeChannel() == 0x7FFFFFFF) ) goodT =
false;
168
169
170 if(!goodT) continue;
171
172 hitCel[lay][cel] = true;
173 }
174
175 SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc, "/Event/Recon/RecMdcTrackCol");
176 if(!newtrkCol){
177 log << MSG::ERROR << "Could not find RecMdcTrackCol" << endreq;
178 return -1;
179 }
180
181 int nNoiRange = 4;
182 double dphi = 1.0;
183 RecMdcTrackCol::iterator it_trk = newtrkCol->begin();
184 for(it_trk = newtrkCol->begin(); it_trk != newtrkCol->end(); it_trk++){
185 HitRefVec gothits = (*it_trk) -> getVecHits();
186 HitRefVec::iterator it_hit = gothits.begin();
187 HepVector helix = (*it_trk)->helix();
188 HepSymMatrix helixErr = (*it_trk)->err();
189 double phi0Rec = (*it_trk)->helix(1);
190 double delphi;
191 if(phi0Rec>phi0) delphi = phi0Rec - phi0;
192 else delphi = phi0 - phi0Rec;
193 if(delphi > CLHEP::pi) delphi -= CLHEP::pi;
194 if(delphi > (CLHEP::pi*0.17)) continue;
195
196 int nhitLay1 = 0;
197 int nhitLay2 = 0;
198 int nhitLay3 = 0;
199 int nhitLay4 = 0;
200 int nhitT1 = 0;
201 int nhitT2 = 0;
202 int nhitT3 = 0;
203 int nhitT4 = 0;
204 for(lay=0; lay<8; lay++){
205 double docamin = 0.7;
206 if(lay>7) docamin = 0.9;
207 int celmin = -1;
209 for(cel=0; cel<ncel; cel++){
210 double wphi;
211 const MdcGeoWire* pWire = m_mdcGeomSvc -> Wire(lay, cel);
214 double rr = sqrt( (xx * xx) + (yy * yy) );
215 if( yy >= 0 ) wphi = acos(xx / rr);
216 else wphi = CLHEP::twopi - acos(xx / rr);
217
218 double phiTrk = phi0 + CLHEP::halfpi;
219 if(phiTrk > CLHEP::twopi) phiTrk -= CLHEP::twopi;
220 if( !( (fabs(wphi-phiTrk) < dphi) || (fabs(wphi+CLHEP::twopi-phiTrk) < dphi) ||
221 (fabs(wphi-CLHEP::twopi-phiTrk) < dphi) ) ){
222 continue;
223 }
224
225 double doca = m_mdcUtilitySvc->
doca(lay, cel, helix, helixErr);
226
227 if(fabs(doca) < fabs(docamin)){
228 docamin = doca;
229 celmin = cel;
230 }
231 }
232 if(celmin > -1){
233 if(lay<4) nhitLay1++;
234 else if(lay<8) nhitLay2++;
235 else if(lay<20) nhitLay3++;
236 else nhitLay4++;
237 for(int ii=(-nNoiRange); ii<=nNoiRange; ii++){
238 if(0==ii) continue;
239 int icell = celmin + ii;
240 if(icell >= ncel) icell -= ncel;
241 if(icell < 0) icell += ncel;
242
243 if(hitCel[lay][icell]){
244 if(lay<4) nhitT1++;
245 else if(lay<8) nhitT2++;
246 else if(lay<20) nhitT3++;
247 else nhitT4++;
248 }
249 }
250 }
251 }
252 if((nhitLay1<=0) || (nhitLay2<=0) || (nhitLay3<=0) || (nhitLay4<=0)) return false;
253 double ratio1 = (double)nhitT1 / ((double)nhitLay1 * (double)nNoiRange*2.0);
254 double ratio2 = (double)nhitT2 / ((double)nhitLay2 * (double)nNoiRange*2.0);
255 double ratio3 = (double)nhitT3 / ((double)nhitLay3 * (double)nNoiRange*2.0);
256 double ratio4 = (double)nhitT4 / ((double)nhitLay4 * (double)nNoiRange*2.0);
257
258 if((ratio1>0.08) || (ratio2>0.08) || (ratio3>0.03) || (ratio4>0.03)) return false;
259 else return true;
260 }
261 return false;
262}
SmartRefVector< RecMdcHit > HitRefVec
virtual const MdcGeoLayer *const Layer(unsigned id)=0
virtual double doca(int layer, int cell, const HepVector helix, const HepSymMatrix errMat, bool passCellRequired=true, bool doSag=true) const =0
HepPoint3D Backward(void) const
static int layer(const Identifier &id)
Values of different levels (failure returns 0)
static int wire(const Identifier &id)
unsigned int getChargeChannel() const
unsigned int getTimeChannel() const