Loading...
Searching...
No Matches
CgemSegmentFitAlg Class Reference
#include <CgemSegmentFitAlg.h>
Inheritance diagram for CgemSegmentFitAlg:
Public Member Functions | |
| CgemSegmentFitAlg (const std::string &name, ISvcLocator *pSvcLocator) | |
| StatusCode | execute () |
| StatusCode | initialize () |
| StatusCode | finalize () |
| StatusCode | exe_v1 () |
| StatusCode | exe_v2 () |
Detailed Description
Definition at line 31 of file CgemSegmentFitAlg.h.
Constructor & Destructor Documentation
◆ CgemSegmentFitAlg()
| CgemSegmentFitAlg::CgemSegmentFitAlg | ( | const std::string & | name, |
| ISvcLocator * | pSvcLocator ) |
Definition at line 64 of file CgemSegmentFitAlg.cxx.
64 :
65 Algorithm(name, pSvcLocator)
66{
67 declareProperty("debug", m_debug = 0);
68 declareProperty("check", m_checkIdx = 0);
69 declareProperty("version", m_ver = 2);
70}
Member Function Documentation
◆ exe_v1()
| StatusCode CgemSegmentFitAlg::exe_v1 | ( | ) |
Definition at line 164 of file CgemSegmentFitAlg.cxx.
165{
166 MsgStream log(msgSvc(), name());
167 log << MSG::INFO << "in excute()" << endreq;
168
170 n_layer = 3;
171
172 // SmartDataPtr<CgemGeomSvc> cgemGeom(eventSvc(), "CgemGeomSvc");
173 // if (!cgemGeom){
174 // log << MSG::WARNING << "Could not retrieve Cgem Geomtry list" << endreq;
175 // return StatusCode::FAILURE;
176 // }
177 // CgemGeomSvc::iterator itGeom = cgemGeom->begin();
178
179 int vec_clusterid[3];
180 m_seg_map.clear();
181
182 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
183 if (!eventHeader) {
184 log << MSG::FATAL << "Could not find Event Header" << endreq;
185 return StatusCode::FAILURE;
186 }
187 cgem_evt = eventHeader->eventNumber();
188 cgem_run = eventHeader->runNumber();
189
190 // if(cgem_evt<=221) return StatusCode::FAILURE;
191 // cout<<"event_ID = "<<cgem_evt<<endl;
192 //if(cgem_evt%1000==0)cout<<"event_ID = "<<cgem_evt<<endl;
193 // if(cgem_evt > 62)return StatusCode::SUCCESS;
194
195 SmartDataPtr<RecCgemSegmentCol> recCgemSegmentCol(eventSvc(), "/Event/Recon/RecCgemSegmentCol");
196 if (!recCgemSegmentCol){
197 log << MSG::WARNING << "Could not retrieve Cgem segment list" << endreq;
198 check_segNum++;
199 // return StatusCode::FAILURE;
200 return StatusCode::SUCCESS;
201 }
202 RecCgemSegmentCol::iterator itSeg;
203 /*
204 itSeg = recCgemSegmentCol->begin();
205 for(;itSeg != recCgemSegmentCol->end(); itSeg++){
206 cout<<"event+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
207 cout<<"event ID = "<<(*itSeg)->geteventid()<<endl;;
208 cout<<"segment ID = "<<(*itSeg)->getsegmentid()<<endl;;
209 cout<<"cluster ID_1 = "<<(*itSeg)->getclusterid_1()<<endl;;
210 cout<<"cluster ID_2 = "<<(*itSeg)->getclusterid_2()<<endl;;
211 cout<<"cluster ID_3 = "<<(*itSeg)->getclusterid_3()<<endl;;
212 cout<<"event+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
213 }
214 */
215 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
216 if (!recCgemClusterCol){
217 log << MSG::WARNING << "Could not retrieve Cgem cluster list" << endreq;
218 return StatusCode::FAILURE;
219 }
220
221 SmartDataPtr<Event::CgemMcHitCol> recCgemMcHitCol(eventSvc(), "/Event/MC/CgemMcHitCol");
222 if (!recCgemMcHitCol)
223 {
224 log << MSG::WARNING << "Could not find event" << endreq;
225 return StatusCode::FAILURE;
226 }
227
228 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
229 if (!mcParticleCol) {
230 log << MSG::ERROR<< "Could not find McParticle" << endreq;
231 }
232 else{
233
235 McParticleCol::iterator iter_mc = mcParticleCol->begin();
236 for (;iter_mc != mcParticleCol->end(); ++iter_mc ) {
237 int pdg=(*iter_mc)->particleProperty();
238 if(fabs(pdg)!=13) continue;
239 if(truth_ntrack >=10)continue;
240 if(pdg<0)mc_charge = 1.;
241 else mc_charge =-1.;
242 double truth_x = (*iter_mc)->initialPosition().x();
243 double truth_y = (*iter_mc)->initialPosition().y();
244 double truth_z = (*iter_mc)->initialPosition().z();
245
246 double truth_px = (*iter_mc)->initialFourMomentum().px();
247 double truth_py = (*iter_mc)->initialFourMomentum().py();
248 double truth_pz = (*iter_mc)->initialFourMomentum().pz();
249 double truth_P = sqrt(truth_px*truth_px+truth_py*truth_py+truth_pz*truth_pz);
250 truth_CosTheta[truth_ntrack] = truth_pz/truth_P;
251
252 truth_position = HepPoint3D(truth_x, truth_y, truth_z);
253 truth_p = Hep3Vector(truth_px, truth_py, truth_pz);
254
256 HepPoint3D tmp_pivot(0., 0., 0.);
257 tmp_helix->pivot(tmp_pivot);
258 truth_dr[truth_ntrack] = tmp_helix->dr();
259 truth_dz[truth_ntrack] = tmp_helix->dz();
260 truth_phi0[truth_ntrack] = tmp_helix->phi0();
261 truth_kappa[truth_ntrack] = tmp_helix->kappa();
262 truth_tanl[truth_ntrack] = tmp_helix->tanl();
263 truth_charge[truth_ntrack] = mc_charge;
264 if(m_debug){
265 cout<<"Track Helix : "<<endl
267 <<"\ndz = "<<tmp_helix->dz()<<"\ttanl = "<<tmp_helix->tanl()<<"\tpt = "<<1./tmp_helix->kappa()<<endl<<endl;
268 }
269
270 truth_ntrack++;
271 delete tmp_helix;
272 }
273 }
274 bool use_truth_hits = true;
275 if(use_truth_hits){
276 SmartDataPtr<Event::CgemMcHitCol> lv_CgemMcHitCol(eventSvc(), "/Event/MC/CgemMcHitCol");
277 if (!lv_CgemMcHitCol)
278 {
279 log << MSG::WARNING << "Could not find event" << endreq;
280 return StatusCode::FAILURE;
281 }
282 Event::CgemMcHitCol::iterator iter_truth = lv_CgemMcHitCol->begin();
283 for( ; iter_truth != lv_CgemMcHitCol->end(); ++iter_truth){
284 if(fabs((*iter_truth)->GetPDGCode())!=13) continue;
285 double tmp_x = ((*iter_truth)->GetPositionXOfPrePoint()+(*iter_truth)->GetPositionXOfPostPoint())/2;
286 double tmp_y = ((*iter_truth)->GetPositionYOfPrePoint()+(*iter_truth)->GetPositionYOfPostPoint())/2;
287 double tmp_z = ((*iter_truth)->GetPositionZOfPrePoint()+(*iter_truth)->GetPositionZOfPostPoint())/2;
288 double tmp_phi = atan2(tmp_y,tmp_x);
290 mc_phi[(*iter_truth)->GetLayerID()] = tmp_phi;
291 double tmp_v = ((*iter_truth)->GetPositionZOfPostPoint()+(*iter_truth)->GetPositionZOfPrePoint())/2./10.*sin(CgemSegmentFun::a_stero[(*iter_truth)->GetLayerID()])*CgemSegmentFun::r_layer[(*iter_truth)->GetLayerID()]/CgemSegmentFun::rl[(*iter_truth)->GetLayerID()]+CgemSegmentFun::r_layer[(*iter_truth)->GetLayerID()]*tmp_phi*cos(CgemSegmentFun::a_stero[(*iter_truth)->GetLayerID()]);
292 mc_v[(*iter_truth)->GetLayerID()] = tmp_v;
293 CgemSegmentFun::vec_z[(*iter_truth)->GetLayerID()] = ((*iter_truth)->GetPositionZOfPostPoint()+(*iter_truth)->GetPositionZOfPrePoint())/2./10.;//cm
294 CgemSegmentFun::vec_phi[(*iter_truth)->GetLayerID()] = tmp_phi;
295 // cout<<"truth_z = "<<tmp_z<<"\ttruth_v = "<<tmp_v<<"\ttruth_phi = "<<tmp_phi<<endl;
296 // cout<<"truth_z = "<<tmp_z<<"\ttruth_v = "<<tmp_v<<"\ttruth_phi = "<<tmp_phi<<endl;
297 }
298 double Ini_d0(0.), Ini_phi0(0.), Ini_kappa(0.), Ini_z0(0.), Ini_tanl(0.);
299 CgemSegmentFun::GetHelixVarBeforeFit( -1, Ini_d0, Ini_phi0, Ini_kappa, Ini_z0, Ini_tanl);
300 }
301
302 itSeg = recCgemSegmentCol->begin();
303 cgem_nsegment=0;
304 for(;itSeg != recCgemSegmentCol->end(); ++itSeg){
305 if(cgem_nsegment >= 1000){check_segOver++;continue;}
306 cgem_segmentid[cgem_nsegment] = (*itSeg)->getsegmentid();
307 vec_clusterid[0] = (*itSeg)->getclusterid_1();
308 vec_clusterid[1] = (*itSeg)->getclusterid_2();
309 vec_clusterid[2] = (*itSeg)->getclusterid_3();
310 int IsTruth = (*itSeg)->getmatch();
311
312 // initialize minuit
313 Int_t ierflg;
314 Double_t arglist[10];
315 m_failed = 0;
316
317 m_pMnTrk = new TMinuit(5);
318 m_pMnTrk -> SetPrintLevel(-1);
319 m_pMnTrk -> SetFCN(CgemSegmentFun::fcnTrk);
320 m_pMnTrk -> SetErrorDef(1.0); // For Chisq
321 // m_pMnTrk -> SetErrorDef(0.5); // For likelihood
322 m_pMnTrk -> mnparm(0, "dr" , 0, 0.1 , -10, 10, ierflg);//Set starting values and step sizes for parameters
323 m_pMnTrk -> mnparm(1, "phi0" , 0, 0.16, 0, 6.3, ierflg);
324 m_pMnTrk -> mnparm(2, "kappa", 0, 0.4, -100, 100, ierflg);
325 m_pMnTrk -> mnparm(3, "dz" , 0, 0.5 , -50, 50, ierflg);
326 m_pMnTrk -> mnparm(4, "tanL" , 0, 0.005 , -10, 10, ierflg);
327 //m_pMnTrk -> mnparm(0, "dr" , 0, 0.09 , -10, 10, ierflg);//Set starting values and step sizes for parameters
328 //m_pMnTrk -> mnparm(1, "phi0" , 0, 0.016, 0, 6.3, ierflg);
329 //m_pMnTrk -> mnparm(2, "kappa", 0, 0.01, -30, 30, ierflg);
330 //m_pMnTrk -> mnparm(3, "dz" , 0, 0.05 , -50, 50, ierflg);
331 //m_pMnTrk -> mnparm(4, "tanL" , 0, 0.004 , -10, 10, ierflg);
332 arglist[0] = 0;
333 m_pMnTrk -> mnexcm("Set NOW", arglist, 0, ierflg);
334
335
336 RecCgemClusterCol::iterator itTrk=recCgemClusterCol->begin();
337 for(; itTrk != recCgemClusterCol->end(); ++itTrk){
338 if((*itTrk)->getflag() == 0){
339 n_stripX[(*itTrk)->getlayerid()] = (*itTrk)->getclusterflage()-(*itTrk)->getclusterflagb()+1;
340 CgemSegmentFun::nStrip_X[(*itTrk)->getlayerid()] = (*itTrk)->getclusterflage()-(*itTrk)->getclusterflagb()+1;//cm
341 }
342 if((*itTrk)->getflag() == 1){
343 n_stripV[(*itTrk)->getlayerid()] = (*itTrk)->getclusterflage()-(*itTrk)->getclusterflagb()+1;
344 CgemSegmentFun::nStrip_V[(*itTrk)->getlayerid()] = (*itTrk)->getclusterflage()-(*itTrk)->getclusterflagb()+1;
345 }
346 if((*itTrk)->getflag() != 2) continue;
347 if(!CheckClusterID((*itTrk)->getlayerid(), (*itTrk)->getclusterid(), vec_clusterid)) continue;
348 //CgemSegmentFun::sheet_ID[(*itTrk)->getlayerid()] = (*itTrk)->getsheetid();
349 double tmp_phi = (*itTrk)->getrecphi();
350
352 CgemSegmentFun::vec_phi[(*itTrk)->getlayerid()] = (*itTrk)->getrecphi();
353 //CgemSegmentFun::vec_v[(*itTrk)->getlayerid()] = (*itTrk)->getrecv()/10.;//cm
354 rec_phi[(*itTrk)->getlayerid()] = (*itTrk)->getrecphi();
355 rec_v[(*itTrk)->getlayerid()] = (*itTrk)->getRecZ()/10.*sin(CgemSegmentFun::a_stero[(*itTrk)->getlayerid()])*CgemSegmentFun::r_layer[(*itTrk)->getlayerid()]/CgemSegmentFun::rl[(*itTrk)->getlayerid()]+CgemSegmentFun::r_layer[(*itTrk)->getlayerid()]*tmp_phi*cos(CgemSegmentFun::a_stero[(*itTrk)->getlayerid()]);//cm
357 //cout<<"vec_Z = "<<(*itTrk)->getRecZ()/10.<<"\tvec_phi = "<<(*itTrk)->getrecphi()<<endl;
358 //cout<<"layerID = "<<(*itTrk)->getlayerid()<<"\tX_n_strip = "<<n_stripX[(*itTrk)->getlayerid()] <<"\tV_n_strip = "<<n_stripV[(*itTrk)->getlayerid()] <<endl;
359 //cout<<"vec_v["<<(*itTrk)->getlayerid()<<"] = "<<(*itTrk)->getrecv()/10.<<"\t"<<rec_v[(*itTrk)->getlayerid()]<<endl;
360 //cout<<"sheetID["<<(*itTrk)->getlayerid()<<"] = "<<(*itTrk)->getsheetid()<<endl;
361 }
362
363 double Ini_d0(0.), Ini_phi0(0.), Ini_kappa(0.), Ini_z0(0.), Ini_tanl(0.);
364 double chisq(0.);
365
366 int charge = GetChargeOfTrack();
367 cgem_charge[cgem_nsegment] = charge;
368 CgemSegmentFun::GetHelixVarBeforeFit( charge, Ini_d0, Ini_phi0, Ini_kappa, Ini_z0, Ini_tanl);
369
370 double ErrMatrix[5][5] ;
371 double TrkPar[5] = {Ini_d0, Ini_phi0, Ini_kappa, Ini_z0, Ini_tanl};
372 //double TrkPar[5] = {truth_dr[0], truth_phi0[0], truth_kappa[0], truth_dz[0], truth_tanl[0]};
373 double TrkParErr[5] ;
374 double helixErr[15];
375 int fit = mnTrkFit(TrkPar, TrkParErr, &ErrMatrix[0][0], chisq);
376 int fitAgain = 0; double kappa_err;
377 if(!fit){
378 // cout<< "fit failed"<<"\tevent ID = "<<cgem_evt<<"\trun = "<<cgem_run << endl;
379 m_failed = 1;
380 double tmp_recphi[3], tmp_kappa[8];
381 for(int ii = 0, mm = 0; ii < 2; ii++)
382 for(int jj = 0; jj < 2; jj++)
383 for(int kk = 0; kk < 2; kk++, mm++){
384 tmp_recphi[0] = CgemSegmentFun::vec_phi[0] + pow(-1.,ii)*(wide_delta*n_stripX[0]+pitch*(n_stripX[0]-1))/CgemSegmentFun::r_layer[0];
385 tmp_recphi[1] = CgemSegmentFun::vec_phi[1] + pow(-1.,jj)*(wide_delta*n_stripX[1]+pitch*(n_stripX[1]-1))/CgemSegmentFun::r_layer[1];
386 tmp_recphi[2] = CgemSegmentFun::vec_phi[2] + pow(-1.,kk)*(wide_delta*n_stripX[2]+pitch*(n_stripX[2]-1))/CgemSegmentFun::r_layer[2];
387 charge = GetChargeOfTrack(tmp_recphi);
388 tmp_kappa[mm] = CgemSegmentFun::GetKappaAfterFit(charge, tmp_recphi);
389 //cout<<"tmp kappa["<<mm<<"] = "<<tmp_kappa[mm]<<endl;
390 }
391 double kappa_s = GetMin(tmp_kappa);
392 double kappa_l = GetMax(tmp_kappa);
393 kappa_err = (kappa_s - kappa_l)/2.;
394 double TrkPar[5] = {Ini_d0, Ini_phi0, (kappa_l+kappa_s)/2., Ini_z0, Ini_tanl};
395 TrkPar[2] = kappa_l;
396 fitAgain = mnTrkFit(TrkPar, TrkParErr, &ErrMatrix[0][0], chisq);
397
398 }
399 //cout<<"m_failed = "<<m_failed<<endl;
400
401 //cout<<"charge = "<<charge<<"\t"<<fit<<"\t"<<fitAgain<<endl<<endl;
402 // cout<<
403 // "delta_dr = "<<Ini_d0-truth_dr[0]<<endl<<
404 // "delta_phi0 = "<<Ini_phi0-truth_phi0[0]<<endl<<
405 // "delta_kappa = "<<Ini_kappa-truth_kappa[0]<<endl<<
406 // "delta_dz = "<<Ini_z0-truth_dz[0]<<endl<<
407 // "delta_tanl = "<<Ini_tanl-truth_tanl[0]<<endl<<endl;
408
409 //RecCgemSegment* recsegment = new RecCgemSegment();
410 cgem_fitFlag[cgem_nsegment] = fit;
411 if(fit==1||fitAgain){
412 if(fit){
413 for (int ie = 0, k = 0 ; ie < 5 ; ie++){
414 for (int je = ie ; je < 5 ; je ++, k++) {
415 helixErr[k] = ErrMatrix[ie][je];
416 }
417 }
418 }
419
420 if(fitAgain){
421 for (int ie = 0, k = 0 ; ie < 5 ; ie++){
422 for (int je = ie ; je < 5 ; je ++, k++) {
423 if(ie==2&&je<2)helixErr[k] = ErrMatrix[ie+2][je]; //k = 2,6
424 else if(je<2&&ie>2) helixErr[k] = ErrMatrix[ie-1][je]; //k = 3,4,7,8
425 else if(je>2&&ie>2) helixErr[k] = ErrMatrix[ie-1][je-1]; //k = 12,13,14
426 else if(je==2) helixErr[k] = 0.; //k = 9,10,11
427 else helixErr[k] = ErrMatrix[ie][je]; // k = 0,1,5
428
429 if(k==9)helixErr[k] = pow(kappa_err,2);
430 // cout<<"ErrMatrix["<<ie<<"]["<<je<<"] = "<<ErrMatrix[ie][je]<<"\t"<<helixErr[k]<<endl;
431 }
432 }
433 TrkParErr[0] = 0.8812/0.0677*TrkParErr[0];
434 TrkParErr[1] = 0.266/0.004*TrkParErr[1];
435 TrkParErr[2] = kappa_err;
436 helixErr[0] = 0.8812/0.0677*helixErr[0];
437 helixErr[5] = 0.266/0.004*helixErr[5];
438 }
439 // cout<<"helixErr[0] = "<<helixErr[0]<<"\tTrkParErr[0] = "<<TrkParErr[0]<<"\n"<<
440 // "helixErr[5] = "<<helixErr[5]<<"\tTrkParErr[1] = "<<TrkParErr[1]<<"\n"<<
441 // "helixErr[9] = "<<helixErr[9]<<"\tTrkParErr[2] = "<<TrkParErr[2]<<"\n"<<
442 // "helixErr[12] = "<<helixErr[12]<<"\tTrkParErr[3] = "<<TrkParErr[3]<<"\n"<<
443 // "helixErr[14] = "<<helixErr[14]<<"\tTrkParErr[4] = "<<TrkParErr[4]<<"\n"<<endl<<endl;
444 //recsegment->sethelix(TrkPar);
445 //recsegment->sethelix_err(helixErr);
446 //recsegment->setmatch(IsTruth);
447 (*itSeg)->sethelix(TrkPar);
448 (*itSeg)->sethelix_err(helixErr);
449 (*itSeg)->setmatch(IsTruth);
450 cgem_dr[cgem_nsegment] = TrkPar[0];
451 cgem_phi0[cgem_nsegment] = TrkPar[1];
452 cgem_kappa[cgem_nsegment] = TrkPar[2];
453 cgem_dz[cgem_nsegment] = TrkPar[3];
454 cgem_tanl[cgem_nsegment] = TrkPar[4];
455 cgem_chisq[cgem_nsegment] = chisq;
456 if(m_debug){
457 cout<< "fit success" << endl
458 <<"dr = " <<TrkPar[0]<<" +- "<<TrkParErr[0]<<"\t\t"
459 <<"phi0 = "<<TrkPar[1]<<" +- "<<TrkParErr[1]<<"\t\t"
460 <<"kappa= "<<TrkPar[2]<<" +- "<<TrkParErr[2]<<"\n"
461 <<"dz = " <<TrkPar[3]<<" +- "<<TrkParErr[3]<<"\t\t"
462 <<"tanl = "<<TrkPar[4]<<" +- "<<TrkParErr[4]<<"\n\n"
463 <<"kappa_err= "<<TrkParErr[2]<<endl;
464 cout<<"chisq = "<<chisq<<endl;
465 }
466 }
467 else{
468 //cout<< "fit failed"<<"\tevent ID = "<<cgem_evt<<"\trun = "<<cgem_run << endl;
469
470 double fail_helix[5], fail_helix_err[15] ;
471 for(int i = 0; i <15; i++){
472 if(i<5)fail_helix[i] = -999.;
473 fail_helix_err[i] = -999.;
474 }
475 (*itSeg)->sethelix(fail_helix);
476 (*itSeg)->sethelix_err(fail_helix_err);
477 (*itSeg)->setmatch(IsTruth);
478 cgem_dr[cgem_nsegment] = -999.;
479 cgem_phi0[cgem_nsegment] = -999.;
480 cgem_kappa[cgem_nsegment] = -999.;
481 cgem_dz[cgem_nsegment] = -999.;
482 cgem_tanl[cgem_nsegment] = -999.;
483 cgem_chisq[cgem_nsegment] = -999.;
484 }
485 //m_seg_map.insert(pair<int,RecCgemSegment*>(cgem_nsegment,recsegment));
486
487 cgem_nsegment++;
488
489 //cout<<"***************************#######*******************************"<<endl<<endl;
490 delete m_pMnTrk;
491 //delete recsegment;
492 }
493 /*
494 RecCgemSegmentCol* segmentcol = new RecCgemSegmentCol;
495
496 IDataProviderSvc* evtSvc = NULL;
497 Gaudi::svcLocator()->service("EventDataSvc",evtSvc);
498 if (evtSvc) {
499 log << MSG::INFO << "makeTds:event Svc has been found" << endreq;
500 } else {
501 log << MSG::FATAL << "makeTds:Could not find eventSvc" << endreq;
502 return StatusCode::SUCCESS;
503 }
504
505 StatusCode segmentsc;
506 IDataManagerSvc *dataManSvc;
507 dataManSvc= dynamic_cast<IDataManagerSvc*>(evtSvc);
508 DataObject *aSegmentCol;
509 evtSvc->findObject("/Event/Recon/RecCgemSegmentCol",aSegmentCol);
510 if(aSegmentCol != NULL) {
511 dataManSvc->clearSubTree("/Event/Recon/RecCgemSegmentCol");
512 evtSvc->unregisterObject("/Event/Recon/RecCgemSegmentCol");
513 }
514
515
516 //push back the segment from the map//
517 for(m_seg_map_it = m_seg_map.begin();m_seg_map_it!=m_seg_map.end();++m_seg_map_it){
518 RecCgemSegment* recsegment = (*m_seg_map_it).second;
519 segmentcol->push_back(recsegment);
520 cout<<">>>>>>> push recsegment, clusterid_1="<<recsegment->getclusterid_1()<<"kappa: "<<recsegment->gethelix(2)<<endl;
521
522 }
523
524 segmentsc = evtSvc->registerObject("/Event/Recon/RecCgemSegmentCol", segmentcol);
525 if( segmentsc.isFailure() ) {
526 log << MSG::FATAL << "Could not register RecCgemSegment" << endreq;
527 return StatusCode::SUCCESS;
528 }
529 log << MSG::INFO << "RecCgemSegmentCol registered successfully!" <<endreq;
530 */
531
532
533
534 helix_ntuple->write();
535 /*
536 SmartDataPtr<RecCgemSegmentCol> recCgemSegmentCol2(eventSvc(), "/Event/Recon/RecCgemSegmentCol");
537 itSeg = recCgemSegmentCol2->begin();
538 for(;itSeg != recCgemSegmentCol2->end(); itSeg++){
539 cout<<"event+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
540 cout<<"event ID = "<<(*itSeg)->geteventid()<<endl;;
541 cout<<"segment ID = "<<(*itSeg)->getsegmentid()<<endl;;
542 cout<<"cluster ID_1 = "<<(*itSeg)->getclusterid_1()<<endl;;
543 cout<<"cluster ID_2 = "<<(*itSeg)->getclusterid_2()<<endl;;
544 cout<<"cluster ID_3 = "<<(*itSeg)->getclusterid_3()<<endl;;
545 cout<<"event+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
546 }
547 */
548
549
550 return StatusCode::SUCCESS;
551}
Helix parameter class.
double phi0(void) const
double tanl(void) const
double dz(void) const
double dr(void) const
returns an element of parameters.
const HepPoint3D & pivot(void) const
returns pivot position.
double kappa(void) const
void GetHelixVarBeforeFit(int charge, double &d0, double &phi0, double &omega, double &z0, double &tanl)
Definition CgemSegmentFun.cxx:255
double GetKappaAfterFit(int charge, double *rec_phi)
Definition CgemSegmentFun.cxx:240
void fcnTrk(int &npar, double *gin, double &f, double *par, int iflag)
Definition CgemSegmentFun.cxx:205
Referenced by execute().
◆ exe_v2()
| StatusCode CgemSegmentFitAlg::exe_v2 | ( | ) |
Definition at line 553 of file CgemSegmentFitAlg.cxx.
554{
555 MsgStream log(msgSvc(), name());
556 log << MSG::INFO << "in excute()" << endreq;
557
558 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
559 if (!eventHeader) {
560 log << MSG::FATAL << "Could not find Event Header" << endreq;
561 return StatusCode::FAILURE;
562 }
563 m_run = eventHeader->runNumber();
564 m_evt = eventHeader->eventNumber();
565 getMcPar();
566
567 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
568 if (!recCgemClusterCol){
569 log << MSG::WARNING << "Could not retrieve Cgem cluster list" << endreq;
570 return StatusCode::FAILURE;
571 }
572 RecCgemClusterCol::iterator itCluBegin=recCgemClusterCol->begin();
573
574 SmartDataPtr<RecCgemSegmentCol> recCgemSegmentCol(eventSvc(), "/Event/Recon/RecCgemSegmentCol");
575 if (!recCgemSegmentCol){
576 log << MSG::WARNING << "Could not retrieve Cgem segment list" << endreq;
577 check_segNum++;
578 // return StatusCode::FAILURE;
579 return StatusCode::SUCCESS;
580 }
581 if(m_checkIdx) {
582 cgem_run = m_run;
583 cgem_evt = m_evt;
584 cgem_nsegment=0;
585 }
586 int cluId[3];
587 RecCgemSegmentCol::iterator itSeg = recCgemSegmentCol->begin();
588 for( ;itSeg != recCgemSegmentCol->end(); ++itSeg)
589 {
590 cluId[0] = (*itSeg)->getclusterid_1();
591 cluId[1] = (*itSeg)->getclusterid_2();
592 cluId[2] = (*itSeg)->getclusterid_3();
593 if(m_debug)
594 {
595 cout<<"input clusters: "<<endl;
596 }
597 double vec_Q[3], vec_T[3];
598 for(int i=0; i<3; i++)
599 {
600 RecCgemClusterCol::iterator itClu = itCluBegin+cluId[i];
603 vec_Q[i]=(*itClu)->getenergydeposit();
604 vec_T[i]=100;
605 if(m_debug)
606 {
607 cout <<"layer "<<i<<", cluster id = "<<cluId[i]<<", itCluId = "<<(*itClu)->getclusterid()
609 }
610 }
611 //cout<<"charge = "<<GetChargeOfTrack()<<", "<<estiChargeOfSeg()<<endl;
612 int charge = estiChargeOfSeg();
613
614 // initial Helix parameters
615 double dr_ini, phi0_ini, kappa_ini, dZ_ini, tanL_ini;
616 CgemSegmentFun::GetHelixVarBeforeFit(charge, dr_ini, phi0_ini, kappa_ini, dZ_ini, tanL_ini);
617 double dr_ini0 = dr_ini;
618 double phi0_ini0 = phi0_ini;
619 double kappa_ini0 = kappa_ini;
620 double dZ_ini0 = dZ_ini;
621 double tanL_ini0 = tanL_ini;
622
623 double helix[5] = {dr_ini, phi0_ini, kappa_ini, dZ_ini, tanL_ini};
624 double helixErr[15];
625 HepPoint3D pivot(0., 0., 0.);
626 HepVector vec_trkpar(5);
627 for(int i = 0; i < 5; i++){
628 vec_trkpar[i] = helix[i];
629 }
630 //KalmanFit::Helix* cgem_helix = new KalmanFit::Helix(pivot, vec_trkpar);
631 KalmanFit::Helix cgem_helix(pivot, vec_trkpar);
632
633 /* initialize error matrix */
634 int iView(0), mode(2);
635 for(int i=0; i<3; i++)
636 {
638 HepPoint3D pos = cgem_helix.x(dphi);
639 double phi_pos = pos.phi();
640 Hep3Vector p3 = cgem_helix.momentum(dphi);
641 double phi_p = p3.phi();
642 double dPhi = phi_p-phi_pos;
643 while(dPhi>CLHEP::pi) dPhi-=CLHEP::twopi;
644 while(dPhi<-CLHEP::pi) dPhi+=CLHEP::twopi;
645 CgemSegmentFun::x_reso[i]= (myCgemCalibSvc->getSigma(i,0,mode,dPhi,vec_Q[i],vec_T[i]))/10.;// mm -> cm
646 CgemSegmentFun::v_reso[i]= (myCgemCalibSvc->getSigma(i,1,mode,dPhi,vec_Q[i],vec_T[i]))/10.;// mm -> cm
647 }
648
649 // initialize minuit
650 Int_t ierflg, istat;
651 Double_t arglist[10];
652 m_failed = 0;
653
654 m_pMnTrk = new TMinuit(5);// num of par
655 m_pMnTrk -> SetPrintLevel(-1);// no print out
656 m_pMnTrk -> mnexcm("Set NOW", arglist, 1, ierflg);// no warning
657 m_pMnTrk -> SetFCN(CgemSegmentFun::fcnTrk);
658 m_pMnTrk -> SetErrorDef(1.0); // For Chisq
659 // m_pMnTrk -> SetErrorDef(0.5); // For likelihood
660 bool fitSucceed = false;
661 bool fitUsable = false;
662 int istat_last = 0;
663 int nFitTried = 0; int iFitStored = 0;
664 double chi2_min = 999;
665 double chi2_saved = 999;
666 double chi2_fixKapa_min = 990;
667 while((nFitTried<3&&!fitSucceed)||(nFitTried<5&&!fitUsable)) {
668 m_pMnTrk-> mnparm(0, "dr" , dr_ini, 0.001 , dr_ini-10, dr_ini+10, ierflg);//Set starting values and step sizes for parameters
669 m_pMnTrk-> mnparm(1, "phi0" , phi0_ini, 0.001, phi0_ini-1.0, phi0_ini+1.0, ierflg);
670 m_pMnTrk-> mnparm(2, "kappa", kappa_ini, max(fabs(kappa_ini)*0.01,0.01), kappa_ini-5.0, kappa_ini+5.0,ierflg);
671 if(nFitTried>=3) m_pMnTrk->FixParameter(2);
672 m_pMnTrk -> mnparm(3, "dz" , dZ_ini, 0.01 , dZ_ini-20.0, dZ_ini+20.0, ierflg);
673 m_pMnTrk -> mnparm(4, "tanL" , tanL_ini, 0.001 , tanL_ini-1.0, tanL_ini+1.0, ierflg);
674 arglist[0] = 2000; //arglist[1] = 0.1;
675 m_pMnTrk -> mnexcm("MIGRAD", arglist, 1, ierflg);
676 nFitTried++;
677 double trkpar[5], trkparErr[5], emat[25], chisq;
678 Double_t fmin, edm, errdef;
679 Int_t nvpar, nparx;
680 m_pMnTrk -> mnstat(fmin, edm, errdef, nvpar, nparx, istat);
681 for(int i=0; i<5; i++){
682 m_pMnTrk -> GetParameter(i, trkpar[i], trkparErr[i]);
683 }
684 // error matrix
685 m_pMnTrk -> mnemat(emat, 5);
686 chisq = fmin;
687
688 //HepVector aNew = CgemSegmentFun::CalculateHelix(vec_trkpar, pivot);
689
690 if(m_debug) {
691 cout<<"Fit ierflg="<<ierflg<<", istat="<<istat<<endl;
692 cout.precision(4);
693 cout<<setw(10)<<" "<<setw(10)<<"dr"<<setw(10)<<"phi0"<<setw(10)<<"kappa"<<setw(10)<<"dz"<<setw(10)<<"tanL"<<endl;
694 cout<<setw(10)<<"MC truth"<<setw(10)<<m_helix_mc[0]<<setw(10)<<m_helix_mc[1]<<setw(10)<<m_helix_mc[2]<<setw(10)<<m_helix_mc[3]<<setw(10)<<m_helix_mc[4]<<endl;
695 cout<<setw(10)<<"ini"<<setw(10)<<dr_ini<<setw(10)<<phi0_ini<<setw(10)<<kappa_ini<<setw(10)<<dZ_ini<<setw(10)<<tanL_ini<<endl;
696 //cout<<setw(10)<<"calcu"<<setw(10)<<aNew[0]<<setw(10)<<aNew[1]<<setw(10)<<aNew[2]<<setw(10)<<aNew[3]<<setw(10)<<aNew[4]<<"chi2="<<chisq<<endl;
697 cout<<setw(10)<<"fit"<<setw(10)<<trkpar[0]<<setw(10)<<trkpar[1]<<setw(10)<<trkpar[2]<<setw(10)<<trkpar[3]<<setw(10)<<trkpar[4]<<"chi2="<<chisq<<endl;
698 cout<<setw(10)<<"fit err"<<setw(10)<<trkparErr[0]<<setw(10)<<trkparErr[1]<<setw(10)<<trkparErr[2]<<setw(10)<<trkparErr[3]<<setw(10)<<trkparErr[4]<<endl;
699 cout.precision(6);
700 }
701
702 /*
703 *-* IERFLG is now (94.5) defined the same as ICONDN in MNCOMD
704 *-* = 0: command executed normally
705 *-* 1: command is blank, ignored
706 *-* 2: command line unreadable, ignored
707 *-* 3: unknown command, ignored
708 *-* 4: abnormal termination (e.g., MIGRAD not converged)
709 *-* 9: reserved
710 *-* 10: END command
711 *-* 11: EXIT or STOP command
712 *-* 12: RETURN command
713 *
714 * ISTAT: a status integer indicating how good is the covariance
715 *-* matrix: 0= not calculated at all
716 *-* 1= approximation only, not accurate
717 *-* 2= full matrix, but forced positive-definite
718 *-* 3= full accurate covariance matrix
719 */
720 if( (0==ierflg) && (istat==3) ){
721 fitSucceed=true;
722 }
723 else {
724 //cout<<"TMinuit failed"<<endl;
725 if(m_debug) {
726 cout<<"TMinuit tried "<<nFitTried<<" times!"<<endl;
727 }
728 dr_ini=trkpar[0];
729 phi0_ini=trkpar[1];
730 //kappa_ini=trkpar[2]+pow(-1,nFitTried);
731 kappa_ini=trkpar[2]+0.5*pow(-1,nFitTried);
732 dZ_ini=trkpar[3];
733 tanL_ini=trkpar[4];
734 }
735 if(nFitTried>=3) kappa_ini=0.5*pow(-1,nFitTried);
736
737 if( (0==ierflg) && (istat>0) ){
738 if(nFitTried<4) fitUsable = true;
739 if((nFitTried<4&&istat>istat_last)||(nFitTried>3&&chisq<chi2_fixKapa_min)) {
740 int k=0;
741 chi2_fixKapa_min = chisq;
742 chi2_saved = chisq;
743 for(int ii=0; ii<5; ii++)
744 {
745 helix[ii]=trkpar[ii];
746 for(int jj=ii; jj<5; jj++)
747 {
748 helixErr[k]=emat[ii*5+jj];
749 k++;
750 }
751 }
752 istat_last = istat;
753 iFitStored = nFitTried;
754
755 }
756 }
757 }
758 //if(!fitUsable) cout<<"TMinuit failed after "<<nFitTried<<" attemps!"<<endl;
759 if(m_debug)
760 cout<<"Fit results "<<iFitStored<<" are saved!"<<endl;
761
762 if(m_checkIdx&&cgem_nsegment<1000) {
763 cgem_segmentid[cgem_nsegment]=(*itSeg)->getsegmentid();
764 cgem_dr_ini[cgem_nsegment]=dr_ini0;
765 cgem_dz_ini[cgem_nsegment]=dZ_ini0;
766 cgem_phi0_ini[cgem_nsegment]=phi0_ini0;
767 cgem_kappa_ini[cgem_nsegment]=kappa_ini0;
768 cgem_tanl_ini[cgem_nsegment]=tanL_ini0;
769 cgem_dr[cgem_nsegment]=helix[0];
770 cgem_dz[cgem_nsegment]=helix[3];
771 cgem_phi0[cgem_nsegment]=helix[1];
772 cgem_kappa[cgem_nsegment]=helix[2];
773 cgem_tanl[cgem_nsegment]=helix[4];
774 cgem_chisq[cgem_nsegment]=chi2_saved;
775 cgem_fitFlag[cgem_nsegment]=iFitStored;
776 cgem_charge[cgem_nsegment]=charge;
777 cgem_nsegment++;
778 }
779
780 // fill out
781 (*itSeg)->sethelix(helix);
782 (*itSeg)->sethelix_err(helixErr);
783
784 delete m_pMnTrk;
785
786 }// loop CGEM segments
787 //cout<<"end exe_v2()"<<endl;
788 if(m_checkIdx) helix_ntuple->write();
789
790 return StatusCode::SUCCESS;
791
792}
virtual double getSigma(int layer, int xvFlag, int readoutMode, double angle, double Q, double T) const =0
double IntersectCylinder(double r, KalmanFit::Helix helix)
Definition CgemSegmentFun.cxx:225
Referenced by execute().
◆ execute()
| StatusCode CgemSegmentFitAlg::execute | ( | ) |
Definition at line 156 of file CgemSegmentFitAlg.cxx.
157{
160
161 return StatusCode::SUCCESS;
162}
◆ finalize()
| StatusCode CgemSegmentFitAlg::finalize | ( | ) |
Definition at line 1000 of file CgemSegmentFitAlg.cxx.
1000 {
1001
1002 MsgStream log(msgSvc(), name());
1003 log << MSG::INFO << "in finalize()" << endreq;
1004 //cout<<"leave segment number by dZ dPhi = "<<check_segNum<<endl
1005 // <<"leave segment number over 1000 = "<<check_segOver<<endl;
1006
1007 //cout<<endl<<"******* CgemSegmentFitAlg *******"<<endl;
1008
1009 return StatusCode::SUCCESS;
1010}
◆ initialize()
| StatusCode CgemSegmentFitAlg::initialize | ( | ) |
Definition at line 73 of file CgemSegmentFitAlg.cxx.
73 {
74
75 MsgStream log(msgSvc(), name());
76 log << MSG::INFO << "in initialize()" << endreq;
77 if(m_checkIdx>0) {
78 /////////////////////////////// book ntuple //////////////////////////////////////
79 INTupleSvc* ntupleSvc;
82 if ( mt ) helix_ntuple = mt;
83 else{
85 StatusCode status;
86 if(helix_ntuple){
87 status = helix_ntuple->addItem("run", cgem_run);
88 status = helix_ntuple->addItem("evt", cgem_evt);
89 status = helix_ntuple->addItem("ntrack",truth_ntrack,0,10);
90 status = helix_ntuple->addIndexedItem("truth_dr", truth_ntrack , truth_dr);
91 status = helix_ntuple->addIndexedItem("truth_dz", truth_ntrack , truth_dz);
92 status = helix_ntuple->addIndexedItem("truth_phi0", truth_ntrack , truth_phi0);
93 status = helix_ntuple->addIndexedItem("truth_kappa", truth_ntrack , truth_kappa);
94 status = helix_ntuple->addIndexedItem("truth_tanl", truth_ntrack , truth_tanl);
95 status = helix_ntuple->addIndexedItem("truth_charge", truth_ntrack , truth_charge);
96 status = helix_ntuple->addIndexedItem("truth_CosTheta", truth_ntrack , truth_CosTheta);
97 status = helix_ntuple->addItem("nsegment",cgem_nsegment,0,1000);
98 status = helix_ntuple->addIndexedItem("segmentid", cgem_nsegment , cgem_segmentid);
99 status = helix_ntuple->addIndexedItem("dr_ini", cgem_nsegment , cgem_dr_ini);
100 status = helix_ntuple->addIndexedItem("dz_ini", cgem_nsegment , cgem_dz_ini);
101 status = helix_ntuple->addIndexedItem("phi0_ini", cgem_nsegment , cgem_phi0_ini);
102 status = helix_ntuple->addIndexedItem("kappa_ini", cgem_nsegment , cgem_kappa_ini);
103 status = helix_ntuple->addIndexedItem("tanl_ini", cgem_nsegment , cgem_tanl_ini);
104 status = helix_ntuple->addIndexedItem("dr", cgem_nsegment , cgem_dr);
105 status = helix_ntuple->addIndexedItem("dz", cgem_nsegment , cgem_dz);
106 status = helix_ntuple->addIndexedItem("phi0", cgem_nsegment , cgem_phi0);
107 status = helix_ntuple->addIndexedItem("kappa", cgem_nsegment , cgem_kappa);
108 status = helix_ntuple->addIndexedItem("tanl", cgem_nsegment , cgem_tanl);
109 status = helix_ntuple->addIndexedItem("chisq", cgem_nsegment , cgem_chisq);
110 status = helix_ntuple->addIndexedItem("fitFlag", cgem_nsegment , cgem_fitFlag);
111 status = helix_ntuple->addIndexedItem("charge", cgem_nsegment , cgem_charge);
112
113 status = helix_ntuple->addItem("n_layer",n_layer,0,3);
114 status = helix_ntuple->addIndexedItem("n_stripX", n_layer , n_stripX);
115 status = helix_ntuple->addIndexedItem("n_stripV", n_layer , n_stripV);
116 status = helix_ntuple->addIndexedItem("mc_phi", n_layer , mc_phi);
117 status = helix_ntuple->addIndexedItem("mc_v", n_layer , mc_v);
118 status = helix_ntuple->addIndexedItem("rec_phi", n_layer , rec_phi);
119 status = helix_ntuple->addIndexedItem("rec_v", n_layer , rec_v);
120 }
121 }
122 }
123 check_segNum = 0;
124 check_segOver = 0;
125 wide_delta = 570./2/10000; //cm
126 pitch = (650.-570.)/2/10000; //cm
127
128 ISvcLocator* svcLocator = Gaudi::svcLocator();
129 ICgemGeomSvc* ISvc;
130 StatusCode sc=svcLocator->service("CgemGeomSvc", ISvc);
133 for(int i=0; i<nCgemLayer; i++)
134 {
141 if(m_debug) {
142 cout<<"iLayer, r, ang_stereo, Rx, Rv = "<<i<<", "<<CgemSegmentFun::rl[i]<<", "<<CgemSegmentFun::a_stero[i]<<", "<<CgemSegmentFun::r_X[i]<<", "<<CgemSegmentFun::r_V[i]<<endl;
143 }
144 }
145
146 sc = service ("CgemCalibFunSvc", myCgemCalibSvc);
147 if ( sc.isFailure() ){
148 cout<< name() << "Could not load MdcCalibFunSvc!" << endl;
149 return sc;
150 }
151
152 return StatusCode::SUCCESS;
153
154}
Definition CgemGeoLayer.h:9
Definition CgemGeoReadoutPlane.h:14
Definition CgemGeomSvc.h:25
CgemGeoReadoutPlane * getReadoutPlane(int iLayer, int iSheet) const
Definition CgemGeomSvc.h:140
Definition ICgemGeomSvc.h:22
The documentation for this class was generated from the following files:
