Scripts

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"});