261{
262 MsgStream log(
msgSvc(), name());
263 log << MSG::INFO << "in excute()" << endreq;
264
265
266
267 if(0 == (cgem_totevt % 10000) || m_debug>0) cout << "Event " << cgem_totevt << endl;
268 cgem_totevt++;
269
270 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
271 if (!eventHeader) {
272 log << MSG::FATAL << "Could not find Event Header" << endreq;
273 return( StatusCode::FAILURE);
274 }
275 int iEvt = eventHeader->eventNumber();
276 int iRun = eventHeader->runNumber();
277
278 bool t_ismc = false;
279 if(iRun<0){
280 t_ismc=true;
282 }
283
284 double r[3] = {79.838 ,125.104 ,167.604};
285 double pi=acos(-1.0);
286 vector<double>sim_phi[3];
287 vector<double>sim_Z[3];
288 vector<int>sim_pdg[3];
289 vector<int>sim_track[3];
290 vector<int>sim_hit[3];
291
292 vector<double>phi[3];
293 vector<double>Z[3];
294 vector<int>Truth[3];
295 vector<int>rec_pdg[3];
296
297 vector<int>clusterid[3];
298 vector<int>hit_id[1000];
299 vector<int>track_id[3];
300 vector<int>trkID_S3L6;
301
302
303 map<int,int>isATrk;
304 vector<int>pdg_track;
305 map<int,int>pdg_track_size;
306
307 m_seg_map.clear();
308 N_Seg = 0;
309
310 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
311 if (!recCgemClusterCol){
312 log << MSG::WARNING << "Could not retrieve Cgem cluster list" << endreq;
313 return StatusCode::FAILURE;
314 }
315
316 int count(0);
int ss = 0;
317 int n_1DCluster=0;
318 RecCgemClusterCol::iterator itTrk=recCgemClusterCol->begin();
319 for( cgem_test_nStrip = 0, cgem_XV_ncluster = 0; itTrk != recCgemClusterCol->end(); ++itTrk)
320 {
321 if(ss>=999){N_cluster++;continue;}
322 if(cgem_test_nStrip >= 2999) {N_strip++;continue;}
323 if(cgem_XV_ncluster >= 1999) {N_XV_cluster++;continue;}
324 if((*itTrk)->getflag() == 3) N_flag3++;
325 if(cgem_ntuple){
326
327 if((*itTrk)->getflag() == 0 || (*itTrk)->getflag() == 1)
328 {
329 int cgem_Num_Strip = (*itTrk)->getclusterflage() - (*itTrk)->getclusterflagb() + 1;
330 cgem_nStrip_flagb[cgem_XV_ncluster]=cgem_test_nStrip;
331 cgem_nStrip[cgem_XV_ncluster]=cgem_Num_Strip;
332 for(int ll = 0 ; ll < cgem_Num_Strip ; ll++){
333 cgem_ident_ID[cgem_test_nStrip] =
CgemID::getIntID((*itTrk)->getlayerid(),(*itTrk)->getsheetid(),(*itTrk)->getflag(),ll + (*itTrk)->getclusterflagb());
334 cgem_layer_ID[cgem_test_nStrip] = (*itTrk)->getlayerid();
335 cgem_sheet_ID[cgem_test_nStrip] = (*itTrk)->getsheetid();
336 cgem_strip_ID[cgem_test_nStrip] = ll + (*itTrk)->getclusterflagb();
337 cgem_strip_type[cgem_test_nStrip] = (*itTrk)->getflag();
338 cgem_test_nStrip++;
339 }
340 cgem_XV_clusterid[cgem_XV_ncluster] = (*itTrk)->getclusterid();
341 cgem_XV_flag[cgem_XV_ncluster] = (*itTrk)->getflag();
342 cgem_XV_layerID[cgem_XV_ncluster] = (*itTrk)->getlayerid();
343 cgem_XV_ncluster++;
344 }
345 if((*itTrk)->getflag() == 2)
346 {
347 cgem_evt = iEvt;
348 cgem_run = iRun;
349
350 cgem_clusterid[ss] = (*itTrk)->getclusterid();
351 cgem_layerid[ss] = (*itTrk)->getlayerid();
352 cgem_sheetid[ss] = (*itTrk)->getsheetid();
353 cgem_flag[ss] = (*itTrk)->getflag();
354 cgem_energydeposit[ss] = (*itTrk)->getenergydeposit();
355 cgem_recphi[ss] = (*itTrk)->getrecphi();
356 cgem_recv[ss] = (*itTrk)->getrecv();
357 cgem_recZ[ss] = (*itTrk)->getRecZ();
358 cgem_clusterflag_first[ss] = (*itTrk)->getclusterflagb();
359 cgem_clusterflag_second[ss] = (*itTrk)->getclusterflage();
360
361 }
362 }
363 if((*itTrk)->getflag() == 2)
364 {
365 phi[(*itTrk)->getlayerid()].push_back((*itTrk)->getrecphi());
366 Z[(*itTrk)->getlayerid()].push_back((*itTrk)->getRecZ());
367 ss++;
368 }
369 }
370 if(cgem_ntuple) cgem_ncluster = ss;
371
372
373 cgem_Num21 = 0; cgem_Num23 = 0;
374 for(int ii = 0; ii < Z[1].size(); ii++){
375 for(int jj = 0; jj < Z[0].size(); jj++){
376
377 if(cgem_Num21 == 999)continue;
378 cgem_deltaZ21[cgem_Num21] = Z[0][jj] - Z[1][ii];
379 cgem_deltaPhi21[cgem_Num21] = phi[0][jj] - phi[1][ii];
380 if(cgem_deltaPhi21[cgem_Num21]>-2*
pi && cgem_deltaPhi21[cgem_Num21]<-
pi)cgem_deltaPhi21[cgem_Num21] = cgem_deltaPhi21[cgem_Num21] + 2*
pi;
381 if(cgem_deltaPhi21[cgem_Num21]<2*
pi && cgem_deltaPhi21[cgem_Num21]>
pi) cgem_deltaPhi21[cgem_Num21] = cgem_deltaPhi21[cgem_Num21] - 2*
pi;
382 cgem_Num21++;
383 }
384 for(int ll = 0; ll < Z[2].size(); ll++){
385
386 if(cgem_Num23 == 999)continue;
387 cgem_deltaZ23[cgem_Num23] = Z[2][ll] - Z[1][ii];
388 cgem_deltaPhi23[cgem_Num23] = phi[2][ll] - phi[1][ii];
389 if(cgem_deltaPhi23[cgem_Num23]>-2*
pi && cgem_deltaPhi23[cgem_Num23]<-
pi)cgem_deltaPhi23[cgem_Num23] = cgem_deltaPhi23[cgem_Num23] + 2*
pi;
390 if(cgem_deltaPhi23[cgem_Num23]<2*
pi && cgem_deltaPhi23[cgem_Num23]>
pi) cgem_deltaPhi23[cgem_Num23] = cgem_deltaPhi23[cgem_Num23] - 2*
pi;
391 cgem_Num23++;
392 }
393 }
394
395
396 if(t_ismc){
397 MsgStream log(
msgSvc(), name());
398 log << MSG::INFO << "in TruthExecute()" << endreq;
399
400 SmartDataPtr<Event::CgemMcHitCol> lv_CgemMcHitCol(eventSvc(), "/Event/MC/CgemMcHitCol");
401 if (!lv_CgemMcHitCol)
402 {
403 log << MSG::WARNING << "Could not find event" << endreq;
404 return StatusCode::FAILURE;
405 }
406
407 vector<int> hit_ID;
408 vector<int> Track_ID;
409
410 Event::CgemMcHitCol::iterator iter_truth = lv_CgemMcHitCol->begin();
411 int i=0;
412 for( mc_test_nStrip= 0; iter_truth != lv_CgemMcHitCol->end(); ++iter_truth)
413 {
414 if(i>=1000){N_hit++;continue;}
415
416 m_mc_ID_track[i]=(*iter_truth)->GetTrackID();
417 m_mc_ID_layer[i]=(*iter_truth)->GetLayerID();
418 m_mc_pdg_code[i]=(*iter_truth)->GetPDGCode();
419 m_mc_ID_parent[i]=(*iter_truth)->GetParentID();
420 m_mc_E_deposit[i]=(*iter_truth)->GetTotalEnergyDeposit();
421 m_mc_XYZ_pre_X[i] =(*iter_truth)->GetPositionXOfPrePoint() ;
422 m_mc_XYZ_pre_Y[i] =(*iter_truth)->GetPositionYOfPrePoint() ;
423 m_mc_XYZ_pre_Z[i] =(*iter_truth)->GetPositionZOfPrePoint() ;
424 m_mc_XYZ_post_X[i]=(*iter_truth)->GetPositionXOfPostPoint();
425 m_mc_XYZ_post_Y[i]=(*iter_truth)->GetPositionYOfPostPoint();
426 m_mc_XYZ_post_Z[i]=(*iter_truth)->GetPositionZOfPostPoint();
427 m_mc_P_pre_X[i]=(*iter_truth)->GetMomentumXOfPrePoint() ;
428 m_mc_P_pre_Y[i]=(*iter_truth)->GetMomentumYOfPrePoint() ;
429 m_mc_P_pre_Z[i]=(*iter_truth)->GetMomentumZOfPrePoint() ;
430 m_mc_P_post_X[i]=(*iter_truth)->GetMomentumXOfPostPoint();
431 m_mc_P_post_Y[i]=(*iter_truth)->GetMomentumYOfPostPoint();
432 m_mc_P_post_Z[i]=(*iter_truth)->GetMomentumZOfPostPoint();
433 m_mc_rec_z[i]=(m_mc_XYZ_pre_Z[i]+m_mc_XYZ_post_Z[i])/2;
434 m_mc_rec_phi[i]=(m_mc_XYZ_pre_X[i]+m_mc_XYZ_post_X[i])/2/r[m_mc_ID_layer[i]];
435
436 sim_phi[m_mc_ID_layer[i]].push_back(m_mc_rec_phi[i]);
437 sim_Z[m_mc_ID_layer[i]].push_back(m_mc_rec_z[i]);
438 sim_pdg[m_mc_ID_layer[i]].push_back(m_mc_pdg_code[i]);
439 sim_track[m_mc_ID_layer[i]].push_back(m_mc_ID_track[i]);
440 sim_hit[m_mc_ID_layer[i]].push_back(i);
441
442
443 for(int Pp =0; Pp<m_pid.size(); Pp++)
444 {
445
446 if(m_mc_pdg_code[i] == m_pid[Pp] ){
447 if(m_UseTesterSim){
448 if(m_mc_ID_parent[i] == 0){
449 pdg_track.push_back(m_mc_ID_track[i]);
450 pdg_track_size.insert(pair<int,int> (m_mc_ID_track[i],i));
451 }
452 }
453 else{
454 pdg_track.push_back(m_mc_ID_track[i]);
455 pdg_track_size.insert(pair<int,int> (m_mc_ID_track[i],i));
456 }
457
458 }
459 }
460 m_mc_nStripB[i] = mc_test_nStrip;
461 TArrayI ident = (*iter_truth)->GetIdentifier();
462 for(int j= 0;j<ident.GetSize();j++){
463 if(mc_test_nStrip >= 2999){N_sim_strip++; continue;}
464 m_mc_ident_ID[mc_test_nStrip] = ident.GetAt(j);
465 hit_ID.push_back(i);
466 Track_ID.push_back(m_mc_ID_track[i]);
467 mc_test_nStrip++;
468 }
469 m_mc_nStrip[i] = ident.GetSize();
470 i++;
471
472 }
473 m_mc_nHit=i;
474
475 for(m_mc_nTrk = 0; m_mc_nTrk < t_nTrkMC; m_mc_nTrk++){
476 if(m_mc_nTrk>=100){N_trk++;continue;}
477 m_mc_pidTruth[m_mc_nTrk] = t_pidTruth[m_mc_nTrk] ;
478 m_mc_motheridx[m_mc_nTrk] = t_motheridx[m_mc_nTrk] ;
479 m_mc_costaTruth[m_mc_nTrk] = t_costhetaMC[m_mc_nTrk];
480 m_mc_phiTruth[m_mc_nTrk] = t_phiMC[m_mc_nTrk] ;
481 m_mc_Final_vrTruth[m_mc_nTrk]= t_FvrMC[m_mc_nTrk] ;
482 m_mc_Final_vzTruth[m_mc_nTrk]= t_FvzMC[m_mc_nTrk] ;
483 m_mc_vrTruth[m_mc_nTrk] = t_vrMC[m_mc_nTrk] ;
484 m_mc_vzTruth[m_mc_nTrk] = t_vzMC[m_mc_nTrk] ;
485 m_mc_pxTruth[m_mc_nTrk] = t_pxMC[m_mc_nTrk] ;
486 m_mc_pyTruth[m_mc_nTrk] = t_pyMC[m_mc_nTrk] ;
487 m_mc_pzTruth[m_mc_nTrk] = t_pzMC[m_mc_nTrk] ;
488 m_mc_ptTruth[m_mc_nTrk] = t_ptMC[m_mc_nTrk] ;
489 m_mc_pTruth[m_mc_nTrk] = t_pMC[m_mc_nTrk] ;
490 m_mc_qTruth[m_mc_nTrk] = t_qMC[m_mc_nTrk] ;
491
492 for(int Pp =0; Pp<m_pid.size(); Pp++){
493 if(m_mc_pidTruth[m_mc_nTrk] == m_pid[Pp]&&fabs(m_mc_costaTruth[m_mc_nTrk])>0.93&&m_mc_ptTruth[m_mc_nTrk]<0.05)
494 OutDetecter++;
495 }
496 }
497
498 m_mc_Num21 = 0;m_mc_Num23 = 0;bool sim_IsOut21 = false;bool sim_IsOut23 = false;m_mc_Num1 = sim_Z[0].size();m_mc_Num3 = sim_Z[2].size();
499 for(int Xx = 0; Xx != sim_Z[1].size(); Xx++){
500 for(int Yy = 0; Yy != sim_Z[0].size(); Yy++){
501 if(Yy == 999){sim_IsOut21 = true;continue;}
502 m_mc_deltaZ21[m_mc_Num21] = sim_Z[0][Yy] - sim_Z[1][Xx];
503 m_mc_deltaPhi21[m_mc_Num21] = sim_phi[0][Yy] - sim_phi[1][Xx];
504 if(m_mc_deltaPhi21[m_mc_Num21]>-2*
pi && m_mc_deltaPhi21[m_mc_Num21]<-
pi)m_mc_deltaPhi21[m_mc_Num21] = m_mc_deltaPhi21[m_mc_Num21] + 2*
pi;
505 if(m_mc_deltaPhi21[m_mc_Num21]<2*
pi && m_mc_deltaPhi21[m_mc_Num21]>
pi) m_mc_deltaPhi21[m_mc_Num21] = m_mc_deltaPhi21[m_mc_Num21] - 2*
pi;
506
507 m_mc_hit_ID1[m_mc_Num21] = sim_hit[0][Yy];
508 m_mc_hit_ID2[m_mc_Num21] = sim_hit[1][Xx];
509 if(m_mc_pdg_code[sim_hit[1][Xx]] == m_mc_pdg_code[sim_hit[0][Yy]]) m_mc_hit_PDG21[m_mc_Num21] = m_mc_pdg_code[sim_hit[1][Xx]];
510 else m_mc_hit_PDG21[m_mc_Num21] = -999;
511 if(sim_pdg[1][Xx] == sim_pdg[0][Yy]) m_mc_PDG21[m_mc_Num21] = sim_pdg[1][Xx];
512 else m_mc_PDG21[m_mc_Num21] = -999;
513 if(sim_track[1][Xx] == sim_track[0][Yy]) m_mc_Track21[m_mc_Num21] = sim_track[1][Xx];
514 else m_mc_Track21[m_mc_Num21] = -999;
515 m_mc_Num21++;
516 }
517 for(int Zz = 0; Zz != sim_Z[2].size(); Zz++){
518 if(Zz == 999){sim_IsOut23 = true;continue;}
519 m_mc_deltaZ23[m_mc_Num23] = sim_Z[2][Zz] - sim_Z[1][Xx];
520 m_mc_deltaPhi23[m_mc_Num23] = sim_phi[2][Zz] - sim_phi[1][Xx];
521 if(m_mc_deltaPhi23[m_mc_Num23]>-2*
pi && m_mc_deltaPhi23[m_mc_Num23]<-
pi)m_mc_deltaPhi23[m_mc_Num23] = m_mc_deltaPhi23[m_mc_Num23] + 2*
pi;
522 if(m_mc_deltaPhi23[m_mc_Num23]<2*
pi && m_mc_deltaPhi23[m_mc_Num23]>
pi) m_mc_deltaPhi23[m_mc_Num23] = m_mc_deltaPhi23[m_mc_Num23] - 2*
pi;
523
524 m_mc_hit_ID3[m_mc_Num23] = sim_hit[2][Zz];
525 if(m_mc_pdg_code[sim_hit[1][Xx]] == m_mc_pdg_code[sim_hit[2][Zz]]) m_mc_hit_PDG23[m_mc_Num23] = m_mc_pdg_code[sim_hit[1][Xx]];
526 else m_mc_hit_PDG23[m_mc_Num23] = -999;
527 if(sim_pdg[1][Xx] == sim_pdg[2][Zz]) m_mc_PDG23[m_mc_Num23] = sim_pdg[1][Xx];
528 else m_mc_PDG23[m_mc_Num23] = -999;
529 if(sim_track[1][Xx] == sim_track[2][Zz]) m_mc_Track23[m_mc_Num23] = sim_track[1][Xx];
530 else m_mc_Track23[m_mc_Num23] = -999;
531 m_mc_Num23++;
532 }
533 }
534 if(sim_IsOut21) N_sim_N21++;
535 if(sim_IsOut23) N_sim_N23++;
536
537 for(int m = 0; m < ss; m++){
538 int IsTruth = 0;
539 map<int,int>IsEq;
540 vector<int> tmp_hit_ID[1000];
541 int tmp = 0;
542 int hit_num = 0;
543 int tmp_sim = 0;
544
545
546 for(int rr=0;rr<cgem_XV_ncluster;rr++){
547 if(cgem_XV_clusterid[rr] == cgem_clusterflag_first[m] || cgem_XV_clusterid[rr] == cgem_clusterflag_second[m]){
548
549 for(int id = cgem_nStrip_flagb[rr]; id < cgem_nStrip_flagb[rr] + cgem_nStrip[rr]; id++){
550 for(int sim = 0; sim<mc_test_nStrip; sim++){
551 if(cgem_ident_ID[id]==m_mc_ident_ID[sim]){
552 tmp_hit_ID[hit_ID[sim]].push_back(Track_ID[sim]);
553 }
554 }
555 }
556 tmp += cgem_nStrip[rr];
557 }
558 }
559 for(int Mm = 0; Mm < 1000; Mm++){
560 if(tmp_hit_ID[Mm].size() == 0)continue;
561 IsEq.insert(pair<int,int> (Mm,tmp_hit_ID[Mm][0]));
562 hit_num++;
563 if(tmp_hit_ID[Mm].size() >= tmp) tmp_sim++;
564 }
565
566 bool IsPid ; int pid(0);int tmp_track_ID = 0;
567 for(map<int,int>::iterator Mm = IsEq.begin(); Mm != IsEq.end(); Mm++){
568 IsPid = false;
569 for(int Pp =0; Pp<m_pid.size(); Pp++){
570 if(m_UseTesterSim){
571 if(m_mc_pdg_code[Mm->first] == m_pid[Pp] && m_mc_ID_parent[Mm->first] == 0){
572 IsPid = true;
573 pid = m_mc_pdg_code[Mm->first];
574 }
575 }
576 else{
577 if(m_mc_pdg_code[Mm->first] == m_pid[Pp]){
578 IsPid = true;
579 pid = m_mc_pdg_code[Mm->first];
580 }
581 }
582 }
583 hit_id[m].push_back(Mm->first);
584 tmp_track_ID = Mm->second;
585 }
586 track_id[cgem_layerid[m]].push_back(tmp_track_ID);
587 clusterid[cgem_layerid[m]].push_back(m);
588 if(IsPid) rec_pdg[cgem_layerid[m]].push_back(pid);
589 else rec_pdg[cgem_layerid[m]].push_back(-999);
590 if(IsPid && tmp_sim == 1 )IsTruth =1;
591 if(IsTruth == 1)N_Match_all_hit++;
592 if(IsTruth == 1 && hit_num > 1)N_Match_hit++;
593 Truth[cgem_layerid[m]].push_back(IsTruth);
594 }
595
596 cgem_Num21 = 0;cgem_Num23 = 0;bool IsOut21 = false;bool IsOut23 = false;
597 for(int tt = 0; tt < Truth[1].size(); tt++){
598 for(int hh = 0; hh < Truth[0].size(); hh++){
599 if(cgem_Num21 == 999){IsOut21 = true;continue;}
600 if(rec_pdg[1][tt] != rec_pdg[0][hh] ) cgem_pdg21[cgem_Num21] = -9999;
601 else cgem_pdg21[cgem_Num21] = rec_pdg[1][tt];
602 if(Truth[1][tt] == 0 || Truth[0][hh] == 0)cgem_IsTruth21[cgem_Num21] = 0;
603 if(Truth[1][tt] == 1 && Truth[0][hh] == 1){
604 if(track_id[1][tt] == track_id[0][hh])cgem_IsTruth21[cgem_Num21] = 1;
605 else cgem_IsTruth21[cgem_Num21] = 0;
606 }
607 cgem_Num21++;
608 }
609
610 for(int pp = 0; pp < Truth[2].size(); pp++){
611 if(cgem_Num23 == 999){IsOut23 = true;continue;}
612 if(rec_pdg[1][tt] != rec_pdg[2][pp] ) cgem_pdg23[cgem_Num23] = -9999;
613 else cgem_pdg23[cgem_Num23] = rec_pdg[1][tt];
614 if(Truth[1][tt] == 0 || Truth[2][pp] == 0)cgem_IsTruth23[cgem_Num23] = 0;
615 if(Truth[1][tt] == 1 && Truth[2][pp] == 1){
616 if(track_id[1][tt] == track_id[2][pp])cgem_IsTruth23[cgem_Num23] = 1;
617 else cgem_IsTruth23[cgem_Num23] = 0;
618 }
619 cgem_Num23++;
620 }
621 }
622 if(IsOut21) N_N21++;
623 if(IsOut23) N_N23++;
624
625
626 int tmp_pdg[100] = {0};
627 if(pdg_track_size.size() == 0)N_PiE0++;
628 for(int Nn = 0; Nn < pdg_track.size(); Nn++){
629 int Ii = 0;
630 for(map<int,int>::iterator Mm = pdg_track_size.begin(); Mm != pdg_track_size.end(); Mm++){
631 if(pdg_track[Nn] == Mm->first) tmp_pdg[Ii]+=1;
632 Ii++;
633 }
634 }
635
636 int II = 0;
637 for( map<int,int>::iterator Mm = pdg_track_size.begin(); Mm != pdg_track_size.end(); Mm++){
638 if(tmp_pdg[II] < 3){N_PiL2++;trkID_S3L6.push_back(Mm->first);}
639 if(tmp_pdg[II] > 6){N_PiM7++;trkID_S3L6.push_back(Mm->first);}
640 II++;
641 }
642 }
643
644 if(m_UseCut){
645 N_Event++;
646
647 int cgem_Num1= phi[0].size();int cgem_Num3 = phi[2].size();
648 cut_Num = 0;
649 for(int Ii = 0; Ii < cgem_Num21; Ii++){
650 int tmp_2 = Ii/cgem_Num1;
651 int tmp_1 = Ii%cgem_Num1;
652 for(int Jj = 0; Jj < cgem_Num23; Jj++){
653 int tmp2 = Jj/cgem_Num3;
654 int tmp_3 = Jj%cgem_Num3;
655 if(tmp_2 != tmp2)continue;
656 if(!
InZAr(cgem_deltaZ21[Ii],cgem_deltaZ23[Jj]) || !
InPhiAr(cgem_deltaPhi21[Ii],cgem_deltaPhi23[Jj])){N_OutRange++;
continue;}
657
658
665 int IsTruth = cgem_IsTruth21[Ii]*cgem_IsTruth23[Jj];
667 m_seg_map.insert(pair<int,RecCgemSegment*>(N_Seg,recsegment));
668 N_Seg++;
669
670 cut_deltaZ21[cut_Num] = cgem_deltaZ21[Ii];
671 cut_deltaZ23[cut_Num] = cgem_deltaZ23[Jj];
672 cut_deltaPhi21[cut_Num] = cgem_deltaPhi21[Ii];
673 cut_deltaPhi23[cut_Num] = cgem_deltaPhi23[Jj];
674 cut_IsTruth21[cut_Num] = cgem_IsTruth21[Ii];
675 cut_IsTruth23[cut_Num] = cgem_IsTruth23[Jj];
676 cut_pdg21[cut_Num] = cgem_pdg21[Ii];
677 cut_pdg23[cut_Num] = cgem_pdg23[Jj];
678 cut_Num++;
679
680 if(t_ismc){
681 bool UnFit_Hits = false;
682 for(int Mm = 0; Mm <= 1000; Mm++){
683 if(Mm != clusterid[0][tmp_1] && Mm != clusterid[1][tmp_2] && Mm != clusterid[2][tmp_3]) continue;
684 for(int mm = 0; mm < hit_id[Mm].size(); mm++){
685 for(int Pp =0; Pp<trkID_S3L6.size(); Pp++){
686 if(m_mc_ID_track[hit_id[Mm][mm]] == trkID_S3L6[Pp])UnFit_Hits = true;
687 }
688 }
689 }
690 if(UnFit_Hits) continue;
691 int Un_InTrk =0;
692 if(cgem_IsTruth21[Ii]*cgem_IsTruth23[Jj] != 1){
693 int num = 0;
int N_hit_size = 0;
694 for(int Mm = 0; Mm <= 1000; Mm++){
695 if(Mm != clusterid[0][tmp_1] && Mm != clusterid[1][tmp_2] && Mm != clusterid[2][tmp_3]) continue;
696 N_hit_size += hit_id[Mm].size();
697 for(int mm = 0; mm < hit_id[Mm].size(); mm++){
698 for(int Pp =0; Pp<m_pid.size(); Pp++){
699 if(m_mc_pdg_code[hit_id[Mm][mm]] == m_pid[Pp]){
701 if(sqrt(pow(m_mc_XYZ_pre_X[hit_id[Mm][mm]],2)+pow(m_mc_XYZ_pre_Y[hit_id[Mm][mm]],2))>=sqrt(pow(m_mc_XYZ_post_X[hit_id[Mm][mm]],2)+pow(m_mc_XYZ_post_Y[hit_id[Mm][mm]],2)))Un_InTrk++;
702 }
703 }
704 }
705 }
706 if(
num == 0) N_nonPi++;
707 else{
708 if(
num == N_hit_size){
709 if(Un_InTrk == 0 ||
num == Un_InTrk) N_ErrMat++;
710 else N_CrossMat++;
711 }
712 else{
713 if(
num>3)N_OtherPar4++;
714 if(
num==3)N_OtherPar3++;
715 if(
num<3)N_OtherPar++;
716 }
717 if(
num>3) N_4cluster++;
718 if(
num==3) N_3cluster++;
719 }
720 continue;
721 }
722 N_Match++;
723 int InTrk = 0;
724 for(int Mm = 0; Mm <= 1000; Mm++){
725 if(Mm != clusterid[0][tmp_1] && Mm != clusterid[1][tmp_2] && Mm != clusterid[2][tmp_3]) continue;
726 for(int mm = 0; mm < hit_id[Mm].size(); mm++){
727 if(sqrt(pow(m_mc_XYZ_pre_X[hit_id[Mm][mm]],2)+pow(m_mc_XYZ_pre_Y[hit_id[Mm][mm]],2))>=sqrt(pow(m_mc_XYZ_post_X[hit_id[Mm][mm]],2)+pow(m_mc_XYZ_post_Y[hit_id[Mm][mm]],2)))InTrk++;
728 }
729 }
730 if(InTrk == 3)N_inTrk3++;
731 else if(InTrk < 3&& InTrk > 0)N_InTrk++;
732 else N_outTrk++;
733 }
734 }
735 }
736
738
739 IDataProviderSvc* evtSvc = NULL;
740 Gaudi::svcLocator()->service("EventDataSvc",evtSvc);
741 if (evtSvc) {
742 log << MSG::INFO << "makeTds:event Svc has been found" << endreq;
743 } else {
744 log << MSG::FATAL << "makeTds:Could not find eventSvc" << endreq;
745 return StatusCode::SUCCESS;
746 }
747
748 StatusCode segmentsc;
749 IDataManagerSvc *dataManSvc;
750 dataManSvc= dynamic_cast<IDataManagerSvc*>(evtSvc);
751 DataObject *aSegmentCol;
752 evtSvc->findObject("/Event/Recon/RecCgemSegmentCol",aSegmentCol);
753 if(aSegmentCol != NULL) {
754 dataManSvc->clearSubTree("/Event/Recon/RecCgemSegmentCol");
755 evtSvc->unregisterObject("/Event/Recon/RecCgemSegmentCol");
756 }
757
758 segmentsc = evtSvc->registerObject("/Event/Recon/RecCgemSegmentCol", segmentcol);
759 if( segmentsc.isFailure() ) {
760 log << MSG::FATAL << "Could not register RecCgemSegment" << endreq;
761 return StatusCode::SUCCESS;
762 }
763 log << MSG::INFO << "RecCgemSegmentCol registered successfully!" <<endreq;
764
765
766 for(m_seg_map_it = m_seg_map.begin();m_seg_map_it!=m_seg_map.end();++m_seg_map_it){
768 segmentcol->push_back(recsegment);
769 }
770
771
772 if(m_checkSingleTrkMC>0) cut_ntuple->write();
773 }
774
775 if(m_checkSingleTrkMC>0) cgem_ntuple->write();
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792 return StatusCode::SUCCESS;
793}
ObjectVector< RecCgemSegment > RecCgemSegmentCol
static value_type getIntID(unsigned int f_layer, unsigned int f_sheet, unsigned int f_strip_type, unsigned int f_strip)
bool InZAr(double x, double y)
bool InPhiAr(double x, double y)
void setclusterid_2(const int clusterId_2)
void setmatch(const int match)
void setsegmentid(const int segmentId)
void setclusterid_1(const int clusterId_1)
void seteventid(const int eventId)
void setclusterid_3(const int clusterId_3)