@@ -456,10 +456,10 @@ class BlipAnaTreeDataStruct
456456 evtTree->Branch (" hit_ismatch" ,hit_ismatch," hit_ismatch[nhits]/I" );
457457 evtTree->Branch (" hit_trkid" ,hit_trkid," hit_trkid[nhits]/I" );
458458 if ( saveTruthInfo ) {
459- evtTree->Branch (" hit_g4trkid" ,hit_g4trkid," hit_g4trkid[nhits]/I" );
460- evtTree->Branch (" hit_g4frac" ,hit_g4frac," hit_g4frac[nhits]/F" );
461- evtTree->Branch (" hit_g4energy" ,hit_g4energy," hit_g4energy[nhits]/F" );
462- evtTree->Branch (" hit_g4charge" ,hit_g4charge," hit_g4charge[nhits]/F" );
459+ evtTree->Branch (" hit_g4trkid" ,hit_g4trkid," hit_g4trkid[nhits]/I" );
460+ evtTree->Branch (" hit_g4frac" ,hit_g4frac," hit_g4frac[nhits]/F" );
461+ evtTree->Branch (" hit_g4energy" ,hit_g4energy," hit_g4energy[nhits]/F" );
462+ evtTree->Branch (" hit_g4charge" ,hit_g4charge," hit_g4charge[nhits]/F" );
463463 }
464464 evtTree->Branch (" hit_clustid" ,hit_clustid," hit_clustid[nhits]/I" );
465465 evtTree->Branch (" hit_blipid" ,hit_blipid," hit_blipid[nhits]/I" );
@@ -745,22 +745,22 @@ class BlipAna : public art::EDAnalyzer
745745
746746 h_part_process = dir_truth.make <TH1D>(" part_process" ," MCParticle->Process()" ,5 ,0 ,5 );
747747 auto xa = h_part_process->GetXaxis ();
748- xa->SetBinLabel (1 ," primary" );
749- xa->SetBinLabel (2 ," compt" );
750- xa->SetBinLabel (3 ," phot" );
751- xa->SetBinLabel (4 ," conv" );
752- xa->SetBinLabel (5 ," other" );
748+ xa->SetBinLabel (1 ," primary" );
749+ xa->SetBinLabel (2 ," compt" );
750+ xa->SetBinLabel (3 ," phot" );
751+ xa->SetBinLabel (4 ," conv" );
752+ xa->SetBinLabel (5 ," other" );
753753
754754 h_nblips_tm = dir_truth.make <TH1D>(" nblips_tm" ," Truth-matched 3D blips per event" ,blipBins,0 ,blipMax);
755755 h_blip_qcomp = dir_truth.make <TH1D>(" blip_qcomp" ," Fraction of true charge (at anode) reconstructed into 3D blips" ,202 ,0 ,1.01 );
756756 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 );
757- h_blip_reszy ->SetOption (" colz" );
757+ h_blip_reszy ->SetOption (" colz" );
758758 h_blip_resx = dir_truth.make <TH1D>(" blip_res_x" ," Blip position resolution;X_{reco} - X_{true} [cm]" ,200 ,-10 ,10 );
759759 h_blip_resE = dir_truth.make <TH1D>(" blip_res_E" ," Blip energy resolution;(E_{reco} - E_{true})/E_{true} [cm]" ,200 ,-1 .,1 .);
760760 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 .);
761- h_blip_E_vs_resE ->SetOption (" colz" );
761+ h_blip_E_vs_resE ->SetOption (" colz" );
762762 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 );
763- h_clust_qres_vs_q ->SetOption (" colz" );
763+ h_clust_qres_vs_q ->SetOption (" colz" );
764764 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 .);
765765 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 .);
766766 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 );
@@ -801,10 +801,10 @@ class BlipAna : public art::EDAnalyzer
801801 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 );
802802 h_hitqres_scatter[i] = dir_truth.make <TH2D>( Form (" pl%i_hit_qres_scatter" ,i),
803803 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 );
804- h_hitqres_scatter[i] ->SetOption (" colz" );
804+ h_hitqres_scatter[i] ->SetOption (" colz" );
805805 h_hitqres_vs_q[i] = dir_truth.make <TH2D>( Form (" pl%i_hit_qres_vs_q" ,i),
806806 Form (" Plane %i;true hit charge [#times 10^{3} e-];hit charge resolution: (reco-true)/true" ,i),qBins,0 ,qMax/1e3 , 400 ,-2 ,2 );
807- h_hitqres_vs_q[i] ->SetOption (" colz" );
807+ h_hitqres_vs_q[i] ->SetOption (" colz" );
808808
809809 h_chargecomp[i] = dir_truth.make <TH1D>(Form (" pl%i_hit_charge_completeness" ,i),Form (" charge completness, plane %i" ,i),101 ,0 ,1.01 );
810810 h_hitpur[i] = dir_truth.make <TH1D>(Form (" pl%i_hit_purity" ,i),Form (" hit purity, plane %i" ,i),101 ,0 ,1.01 );
@@ -1326,7 +1326,6 @@ void BlipAna::analyze(const art::Event& evt)
13261326 }
13271327
13281328 }
1329-
13301329 if ( fDebugMode ) PrintClusterInfo (clust);
13311330
13321331 }// endloop over 2D hit clusters
@@ -1386,7 +1385,7 @@ void BlipAna::analyze(const art::Event& evt)
13861385 nblips_picky++;
13871386 fNum3DBlipsPicky ++;
13881387 h_blip_charge_picky ->Fill (blp.clusters [fCaloPlane ].Charge );
1389- h_blip_zy_picky ->Fill (blp.Position .Z (), blp.Position .Y ());
1388+ h_blip_zy_picky ->Fill (blp.Position .Z (), blp.Position .Y ());
13901389 h_blip_charge_YU_picky->Fill ( 0.001 *blp.clusters [2 ].Charge , 0.001 *blp.clusters [0 ].Charge );
13911390 h_blip_charge_YV_picky->Fill ( 0.001 *blp.clusters [2 ].Charge , 0.001 *blp.clusters [1 ].Charge );
13921391 h_blip_charge_UV_picky->Fill ( 0.001 *blp.clusters [0 ].Charge , 0.001 *blp.clusters [1 ].Charge );
@@ -1489,27 +1488,27 @@ void BlipAna::endJob(){
14891488 // printf(" picky frac : %5.3f\n", fNum3DBlipsPicky/float(fNum3DBlips));
14901489
14911490 if (fIsMC ){
1492- printf (" MC-matched blips per evt : %.3f\n " , fNum3DBlipsTrue /nEvents);
1493- printf (" MC blip purity : %.3f\n " , fNum3DBlipsTrue /float (fNum3DBlips ));
1494- printf (" MC blip purity, 3 planes : %.3f\n " , fNum3DBlipsTrue3P /float (fNum3DBlips3Plane ));
1495- if ( h_blip_qcomp->GetMean () > 0 )
1496- printf (" Charge completeness, total : %.4f +/- %.4f\n " , h_blip_qcomp->GetMean (), h_blip_qcomp->GetStdDev ()/sqrt (fNumEvents ));
1497- // printf(" < 2MeV : %.4f +/- %.4f\n", h_blip_qcomp_2MeV->GetMean(), h_blip_qcomp_2MeV->GetStdDev()/sqrt(fNumEvents));
1498- // printf(" Blip purity : %.4f\n", h_blip_pur->GetMean());
1491+ printf (" MC-matched blips per evt : %.3f\n " , fNum3DBlipsTrue /nEvents);
1492+ printf (" MC blip purity : %.3f\n " , fNum3DBlipsTrue /float (fNum3DBlips ));
1493+ printf (" MC blip purity, 3 planes : %.3f\n " , fNum3DBlipsTrue3P /float (fNum3DBlips3Plane ));
1494+ if ( h_blip_qcomp->GetMean () > 0 )
1495+ printf (" Charge completeness, total : %.4f +/- %.4f\n " , h_blip_qcomp->GetMean (), h_blip_qcomp->GetStdDev ()/sqrt (fNumEvents ));
1496+ // printf(" < 2MeV : %.4f +/- %.4f\n", h_blip_qcomp_2MeV->GetMean(), h_blip_qcomp_2MeV->GetStdDev()/sqrt(fNumEvents));
1497+ // printf(" Blip purity : %.4f\n", h_blip_pur->GetMean());
14991498 }
15001499 printf (" Mean blip charge : %.0f e-\n " , h_blip_charge->GetMean ());
15011500 printf (" \n " );
15021501 for (size_t i=0 ; i<kNplanes ; i++){
1503- printf (" Plane %lu -------------------------\n " ,i);
1504- printf (" * total hits/evt : %.2f\n " ,fNumHits [i]/(float )fNumEvents );
1505- printf (" * untracked hits/evt : %.2f (%.2f plane-matched)\n " ,fNumHitsUntracked [i]/(float )fNumEvents , fNumHitsMatched [i]/(float )fNumEvents );
1506- // printf(" * plane-matched hits/evt : %.2f\n",fNumHitsMatched[i]/(float)fNumEvents);
1507- if (fIsMC ) {
1508- printf (" * true-matched hits/evt : %.2f (%.2f plane-matched)\n " ,fNumHitsTrue [i]/(float )fNumEvents , fNumHitsMatchedTrue [i]/(float )fNumEvents );
1509- if ( h_chargecomp[i]->GetMean () > 0 )
1510- printf (" * charge completeness : %.4f\n " ,h_chargecomp[i]->GetMean ());
1511- printf (" * hit purity : %.4f\n " ,h_hitpur[i]->GetMean ());
1512- }
1502+ printf (" Plane %lu -------------------------\n " ,i);
1503+ printf (" * total hits/evt : %.2f\n " ,fNumHits [i]/(float )fNumEvents );
1504+ printf (" * untracked hits/evt : %.2f (%.2f plane-matched)\n " ,fNumHitsUntracked [i]/(float )fNumEvents , fNumHitsMatched [i]/(float )fNumEvents );
1505+ // printf(" * plane-matched hits/evt : %.2f\n",fNumHitsMatched[i]/(float)fNumEvents);
1506+ if (fIsMC ) {
1507+ printf (" * true-matched hits/evt : %.2f (%.2f plane-matched)\n " ,fNumHitsTrue [i]/(float )fNumEvents , fNumHitsMatchedTrue [i]/(float )fNumEvents );
1508+ if ( h_chargecomp[i]->GetMean () > 0 )
1509+ printf (" * charge completeness : %.4f\n " ,h_chargecomp[i]->GetMean ());
1510+ printf (" * hit purity : %.4f\n " ,h_hitpur[i]->GetMean ());
1511+ }
15131512 }
15141513 printf (" \n ***********************************************\n " );
15151514
@@ -1582,6 +1581,9 @@ void BlipAna::PrintClusterInfo(const blip::HitClust& hc){
15821581 hc.EdepID ,
15831582 hc.isMatched
15841583 );
1584+ printf (" G4 IDs contib \n " );
1585+ for (int g4ID : hc.G4IDs ) printf (" %i " , g4ID);
1586+ printf (" \n " );
15851587}
15861588
15871589void BlipAna::PrintBlipInfo (const blip::Blip& bl){
0 commit comments