@@ -458,10 +458,10 @@ class BlipAnaTreeDataStruct
458458 evtTree->Branch (" hit_ismatch" ,hit_ismatch," hit_ismatch[nhits]/I" );
459459 evtTree->Branch (" hit_trkid" ,hit_trkid," hit_trkid[nhits]/I" );
460460 if ( saveTruthInfo ) {
461- evtTree->Branch (" hit_g4trkid" ,hit_g4trkid," hit_g4trkid[nhits]/I" );
462- evtTree->Branch (" hit_g4frac" ,hit_g4frac," hit_g4frac[nhits]/F" );
463- evtTree->Branch (" hit_g4energy" ,hit_g4energy," hit_g4energy[nhits]/F" );
464- evtTree->Branch (" hit_g4charge" ,hit_g4charge," hit_g4charge[nhits]/F" );
461+ evtTree->Branch (" hit_g4trkid" ,hit_g4trkid," hit_g4trkid[nhits]/I" );
462+ evtTree->Branch (" hit_g4frac" ,hit_g4frac," hit_g4frac[nhits]/F" );
463+ evtTree->Branch (" hit_g4energy" ,hit_g4energy," hit_g4energy[nhits]/F" );
464+ evtTree->Branch (" hit_g4charge" ,hit_g4charge," hit_g4charge[nhits]/F" );
465465 }
466466 evtTree->Branch (" hit_clustid" ,hit_clustid," hit_clustid[nhits]/I" );
467467 evtTree->Branch (" hit_blipid" ,hit_blipid," hit_blipid[nhits]/I" );
@@ -747,22 +747,22 @@ class BlipAna : public art::EDAnalyzer
747747
748748 h_part_process = dir_truth.make <TH1D>(" part_process" ," MCParticle->Process()" ,5 ,0 ,5 );
749749 auto xa = h_part_process->GetXaxis ();
750- xa->SetBinLabel (1 ," primary" );
751- xa->SetBinLabel (2 ," compt" );
752- xa->SetBinLabel (3 ," phot" );
753- xa->SetBinLabel (4 ," conv" );
754- xa->SetBinLabel (5 ," other" );
750+ xa->SetBinLabel (1 ," primary" );
751+ xa->SetBinLabel (2 ," compt" );
752+ xa->SetBinLabel (3 ," phot" );
753+ xa->SetBinLabel (4 ," conv" );
754+ xa->SetBinLabel (5 ," other" );
755755
756756 h_nblips_tm = dir_truth.make <TH1D>(" nblips_tm" ," Truth-matched 3D blips per event" ,blipBins,0 ,blipMax);
757757 h_blip_qcomp = dir_truth.make <TH1D>(" blip_qcomp" ," Fraction of true charge (at anode) reconstructed into 3D blips" ,202 ,0 ,1.01 );
758758 h_blip_reszy = dir_truth.make <TH2D>(" blip_res_zy" ," Blip position resolution;Z_{reco} - Z_{true} [cm];Y_{reco} - Y_{true} [cm]" ,150 ,-15 ,15 ,150 ,-15 ,15 );
759- h_blip_reszy ->SetOption (" colz" );
759+ h_blip_reszy ->SetOption (" colz" );
760760 h_blip_resx = dir_truth.make <TH1D>(" blip_res_x" ," Blip position resolution;X_{reco} - X_{true} [cm]" ,200 ,-10 ,10 );
761761 h_blip_resE = dir_truth.make <TH1D>(" blip_res_E" ," Blip energy resolution;(E_{reco} - E_{true})/E_{true} [cm]" ,200 ,-1 .,1 .);
762762 h_blip_E_vs_resE = dir_truth.make <TH2D>(" blip_res_energy" ," Energy resolution of 3D blips;Energy [MeV];#deltaE/E_{true}" ,100 ,0 ,5 ,200 ,-1 .,1 .);
763- h_blip_E_vs_resE ->SetOption (" colz" );
763+ h_blip_E_vs_resE ->SetOption (" colz" );
764764 h_clust_qres_vs_q = dir_truth.make <TH2D>(" qres_vs_q" ," Clusters on collection plane;True charge deposited [ #times 10^{3} e- ];Reco resolution" ,160 ,0 ,80 ,200 ,-1 ,1 );
765- h_clust_qres_vs_q ->SetOption (" colz" );
765+ h_clust_qres_vs_q ->SetOption (" colz" );
766766 h_clust_qres_anode = dir_truth.make <TH1D>(" qres_anode" ," Reco charge vs true charge collected;( reco-true ) / true;Area-normalized entries" ,200 ,-1 .,1 .);
767767 h_clust_qres_dep = dir_truth.make <TH1D>(" qres_dep" ," Reco charge vs true charge deposited;( reco-true ) / true;Area-normalized entries" ,200 ,-1 .,1 .);
768768 h_qratio_vs_time_sim = dir_truth.make <TH2D>(" qratio_vs_time_sim" ," ;Drift time [#mus]; Q_{anode} / Q_{dep}" ,44 ,100 ,2300 , 1000 ,0.50 ,1.50 );
@@ -803,10 +803,10 @@ class BlipAna : public art::EDAnalyzer
803803 h_hitqres[i] = dir_truth.make <TH1D>(Form (" pl%i_hit_q_res" ,i), Form (" Plane %i hits;charge resolution: (reco-true)/true" ,i), 400 ,-1 ,1 );
804804 h_hitqres_scatter[i] = dir_truth.make <TH2D>( Form (" pl%i_hit_qres_scatter" ,i),
805805 Form (" Plane %i;true hit charge [#times 10^{3} e-];Reconstructed hit charge [#times 10^{3} e-]" ,i),qBins,0 ,qMax/1e3 ,qBins,0 ,qMax/1e3 );
806- h_hitqres_scatter[i] ->SetOption (" colz" );
806+ h_hitqres_scatter[i] ->SetOption (" colz" );
807807 h_hitqres_vs_q[i] = dir_truth.make <TH2D>( Form (" pl%i_hit_qres_vs_q" ,i),
808808 Form (" Plane %i;true hit charge [#times 10^{3} e-];hit charge resolution: (reco-true)/true" ,i),qBins,0 ,qMax/1e3 , 400 ,-2 ,2 );
809- h_hitqres_vs_q[i] ->SetOption (" colz" );
809+ h_hitqres_vs_q[i] ->SetOption (" colz" );
810810
811811 h_chargecomp[i] = dir_truth.make <TH1D>(Form (" pl%i_hit_charge_completeness" ,i),Form (" charge completness, plane %i" ,i),101 ,0 ,1.01 );
812812 h_hitpur[i] = dir_truth.make <TH1D>(Form (" pl%i_hit_purity" ,i),Form (" hit purity, plane %i" ,i),101 ,0 ,1.01 );
@@ -1328,7 +1328,6 @@ void BlipAna::analyze(const art::Event& evt)
13281328 }
13291329
13301330 }
1331-
13321331 if ( fDebugMode ) PrintClusterInfo (clust);
13331332
13341333 }// endloop over 2D hit clusters
@@ -1388,7 +1387,7 @@ void BlipAna::analyze(const art::Event& evt)
13881387 nblips_picky++;
13891388 fNum3DBlipsPicky ++;
13901389 h_blip_charge_picky ->Fill (blp.clusters [fCaloPlane ].Charge );
1391- h_blip_zy_picky ->Fill (blp.Position .Z (), blp.Position .Y ());
1390+ h_blip_zy_picky ->Fill (blp.Position .Z (), blp.Position .Y ());
13921391 h_blip_charge_YU_picky->Fill ( 0.001 *blp.clusters [2 ].Charge , 0.001 *blp.clusters [0 ].Charge );
13931392 h_blip_charge_YV_picky->Fill ( 0.001 *blp.clusters [2 ].Charge , 0.001 *blp.clusters [1 ].Charge );
13941393 h_blip_charge_UV_picky->Fill ( 0.001 *blp.clusters [0 ].Charge , 0.001 *blp.clusters [1 ].Charge );
@@ -1491,27 +1490,27 @@ void BlipAna::endJob(){
14911490 // printf(" picky frac : %5.3f\n", fNum3DBlipsPicky/float(fNum3DBlips));
14921491
14931492 if (fIsMC ){
1494- printf (" MC-matched blips per evt : %.3f\n " , fNum3DBlipsTrue /nEvents);
1495- printf (" MC blip purity : %.3f\n " , fNum3DBlipsTrue /float (fNum3DBlips ));
1496- printf (" MC blip purity, 3 planes : %.3f\n " , fNum3DBlipsTrue3P /float (fNum3DBlips3Plane ));
1497- if ( h_blip_qcomp->GetMean () > 0 )
1498- printf (" Charge completeness, total : %.4f +/- %.4f\n " , h_blip_qcomp->GetMean (), h_blip_qcomp->GetStdDev ()/sqrt (fNumEvents ));
1499- // printf(" < 2MeV : %.4f +/- %.4f\n", h_blip_qcomp_2MeV->GetMean(), h_blip_qcomp_2MeV->GetStdDev()/sqrt(fNumEvents));
1500- // printf(" Blip purity : %.4f\n", h_blip_pur->GetMean());
1493+ printf (" MC-matched blips per evt : %.3f\n " , fNum3DBlipsTrue /nEvents);
1494+ printf (" MC blip purity : %.3f\n " , fNum3DBlipsTrue /float (fNum3DBlips ));
1495+ printf (" MC blip purity, 3 planes : %.3f\n " , fNum3DBlipsTrue3P /float (fNum3DBlips3Plane ));
1496+ if ( h_blip_qcomp->GetMean () > 0 )
1497+ printf (" Charge completeness, total : %.4f +/- %.4f\n " , h_blip_qcomp->GetMean (), h_blip_qcomp->GetStdDev ()/sqrt (fNumEvents ));
1498+ // printf(" < 2MeV : %.4f +/- %.4f\n", h_blip_qcomp_2MeV->GetMean(), h_blip_qcomp_2MeV->GetStdDev()/sqrt(fNumEvents));
1499+ // printf(" Blip purity : %.4f\n", h_blip_pur->GetMean());
15011500 }
15021501 printf (" Mean blip charge : %.0f e-\n " , h_blip_charge->GetMean ());
15031502 printf (" \n " );
15041503 for (size_t i=0 ; i<kNplanes ; i++){
1505- printf (" Plane %lu -------------------------\n " ,i);
1506- printf (" * total hits/evt : %.2f\n " ,fNumHits [i]/(float )fNumEvents );
1507- printf (" * untracked hits/evt : %.2f (%.2f plane-matched)\n " ,fNumHitsUntracked [i]/(float )fNumEvents , fNumHitsMatched [i]/(float )fNumEvents );
1508- // printf(" * plane-matched hits/evt : %.2f\n",fNumHitsMatched[i]/(float)fNumEvents);
1509- if (fIsMC ) {
1510- printf (" * true-matched hits/evt : %.2f (%.2f plane-matched)\n " ,fNumHitsTrue [i]/(float )fNumEvents , fNumHitsMatchedTrue [i]/(float )fNumEvents );
1511- if ( h_chargecomp[i]->GetMean () > 0 )
1512- printf (" * charge completeness : %.4f\n " ,h_chargecomp[i]->GetMean ());
1513- printf (" * hit purity : %.4f\n " ,h_hitpur[i]->GetMean ());
1514- }
1504+ printf (" Plane %lu -------------------------\n " ,i);
1505+ printf (" * total hits/evt : %.2f\n " ,fNumHits [i]/(float )fNumEvents );
1506+ printf (" * untracked hits/evt : %.2f (%.2f plane-matched)\n " ,fNumHitsUntracked [i]/(float )fNumEvents , fNumHitsMatched [i]/(float )fNumEvents );
1507+ // printf(" * plane-matched hits/evt : %.2f\n",fNumHitsMatched[i]/(float)fNumEvents);
1508+ if (fIsMC ) {
1509+ printf (" * true-matched hits/evt : %.2f (%.2f plane-matched)\n " ,fNumHitsTrue [i]/(float )fNumEvents , fNumHitsMatchedTrue [i]/(float )fNumEvents );
1510+ if ( h_chargecomp[i]->GetMean () > 0 )
1511+ printf (" * charge completeness : %.4f\n " ,h_chargecomp[i]->GetMean ());
1512+ printf (" * hit purity : %.4f\n " ,h_hitpur[i]->GetMean ());
1513+ }
15151514 }
15161515 printf (" \n ***********************************************\n " );
15171516
@@ -1584,6 +1583,9 @@ void BlipAna::PrintClusterInfo(const blip::HitClust& hc){
15841583 hc.EdepID ,
15851584 hc.isMatched
15861585 );
1586+ printf (" G4 IDs contib \n " );
1587+ for (int g4ID : hc.G4IDs ) printf (" %i " , g4ID);
1588+ printf (" \n " );
15871589}
15881590
15891591void BlipAna::PrintBlipInfo (const blip::Blip& bl){
0 commit comments