#define MuseSelector_cxx //#include "TH1D.h" //#include "TH2D.h" #include "MuseSelector.h" #include "/Users/noahsteinberg/Physics/anamuse-build/Selectors/MyParameters.h" #include "TList.h" void MuseSelector::Init( TTree *tree) { fReader.SetTree(tree); } void MuseSelector::SlaveBegin(TTree *tree) { FillParams(); //Create histograms TH1::SetDefaultSumw2(); TH2::SetDefaultSumw2(); const int nbin = 100; const double dd = 30; h_vertex_dz_ = new TH1D(("h_vertex_dz_"+id).c_str(), "", nbin, -dd, dd); GetOutputList()->Add(h_vertex_dz_); h_vertex_dx_ = new TH1D(("h_vertex_dx_"+id).c_str(), "", nbin, -dd, dd); GetOutputList()->Add(h_vertex_dx_); h_vertex_dy_ = new TH1D(("h_vertex_dy_"+id).c_str(), "", nbin, -dd, dd); GetOutputList()->Add(h_vertex_dy_); h_doca_ = new TH1D(("h_doca_"+id).c_str(), "", 100, 0, dd); GetOutputList()->Add(h_doca_); h_doca_tgt_ = new TH1D(("h_doca_tgt_"+id).c_str(), "", 100, 0, dd); GetOutputList()->Add(h_doca_tgt_); h_doca_no_tgt_ = new TH1D(("h_doca_no_tgt_"+id).c_str(), "", 100, 0, dd); GetOutputList()->Add(h_doca_no_tgt_); h_dtheta_ = new TH1D(("h_dtheta_"+id).c_str(), "", nbin, -80, 80); GetOutputList()->Add(h_dtheta_); const int nvz = int(2*vz); h_vertex_z_ = new TH1D(("h_vertex_z_"+id).c_str(), "", nvz, -vz, vz); GetOutputList()->Add(h_vertex_z_); h_vertex_z_noBLSC_ = new TH1D(("h_vertex_z_noBLSC_"+id).c_str(), "", nvz, -vz, vz); GetOutputList()->Add(h_vertex_z_noBLSC_); h_vertex_z_tgt_ = new TH1D(("h_vertex_z_tgt_"+id).c_str(), "", nvz, -vz, vz); GetOutputList()->Add(h_vertex_z_tgt_); h_vertex_z_no_tgt_ = new TH1D(("h_vertex_z_no_tgt_"+id).c_str(), "", nvz, -vz, vz); GetOutputList()->Add(h_vertex_z_no_tgt_); h_theta_z_ = new TH2D(("h_theta_z_"+id).c_str(), "", nvz, -vz, vz, 80, 20, 100); GetOutputList()->Add(h_theta_z_); h_theta_z_tgt_ = new TH2D(("h_theta_z_tgt_"+id).c_str(), "", nvz, -vz, vz, 80, 20, 100); GetOutputList()->Add(h_theta_z_tgt_); h_theta_z_no_tgt_ = new TH2D(("h_theta_z_no_tgt_"+id).c_str(), "", nvz, -vz, vz, 80, 20, 100); GetOutputList()->Add(h_theta_z_no_tgt_); h_theta_ = new TH1D(("h_theta_"+id).c_str(), "", 80, 20, 100); GetOutputList()->Add(h_theta_); h_theta_no_rescatter_ = new TH1D(("h_theta_no_rescatter_"+id).c_str(), "", 80, 20, 100); GetOutputList()->Add(h_theta_no_rescatter_); h_theta_tgt_ = new TH1D(("h_theta_tgt_"+id).c_str(), "", 80, 20, 100); GetOutputList()->Add(h_theta_tgt_); h_theta_tgt_no_rescatter_ = new TH1D(("h_theta_tgt_no_rescatter_"+id).c_str(), "", 80, 20, 100); GetOutputList()->Add(h_theta_tgt_no_rescatter_); h_theta_no_tgt_ = new TH1D(("h_theta_no_tgt_"+id).c_str(), "", 80, 20, 100); GetOutputList()->Add(h_theta_no_tgt_); h_theta_no_tgt_no_rescatter_ = new TH1D(("h_theta_no_tgt_no_rescatter_"+id).c_str(), "", 80, 20, 100); GetOutputList()->Add( h_theta_no_tgt_no_rescatter_); h_window_xy_good_ = new TH2D(("h_window_xy_good_"+id).c_str(), "", 100, -250, 250, 200, -400, 400); GetOutputList()->Add(h_window_xy_good_); h_window_xy_bad_ = new TH2D(("h_window_xy_bad_"+id).c_str(), "", 100, -250, 250, 200, -400, 400); GetOutputList()->Add(h_window_xy_bad_); h_vertex_z_intersection_x_good_ = new TH2D(("h_vertex_z_intersection_x_good_"+id).c_str(), "", 100, -100, 100, 100, -100, 100); GetOutputList()->Add(h_vertex_z_intersection_x_good_); h_vertex_z_intersection_x_bad_ = new TH2D(("h_vertex_z_intersection_x_bad_"+id).c_str(), "", 100, -100, 100, 100, -100, 100); GetOutputList()->Add(h_vertex_z_intersection_x_bad_); h_intersection_x_good_ = new TH1D(("h_intersection_x_good_"+id).c_str(), "", 100, -100, 100); GetOutputList()->Add(h_intersection_x_good_); h_intersection_x_bad_ = new TH1D(("h_intersection_x_bad_"+id).c_str(), "", 100, -100, 100); GetOutputList()->Add(h_intersection_x_bad_); h_angle_num_ = new TH1D(Form("h_angle_num_"), "", 80, 20, 100); GetOutputList()->Add(h_angle_num_); h_angle_eff_denom_ = new TH1D(Form("h_angle_eff_denom_"), "", 80, 20, 100); GetOutputList()->Add(h_angle_eff_denom_); h_angle_ratio_ = new TH1D(Form("h_angle_ratio_"), "", 80, 20, 100); GetOutputList()->Add(h_angle_ratio_); h_angle_num_tgt_ = new TH1D(Form("h_angle_num_tgt_m_tgt_"), "", 80, 20, 100); GetOutputList()->Add(h_angle_num_tgt_); h_angle_eff_denom_tgt_ = new TH1D(Form("h_angle_eff_denom_tgt_"), "", 80, 20, 100); GetOutputList()->Add(h_angle_eff_denom_tgt_); h_angle_ratio_tgt_ = new TH1D(Form("h_angle_ratio_tgt_"), "", 80, 20, 100); GetOutputList()->Add(h_angle_ratio_tgt_); h_angle_num_no_tgt_ = new TH1D(Form("h_angle_num_no_tgt_"), "", 80, 20, 100); GetOutputList()->Add(h_angle_num_no_tgt_); h_angle_eff_denom_no_tgt_ = new TH1D(Form("h_angle_eff_denom_no_tgt_"), "", 80, 20, 100); GetOutputList()->Add(h_angle_eff_denom_no_tgt_); h_angle_ratio_no_tgt_ = new TH1D(Form("h_angle_ratio_no_tgt_"), "", 80, 20, 100); GetOutputList()->Add(h_angle_ratio_no_tgt_); h_angle_pur_num_ = new TH1D(Form("h_angle_pur_num_"), "", 80, 20, 100); GetOutputList()->Add(h_angle_pur_num_); h_angle_pur_denom_ = new TH1D(Form("h_angle_pur_denom_"), "", 80, 20, 100); GetOutputList()->Add(h_angle_pur_denom_); h_angle_pur_ratio_ = new TH1D(Form("h_angle_pur_ratio_"), "", 80, 20, 100); GetOutputList()->Add(h_angle_pur_ratio_); h_angle_pur_num_tgt_ = new TH1D(Form("h_angle_pur_num_tgt_m_tgt_"), "", 80, 20, 100); GetOutputList()->Add(h_angle_pur_num_tgt_); h_angle_pur_denom_tgt_ = new TH1D(Form("h_angle_pur_denom_tgt_"), "", 80, 20, 100); GetOutputList()->Add(h_angle_pur_denom_tgt_); h_angle_pur_ratio_tgt_ = new TH1D(Form("h_angle_pur_ratio_tgt_"), "", 80, 20, 100); GetOutputList()->Add(h_angle_pur_ratio_tgt_); h_angle_pur_num_no_tgt_ = new TH1D(Form("h_angle_pur_num_no_tgt_"), "", 80, 20, 100); GetOutputList()->Add(h_angle_pur_num_no_tgt_); h_angle_pur_denom_no_tgt_ = new TH1D(Form("h_angle_pur_denom_no_tgt_"), "", 80, 20, 100); GetOutputList()->Add(h_angle_pur_denom_no_tgt_); h_angle_pur_ratio_no_tgt_ = new TH1D(Form("h_angle_pur_ratio_no_tgt_"), "", 80, 20, 100); GetOutputList()->Add(h_angle_pur_ratio_no_tgt_); //TGT scintillator histograms h_scint_edep_ = new TH1D(("h_scint_edep_"+id).c_str(),"", 50, 0.0, 5.5); GetOutputList()->Add(h_scint_edep_); h_scint_edep_tgt_ = new TH1D(("h_scint_edep_tgt_"+id).c_str(),"", 50, 0.0, 5.5); GetOutputList()->Add(h_scint_edep_tgt_); h_scint_edep_no_tgt_ = new TH1D(("h_scint_edep_no_tgt_"+id).c_str(),"", 50, 0.0, 5.5); GetOutputList()->Add(h_scint_edep_no_tgt_); h_scint_PID_ = new TH1D (("h_scint_PID_"+id).c_str(),"", 25, 0.0, 25.0); GetOutputList()->Add(h_scint_PID_); h_scint_PID_tgt_ = new TH1D (("h_scint_PID_tgt_"+id).c_str(),"", 25, 0.0, 25.0); GetOutputList()->Add(h_scint_PID_tgt_); h_scint_PID_no_tgt_ = new TH1D (("h_scint_PID_no_tgt_"+id).c_str(),"", 25.0, 0.0, 25.0); GetOutputList()->Add(h_scint_PID_no_tgt_); h_scint_CopyID_ = new TH1D(("h_scint_CopyID_"+id).c_str(),"",6, 0, 6); GetOutputList()->Add(h_scint_CopyID_); h_scint_CopyID_tgt_ = new TH1D(("h_scint_CopyID_tgt_"+id).c_str(),"",6, 0, 6); GetOutputList()->Add(h_scint_CopyID_tgt_); h_scint_CopyID_no_tgt_ = new TH1D(("h_scint_CopyID_no_tgt_"+id).c_str(),"",6, 0, 6); GetOutputList()->Add(h_scint_CopyID_no_tgt_); h_edep_num_scints_ = new TH2D(("h_edep_num_scints_"+id).c_str(), "", 60, 0.0, 6.0, 6, 1, 7); GetOutputList()->Add(h_edep_num_scints_); h_edep_vs_theta_ = new TH2D(("h_edep_vs_theta_"+id).c_str(), "", 50, 0.0, 5.5, 32, 20, 48); GetOutputList()->Add(h_edep_vs_theta_); h_zvertex_vs_vtx_momentum_ = new TProfile(("h_zvertex_vs_vtx_momentum_"+id).c_str(),"", 180, -300, 300, ""); GetOutputList()->Add(h_zvertex_vs_vtx_momentum_); h_zvertex_vs_vtx_momentum_tgt_ = new TProfile(("h_zvertex_vs_vtx_momentum_tgt_"+id).c_str(),"", 180, -300, 300, ""); GetOutputList()->Add(h_zvertex_vs_vtx_momentum_tgt_); h_zvertex_vs_vtx_momentum_no_tgt_ = new TProfile(("h_zvertex_vs_vtx_momentum_no_tgt_"+id).c_str(),"", 180, -300, 300, ""); GetOutputList()->Add(h_zvertex_vs_vtx_momentum_no_tgt_); h2_zvertex_vs_vtx_momentum_ = new TH2D(("h2_zvertex_vs_vtx_momentum_"+id).c_str(),"", 180, -300, 300, 120, 0, 120); GetOutputList()->Add(h2_zvertex_vs_vtx_momentum_); h2_zvertex_vs_vtx_momentum_tgt_ = new TH2D(("h2_zvertex_vs_vtx_momentum_tgt_"+id).c_str(),"", 180, -300, 300, 120, 0, 120 ); GetOutputList()->Add(h2_zvertex_vs_vtx_momentum_tgt_); h2_zvertex_vs_vtx_momentum_no_tgt_ = new TH2D(("h2_zvertex_vs_vtx_momentum_no_tgt_"+id).c_str(),"", 180, -300, 300, 120, 0, 120); GetOutputList()->Add(h2_zvertex_vs_vtx_momentum_no_tgt_); h_mom_comp_ = new TH2D(("h_mom_comp_"+id).c_str(), "", 50, 105, 120, 50, 105, 120); GetOutputList()->Add(h_mom_comp_); h_q2_tgt_ = new TH1D(("h_q2_"+id).c_str(), "", 150, .001, .1); GetOutputList()->Add(h_q2_tgt_); h_q2_no_tgt_ = new TH1D(("h_q2_tgt_"+id).c_str(), "", 150, .001, .1); GetOutputList()->Add(h_q2_no_tgt_); h_q2_ = new TH1D(("h_q2_no_tgt_"+id).c_str(), "", 150, .001, .1); GetOutputList()->Add(h_q2_); h_scale_factor_ = new TH1D("scale_factor", "", 80, 20, 100); GetOutputList()->Add(h_scale_factor_); h_overall_eff_num_= new TH1D("overall_eff_num", "", 80, 20, 100); GetOutputList()->Add(h_overall_eff_num_); h_overall_eff_denom_ = new TH1D("overall_eff_denom", "", 80, 20, 100); GetOutputList()->Add(h_overall_eff_denom_); h_subtracted_ = new TH1D("subtracted", "", 80, 20, 100); GetOutputList()->Add(h_subtracted_); h_eff_corrected_ = new TH1D("eff_corected", "", 80, 20, 100); GetOutputList()->Add(h_eff_corrected_); } Bool_t MuseSelector::Process(Long64_t entry) { //Read the next entry fReader.SetEntry(entry); bool scint_hit_over_thresh = false; if (tgt_scint_fired) { double sum = 0; for (int k = 0; k < tgt_scints_hit; k++) { sum += (*tgt_scints_Edep)[k]; } if(sum > scint_threshold) scint_hit_over_thresh = true; } //Reconstructed target events if ( !hit_veto && !frame_hit && !hit_blsc && !MuonDecay && abs(recon_tgt_hit_pos->x()) < 40 && abs(recon_tgt_hit_pos->y()) < 50 && theta < theta_max * TMath::DegToRad() && theta > theta_min * TMath::DegToRad() && abs(delta_bar)<=1 && abs(vertex->z()) < 80 //&& vertex->z() > 200 //&& vertex->z() < 300 && abs(vertex->x()) < 100 && abs(vertex->y()) < 100 && doca < 25 && scint_hit_over_thresh == false ) { int p = 0; if (is_coulomb) {coulomb++; p++;} if (is_eioni) {eioni++; p++;} if (is_ebrem) {ebrem++; p++;} if (p > 1) multiple++; if (p == 0) none++; h_vertex_z_->Fill(vertex->z(), weight); h_doca_->Fill(doca, weight); h_theta_->Fill(theta*TMath::RadToDeg(), weight); h_zvertex_vs_vtx_momentum_->Fill(vertex->z(), p_before_vertex, weight); h2_zvertex_vs_vtx_momentum_->Fill(vertex->z(), p_before_vertex, weight); if (!rescatter_option) h_theta_no_rescatter_->Fill(theta*TMath::RadToDeg(), weight); h_theta_z_->Fill(vertex->z(), theta*TMath::RadToDeg(),weight); h_q2_->Fill(calcQ2(115, theta), weight); // WANT TO TRY AND USE TargetSD FOR EVERYTHING INSTEAD OF STEPPING ACTION if (is_lh2 == true) { h_vertex_z_tgt_->Fill(vertex->z(), weight); //h_doca_tgt_->Fill(doca, weight); h_theta_tgt_->Fill(theta*TMath::RadToDeg(), weight); h_zvertex_vs_vtx_momentum_tgt_->Fill(vertex->z(), p_before_vertex, weight); h2_zvertex_vs_vtx_momentum_tgt_->Fill(vertex->z(), p_before_vertex, weight); if (!rescatter_option) h_theta_tgt_no_rescatter_->Fill(theta*TMath::RadToDeg(), weight); h_q2_tgt_->Fill(calcQ2(115, theta), weight); h_gem_xy_good_->Fill(recon_tgt_hit_pos->x(), recon_tgt_hit_pos->y(), weight); h_theta_z_tgt_->Fill(vertex->z(), theta*TMath::RadToDeg(),weight); } else { h_vertex_z_no_tgt_->Fill(vertex->z(), weight); //h_doca_no_tgt_->Fill(doca, weight); h_theta_no_tgt_->Fill(theta*TMath::RadToDeg(), weight); h_zvertex_vs_vtx_momentum_no_tgt_->Fill(vertex->z(), p_before_vertex, weight); h2_zvertex_vs_vtx_momentum_no_tgt_->Fill(vertex->z(), p_before_vertex, weight); if (!rescatter_option) h_theta_no_tgt_no_rescatter_->Fill(theta*TMath::RadToDeg(), weight); h_q2_no_tgt_->Fill(calcQ2(115, theta), weight); h_gem_xy_bad_->Fill(recon_tgt_hit_pos->x(), recon_tgt_hit_pos->y(), weight); h_theta_z_no_tgt_->Fill(vertex->z(), theta*TMath::RadToDeg(),weight); } if (!hit_blsc) { //if (edep_blsc < 1.0 && wc_p > 0) { h_vertex_z_noBLSC_->Fill(vertex->z(), weight); } //Other optimization plots to separate inscatter events if (!rescatter_option) { h_vertex_z_intersection_x_good_->Fill(ds_intersection_point->x(), vertex->z(), weight); h_intersection_x_good_->Fill(ds_intersection_point->x(), weight); h_window_xy_good_->Fill(window_inter_point->x(), window_inter_point->y(), weight); h_doca_tgt_->Fill(doca, weight); } if (rescatter_option) { h_vertex_z_intersection_x_bad_->Fill(ds_intersection_point->x(), vertex->z(), weight); h_intersection_x_bad_->Fill(ds_intersection_point->x(), weight); h_window_xy_bad_->Fill(window_inter_point->x(), window_inter_point->y(), weight); h_doca_no_tgt_->Fill(doca, weight); } } return kTRUE; } void MuseSelector::Terminate() { h_vertex_z_->Draw("ep"); } void MuseSelector::FillParams() { //Fill Parameters for the run MyParameters *mypar = dynamic_cast(fInput->FindObject("pars")); if (mypar) { mypar->TheNamedInt("vz", vz); mypar->TheNamedDouble("lower_angle", theta_min); mypar->TheNamedDouble("upper_angle", theta_max); mypar->TheNamedDouble("scint_threshold", scint_threshold); } TNamed *identify = dynamic_cast(fInput->FindObject("id")); if (identify) { id = identify->GetTitle(); } TNamed *particle = dynamic_cast(fInput->FindObject("pid")); if (particle) { part = particle->GetTitle(); } TNamed *rescattering = dynamic_cast(fInput->FindObject("rescatter_option")); if (rescattering) { rescatter = rescattering->GetTitle(); } }