diff --git a/ecal/emcal_electrons.sh b/ecal/emcal_electrons.sh index ecaa1f652b2646d1957b73d5e45487e71f729432..9cfc1f3095e2986930f4dc91dcd5b4b1e978a5b9 100644 --- a/ecal/emcal_electrons.sh +++ b/ecal/emcal_electrons.sh @@ -48,6 +48,7 @@ ls -l # run geant4 simulations npsim --runType batch \ -v WARNING \ + --part.minimalKineticEnergy 1000*GeV \ --numberOfEvents ${JUGGLER_N_EVENTS} \ --compactFile ${JUGGLER_DETECTOR}.xml \ --inputFiles ../${JUGGLER_FILE_NAME_TAG}.hepmc \ @@ -64,6 +65,8 @@ root -b -q "ecal/scripts/rec_emcal_electrons_reader.C(${E_start}, ${E_end}, \"${ #root -b -q "ecal/scripts/makeplot.C(${E_start}, ${E_end}, \"${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE}\", \"results/rec_${JUGGLER_FILE_NAME_TAG}.txt\")" #root -b -q "ecal/scripts/makeplot_input.C(\"${JUGGLER_DETECTOR}/${JUGGLER_SIM_FILE}\", \"results/sim_${JUGGLER_FILE_NAME_TAG}.txt\")" +root -b -q "ecal/scripts/crystal_cal_electrons.cxx(\"${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE}\")" + #paste results/sim_${JUGGLER_FILE_NAME_TAG}.txt results/rec_${JUGGLER_FILE_NAME_TAG}.txt > results/eng_${JUGGLER_FILE_NAME_TAG}.txt #root -b -q "ecal/scripts/read_eng.C(\"results/eng_${JUGGLER_FILE_NAME_TAG}.root\", \"results/eng_${JUGGLER_FILE_NAME_TAG}.txt\")" #root -b -q "ecal/scripts/cal_eng_res.C(\"results/eng_${JUGGLER_FILE_NAME_TAG}.root\")" diff --git a/ecal/emcal_pi0s.sh b/ecal/emcal_pi0s.sh index 2387fbf4811194eaedc8156f81f9644006f5deaf..60ef1a1b65e79403945fa5a7c11ca3389714fba2 100644 --- a/ecal/emcal_pi0s.sh +++ b/ecal/emcal_pi0s.sh @@ -41,6 +41,7 @@ echo "JUGGLER_SIM_FILE = ${JUGGLER_SIM_FILE}" # run geant4 simulations npsim --runType batch \ -v WARNING \ + --part.minimalKineticEnergy 1000*GeV \ --numberOfEvents ${JUGGLER_N_EVENTS} \ --compactFile ${JUGGLER_DETECTOR}.xml \ --inputFiles ../${JUGGLER_FILE_NAME_TAG}.hepmc \ diff --git a/ecal/scripts/crystal_cal_electrons.cxx b/ecal/scripts/crystal_cal_electrons.cxx index 959b5150a854ae19ab7fa29a8cb7f40f1b93387c..c1c171eb3c3cb059aee1e988a7ff1a69df9ef2da 100644 --- a/ecal/scripts/crystal_cal_electrons.cxx +++ b/ecal/scripts/crystal_cal_electrons.cxx @@ -50,7 +50,7 @@ auto eta = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) { return result; }; -void crystal_cal_electrons(const char* in_fname = "rec_emcal_uniform_electrons.root") +void crystal_cal_electrons(const char* in_fname = "topside/rec_emcal_uniform_electrons.root") { ROOT::EnableImplicitMT(); ROOT::RDataFrame df("events", in_fname); @@ -62,8 +62,13 @@ void crystal_cal_electrons(const char* in_fname = "rec_emcal_uniform_electrons.r //TH1D* h5 = new TH1D("h5", "Number of Clusters per event", 5, -0.5, 4.5); //TH1D* h6 = new TH1D("h6", "Scattering Angle(#theta) with CUT", 100, 130.0, 180.0); //TH1D* h7 = new TH1D("h7", "Pseudo-rapidity(#eta) with CUT", 100, -5.0, 0.0); - //TH2D* h8 = new TH2D("h8", "Cluster Hit Position", 62, -62.0, 62.0, 62, -62.0, 62.0); - //TH2D* h9 = new TH2D("h9", "All Hit Position", 62, -62.0, 62.0, 62, -62.0, 62.0); + //TH2D *h8 = new TH2D("h8","Cluster Hit Position", 62,-62.0,62.0,62,-62.0,62.0); + //TH2D *h9 = new TH2D("h9","All Hit Position", 62,-62.0,62.0,62,-62.0,62.0); + //TH1D *h10 = new TH1D("hEnergyRes","Energy Resolution", 100,-0.3,0.3); + //TH1D *h11 = new TH1D("hEnergyResCUT","Energy Resolution with CUT", 100,-0.3,0.3); + //TH1D *h12 = new TH1D("h12","Thrown momentum", 61,e_start-0.5,e_end+0.5); + //TH1D *h13 = new TH1D("h13","Thrown momentum for reconstructed particle", 61,e_start-0.5,e_end+0.5); + //TH1D *h14 = new TH1D("h14","Ratio p_{rec}/p_{thr}", 61,e_start-0.5,e_end+0.5); auto d0 = df.Define("isThrown", "mcparticles2.genStatus == 1") .Define("thrownParticles", "mcparticles2[isThrown]") @@ -76,15 +81,34 @@ void crystal_cal_electrons(const char* in_fname = "rec_emcal_uniform_electrons.r for (const auto& i : in) res.push_back(i.energy); return res; - },{"EcalClusters"}); + },{"EcalClusters"}) + .Define("nclusters","EcalClusters.size()") ; + + + auto d1 = d0.Filter("nclusters==1"); auto h_eta_thrown = d0.Histo1D({"h_eta_thrown", " ; #eta ", 100, -5.0, 0.0}, "thrownEta"); auto h_theta_thrown = d0.Histo1D({"h_theta_thrown", "; #theta", 100, 30.0, 180.0},"thrownTheta"); + auto h_nclusters = d0.Histo1D({"h_nclusters", "; N clusters", 6, 0,6},"nclusters"); + auto h_Ecluster = d0.Histo1D({"h_Ecluster", "; cluster E [GeV]",100, 0,30},"Ecluster"); + auto h_Ecluster1 = d1.Histo1D({"h_Ecluster1", "One cluster events; cluster E [GeV]",100, 0,30},"Ecluster"); auto c = new TCanvas(); + h_eta_thrown->DrawCopy(); - c = new TCanvas(); + c->SaveAs("results/crystal_cal_electrons_etaThrown.png"); + h_theta_thrown->DrawCopy(); + c->SaveAs("results/crystal_cal_electrons_thetaThrown.png"); + + h_nclusters->DrawCopy(); + c->SaveAs("results/crystal_cal_electrons_nclusters.png"); + + h_Ecluster->DrawCopy(); + h_Ecluster1->SetLineColor(2); + h_Ecluster1->DrawCopy("same"); + c->SaveAs("results/crystal_cal_electrons_Ecluster.png"); + //std::string outfilename = "rdf_test.root"; //df2.Snapshot("events", outfilename, {"MCParticles_pt", "mcparticles"});