121{
122 MsgStream log(
msgSvc(), name());
123 log << MSG::INFO << "DedxCalibEvent::genNtuple()" << endreq;
124
125 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
126 if (!eventHeader)
127 {
128 log << MSG::INFO << "Could not find Event Header" << endreq;
129 return;
130 }
131 int eventNO = eventHeader->eventNumber();
132 int runNO = eventHeader->runNumber();
133 log << MSG::INFO << "now begin the event: runNO= "<<runNO<<" eventNO= "<< eventNO<< endreq;
134
135 int runFlag=0;
136 if(runNO<
RUN0) runFlag =0;
137 if(runNO>=
RUN1 && runNO<
RUN2) runFlag =1;
138 if(runNO>=
RUN2 && runNO<
RUN3) runFlag =2;
139 if(runNO>=
RUN3 && runNO<
RUN4) runFlag =3;
140 if(runNO>=
RUN4 && runNO<
RUN5) runFlag =4;
141 if(runNO>=
RUN5 && runNO<
RUN6) runFlag =5;
142 if(runNO>=
RUN6 && runNO<
RUN7) runFlag =6;
143 if(runNO>=
RUN7) runFlag =7;
144
145 double tes = -999.0;
146 int esTimeflag = -1;
147 SmartDataPtr<RecEsTimeCol> aevtimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
148 if( (!aevtimeCol) || (aevtimeCol->size()==0) ){
149 tes = -9999.0;
150 esTimeflag = -1;
151 }else{
152 RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
153 for(; iter_evt!=aevtimeCol->end(); iter_evt++){
154 tes = (*iter_evt)->getTest();
155 esTimeflag = (*iter_evt)->getStat();
156 }
157 }
158 if(runFlag ==2) {if( tes>1000 ) return;}
159 else if(runFlag ==3 ){if (tes>700 ) return;}
160 else {if (tes>1400 ) return;}
161
162 SmartDataPtr<RecMdcDedxCol> newexCol(eventSvc(),"/Event/Recon/RecMdcDedxCol");
163 if (!newexCol)
164 {
165 log << MSG::FATAL << "Could not find RecMdcDedxCol" << endreq;
166 return;
167 }
168
170 iGood.clear();
171 int nCharge = 0;
172 double db=0,dz=0,pt0=0,charge0=0;
173 for(RecMdcDedxCol::iterator it = newexCol->begin(); it != newexCol->end(); it++)
174 {
175 HepVector a;
176 if((*it)->getMdcKalTrack())
177 {
180 {
183 }
184 else
186 }
187 else if((*it)->getMdcTrack())
188 {
191 }
192 else continue;
193 db = a[0];
194 dz = a[3];
195 pt0 = fabs(1.0/a[2]);
196 charge0 = ( a[2] > 0 )? 1 : -1;
197
198
199 if(fabs(dz) >=
VZ0CUT)
continue;
200 if(db >=
VR0CUT)
continue;
203 iGood.push_back((*it)->trackId());
204 nCharge += charge0;
205 }
206
208 {
209 if(iGood.size()!=2 || nCharge!=0 ) return;
210 }
211
212
213
214 double poca_x=0,poca_y=0,poca_z=0;
215 float sintheta=0,costheta=0,ptrk=0,ptrk_t=0,charge=0,trackId=0;
216 int Nhit=0,usedhit=0,Nhits=0,Nphlisthit=0;
217 double dEdx_meas_hit=0, dEdx_meas=0,dEdx_meas_esat=0,dEdx_meas_norun=0;
218 double dEdx_hit[100]={0},pathlength_hit[100]={0},wid_hit[100]={0},layid_hit[100]={0},dd_in_hit[100]={0},eangle_hit[100]={0},zhit_hit[100]={0};
219 int trk_recalg = -1;
221
222 for(RecMdcDedxCol::iterator it = newexCol->begin(); it != newexCol->end(); it++)
223 {
224
226 for(int i = 0; i < iGood.size(); i++)
227 {
228 if((*it)->trackId()==iGood[i])
flag=
true;
229 }
230 if(
flag==
false)
continue;
231
232 HepVector a;
233 HepSymMatrix tkErrM;
234 if((*it)->getMdcKalTrack())
235 {
236
237 poca_x = (*it)->getMdcKalTrack()->x();
238 poca_y = (*it)->getMdcKalTrack()->y();
239 poca_z = (*it)->getMdcKalTrack()->z();
240
242
244 {
248 }
249 else
250 {
253 }
254 }
255 else if((*it)->getMdcTrack())
256 {
257 poca_x = (*it)->getMdcTrack()->x();
258 poca_y = (*it)->getMdcTrack()->y();
259 poca_z = (*it)->getMdcTrack()->z();
260
264 }
265 else continue;
266
267 sintheta =
sin(M_PI_2 - atan(a[4]));
268 costheta =
cos(M_PI_2 - atan(a[4]));
269 ptrk_t = 1.0/fabs( a[2] );
270 ptrk = ptrk_t*sqrt(1 + a[4]*a[4]);
271 charge = ( a[2] > 0 )? 1 : -1;
272
273 Nhit = (*it)->numTotalHits();
274 Nhits = ((*it)->getVecDedxHits()).size();
275 usedhit = (*it)->numGoodHits();
276 trk_recalg = (*it)->status();
277 trackId = (*it)->trackId();
278
280 {
281 if(runFlag ==3 &&(ptrk>1.88 || ptrk<1.80)) continue;
282 if(runFlag ==4 &&(ptrk>1.72 || ptrk<1.45)) continue;
283 if(runFlag ==5 &&(ptrk>2.00 || ptrk<1.70)) continue;
284 if(runFlag ==6 &&(ptrk>1.90 || ptrk<1.00)) continue;
285 if(runFlag ==7 &&(ptrk>1.90 || ptrk<1.40)) continue;
286 if(
abs(costheta)>0.9)
continue;
287
288 if(Nhit<20) continue;
289 if(usedhit<6) continue;
290 }
291
292
293 int layid=0,localwid=0,w0id=0,wid=0,lr=0;
294 double p_hit=0,adc_raw=0,tdc_raw=0,zhit=0,driftd=0,driftD=0,driftT=0,dd_in=0,dd_ex=0,eangle=0;
295 double ph=0,pathlength=0,rphi_path=0;
296 long k=0;
297
299 DedxHitRefVec::iterator it_gothit = gothits.begin();
300 for(;it_gothit!=gothits.end(); it_gothit++)
301 {
302
303 if((*it_gothit)->isMdcHitValid())
304 {
305 RecMdcHit* itor = (*it_gothit)->getMdcHit();
310 wid = w0id + localwid;
314
316 if(lr == 2) cout<<"lr = "<<lr<<endl;
319 driftd =
abs(driftD);
322 if(lr == 0 || lr == 2 ) {dd_in = -
abs(dd_in);dd_ex = -
abs(dd_ex);}
323 if(lr == 1 ) {dd_in =
abs(dd_in);dd_ex =
abs(dd_ex);}
326 eangle = eangle/
M_PI;
327 }
328 else if((*it_gothit)->isMdcKalHelixSegValid())
329 {
332 p_hit = 1.0/fabs(a_hit_in(3))*sqrt(1+a_hit_in(5)*a_hit_in(5));
333
338 wid = w0id + localwid;
342
344 if(lr == 2) cout<<"lr = "<<lr<<endl;
345 driftD = itor->
getDD();
346 driftd =
abs(driftD);
347 driftT = itor->
getDT();
350 if(lr==-1 || lr == 2) {dd_in = dd_in; dd_ex = dd_ex;}
351 else if(lr ==1) {dd_in = -dd_in; dd_ex = -dd_ex;}
353 eangle = eangle/
M_PI;
354 }
355 else continue;
356
357 pathlength=(*it_gothit)->getPathLength();
358 rphi_path=pathlength*sintheta;
359 ph = (*it_gothit)->getDedx();
360 if(layid>3)
361 {
362 dEdx_hit[k]=adc_raw;
363 pathlength_hit[k]=pathlength;
364 wid_hit[k]=wid;
365 layid_hit[k]=layid;
366 dd_in_hit[k]=dd_in;
367 eangle_hit[k]=eangle;
368 zhit_hit[k]=zhit;
369
370 k++;
371 }
372
373
374 m_charge1 = charge;
375 m_t01 = tes;
376 m_driftT = driftT;
377 m_tdc_raw = tdc_raw;
378 m_phraw = adc_raw;
379 m_exraw = ph;
380 m_localwid = localwid;
381 m_wire = wid;
382 m_runNO1 = runNO;
383 m_runFlag1 = runFlag;
384 m_doca_in = dd_in;
385 m_doca_ex = dd_ex;
386 m_driftdist = driftD;
387 m_eangle = eangle;
388 m_zhit = zhit;
389 m_costheta1 = costheta;
390 m_pathL = pathlength;
391 m_layer = layid;
392 m_ptrk1 = ptrk;
393 m_ptrk_hit = p_hit;
394 m_trackId1 = trackId;
395
396 StatusCode status = m_nt2->write();
397 if ( status.isFailure() )
398 log << MSG::ERROR << "Cannot fill Ntuple n102" << endreq;
399 }
400
401 Nphlisthit = k;
402 dEdx_meas = (*it)->probPH();
403
404 dEdx_meas_esat = (*it)->getDedxEsat();
405 dEdx_meas_norun = (*it)->getDedxNoRun();
406 dEdx_meas_hit = (*it)->getDedxHit();
407
408
409
410 m_poca_x = poca_x;
411 m_poca_y = poca_y;
412 m_poca_z = poca_z;
413 m_ptrk_t=ptrk_t;
414 m_ptrk=ptrk;
415 m_sintheta=sintheta;
416 m_costheta=costheta;
417 m_charge=charge;
418 m_runNO = runNO;
419 m_runFlag = runFlag;
420 m_evtNO = eventNO;
421 m_t0 = tes;
422 m_trackId = trackId;
423 m_recalg = trk_recalg;
424
425 m_nhit=Nhit;
426 m_nhits=Nhits;
427 m_nphlisthit=Nphlisthit;
428 m_usedhit=usedhit;
429 for(int j=0; j<Nphlisthit; j++)
430 {
431 m_dEdx_hit[j]=dEdx_hit[j];
432 m_pathlength_hit[j]=pathlength_hit[j];
433 m_wid_hit[j]=wid_hit[j];
434 m_layid_hit[j]=layid_hit[j];
435 m_dd_in_hit[j]=dd_in_hit[j];
436 m_eangle_hit[j]=eangle_hit[j];
437 m_zhit_hit[j]=zhit_hit[j];
438 }
439
440
441 m_dEdx_meas = dEdx_meas;
442
443
444
445 m_parttype = (*it)->particleId();
446 m_chidedxe=(*it)->chiE();
447 m_chidedxmu=(*it)->chiMu();
448 m_chidedxpi=(*it)->chiPi();
449 m_chidedxk=(*it)->chiK();
450 m_chidedxp=(*it)->chiP();
451 for(int i=0;i<5;i++)
452 {
453 m_probpid[i]= (*it)->getPidProb(i);
454 m_expectid[i]= (*it)->getDedxExpect(i);
455 m_sigmaid[i]= (*it)->getSigmaDedx(i);
456 }
457 StatusCode status = m_nt1->write();
458 if ( status.isFailure() )
459 {
460 log << MSG::ERROR << "Cannot fill Ntuple n103" << endreq;
461 }
462 }
463
464}
double sin(const BesAngle a)
double cos(const BesAngle a)
SmartRefVector< RecMdcDedxHit > DedxHitRefVec
const HepVector & getZHelix(const int pid) const
const HepSymMatrix & getFError(const int pid) const
const HepVector & getFHelix(const int pid) const
const HepSymMatrix err() const
const HepVector helix() const
......
virtual const MdcGeoLayer *const Layer(unsigned id)=0
static int layer(const Identifier &id)
Values of different levels (failure returns 0)
static int wire(const Identifier &id)
const int getFlagLR(void) const
const double getZhit(void) const
const double getAdc(void) const
const Identifier getMdcId(void) const
const double getTdc(void) const
const double getDriftDistRight(void) const
const double getDriftT(void) const
const double getEntra(void) const
const double getDriftDistLeft(void) const
const double getDoca(void) const
Identifier getMdcId(void) const
int getFlagLR(void) const
HepVector & getHelixIncl(void)
double getEntra(void) const
double getDocaIncl(void) const
double getDocaExcl(void) const
double getTdc(void) const
double getZhit(void) const
double getAdc(void) const
const HepVector & getZHelix() const
const HepSymMatrix & getFError() const
const HepVector & getFHelix() const