Acceptance

• Created unified script to read reconstructed and thrown information
• Added acceptance calcuation

ecal/scripts/rec_emcal_electrons_reader.C

@@ -2,7 +2,6 @@
 // Read reconstruction ROOT output file
 // Plot variables
 ////////////////////////////////////////
-//int rec_emcal_electrons_reader(double e_start = 1.0, double e_end = 1.0, const char* input_fname = "sim_output/sim_emcal_electrons_output.root", const char *fout = "results/rec_1000k_info_1GeV.txt")
 int rec_emcal_electrons_reader(double e_start = 1.0, double e_end = 1.0, const char* input_fname = "sim_output/sim_emcal_electrons_output.root")
 {
   // Setting figures
@@ -16,11 +15,6 @@ int rec_emcal_electrons_reader(double e_start = 1.0, double e_end = 1.0, const c
   gStyle->SetPadLeftMargin(0.14);
   gStyle->SetPadRightMargin(0.14);
 
-  // Output ASCII file
-  //FILE *fp = fopen(fout,"a");
-  //if(!fp)
-  //      fprintf(stderr,"error: can't open %s\n", fout);
-
   // Input ROOT file
   TFile *f = new TFile(input_fname,"READ");
   TTree *t = (TTree *)f->Get("events");
@@ -29,37 +23,37 @@ int rec_emcal_electrons_reader(double e_start = 1.0, double e_end = 1.0, const c
   t->SetMakeClass(1);
   t->SetBranchStatus("*", 0);
   Int_t mcparticles2_;
-  //t->SetBranchStatus("mcparticles2", 1);
-  //t->SetBranchAddress("mcparticles2", &mcparticles2_);
-  //Int_t CrystalEcalHits2_;
-  //t->SetBranchStatus("CrystalEcalHits2", 1);
-  //t->SetBranchAddress("CrystalEcalHits2", &CrystalEcalHits2_);
+  t->SetBranchStatus("mcparticles2", 1);
+  t->SetBranchAddress("mcparticles2", &mcparticles2_);
+  Int_t CrystalEcalHits2_;
+  t->SetBranchStatus("CrystalEcalHits2", 1);
+  t->SetBranchAddress("CrystalEcalHits2", &CrystalEcalHits2_);
   Int_t RecoEcalHits_;
   t->SetBranchStatus("RecoEcalHits", 1);
   t->SetBranchAddress("RecoEcalHits", &RecoEcalHits_);
   Int_t EcalClusters_;
   t->SetBranchStatus("EcalClusters", 1);
   t->SetBranchAddress("EcalClusters", &EcalClusters_);
-  //const Int_t kMaxmcparticles2 = 100000;
-  //Double_t px[kMaxmcparticles2];
-  //Double_t py[kMaxmcparticles2];
-  //Double_t pz[kMaxmcparticles2];
-  //Double_t mass[kMaxmcparticles2];
-  //t->SetBranchStatus("mcparticles2.psx",1);
-  //t->SetBranchStatus("mcparticles2.psy",1);
-  //t->SetBranchStatus("mcparticles2.psz",1);
-  //t->SetBranchStatus("mcparticles2.mass",1);
-  //t->SetBranchAddress("mcparticles2.psx",px);
-  //t->SetBranchAddress("mcparticles2.psy",py);
-  //t->SetBranchAddress("mcparticles2.psz",pz);
-  //t->SetBranchAddress("mcparticles2.mass",mass);
-  //const Int_t kMaxCrystalEcalHits2 = 100000;
-  //Double_t truth_deposit[kMaxCrystalEcalHits2];
-  //Double_t energyDeposit[kMaxCrystalEcalHits2];
-  //t->SetBranchStatus("CrystalEcalHits2.truth.deposit",1);
-  //t->SetBranchStatus("CrystalEcalHits2.energyDeposit",1);
-  //t->SetBranchAddress("CrystalEcalHits2.truth.deposit",truth_deposit);
-  //t->SetBranchAddress("CrystalEcalHits2.energyDeposit",energyDeposit);
+  const Int_t kMaxmcparticles2 = 100000;
+  Double_t px[kMaxmcparticles2];
+  Double_t py[kMaxmcparticles2];
+  Double_t pz[kMaxmcparticles2];
+  Double_t mass[kMaxmcparticles2];
+  t->SetBranchStatus("mcparticles2.psx",1);
+  t->SetBranchStatus("mcparticles2.psy",1);
+  t->SetBranchStatus("mcparticles2.psz",1);
+  t->SetBranchStatus("mcparticles2.mass",1);
+  t->SetBranchAddress("mcparticles2.psx",px);
+  t->SetBranchAddress("mcparticles2.psy",py);
+  t->SetBranchAddress("mcparticles2.psz",pz);
+  t->SetBranchAddress("mcparticles2.mass",mass);
+  const Int_t kMaxCrystalEcalHits2 = 100000;
+  Double_t truth_deposit[kMaxCrystalEcalHits2];
+  Double_t energyDeposit[kMaxCrystalEcalHits2];
+  t->SetBranchStatus("CrystalEcalHits2.truth.deposit",1);
+  t->SetBranchStatus("CrystalEcalHits2.energyDeposit",1);
+  t->SetBranchAddress("CrystalEcalHits2.truth.deposit",truth_deposit);
+  t->SetBranchAddress("CrystalEcalHits2.energyDeposit",energyDeposit);
   const Int_t kMaxRecoEcalHits = 100000;
   Double_t rec_x_pos[kMaxRecoEcalHits];
   Double_t rec_y_pos[kMaxRecoEcalHits];
@@ -147,8 +141,8 @@ int rec_emcal_electrons_reader(double e_start = 1.0, double e_end = 1.0, const c
 	// Read event by event
 	t->GetEntry(ievent);
 	// Read number of hits/clusters
-	//Int_t nmcparticle      = 2;
-	//Int_t nCrystalEcalHits = CrystalEcalHits2_;
+	Int_t nmcparticle      = 2;
+	Int_t nCrystalEcalHits = CrystalEcalHits2_;
 	Int_t nreconhits       = RecoEcalHits_;
 	Int_t ncluster         = EcalClusters_;
 	// Initialize total energy variables
@@ -157,18 +151,18 @@ int rec_emcal_electrons_reader(double e_start = 1.0, double e_end = 1.0, const c
         total_sim_e       = 0.0;
 	total_cluster_e   = 0.0;
 	// Thrown energy, momentum, and mass
-        //momentum2 = px[nmcparticle]*px[nmcparticle]+py[nmcparticle]*py[nmcparticle]+pz[nmcparticle]*pz[nmcparticle];
-        //momentum  = TMath::Sqrt(momentum2);
-	//mass2     = mass[nmcparticle]*mass[nmcparticle];
-        //total_thr_e = sqrt(momentum2 + mass2)/1.e+3;
-	//h12->Fill(momentum, 1.0);
+        momentum2 = px[nmcparticle]*px[nmcparticle]+py[nmcparticle]*py[nmcparticle]+pz[nmcparticle]*pz[nmcparticle];
+        momentum  = TMath::Sqrt(momentum2);
+	mass2     = mass[nmcparticle]*mass[nmcparticle];
+        total_thr_e = sqrt(momentum2 + mass2)/1.e+3;
+	h12->Fill(momentum, 1.0);
 
         // Loop over simulated hit by simulated hit
-        //for (int isimhit=0; isimhit < nCrystalEcalHits; isimhit++)
-        //{
-        //        total_truth_sim_e += truth_deposit[isimhit];
-        //        total_sim_e       += energyDeposit[isimhit]/1.e+3;
-        //}
+        for (int isimhit=0; isimhit < nCrystalEcalHits; isimhit++)
+        {
+                total_truth_sim_e += truth_deposit[isimhit];
+                total_sim_e       += energyDeposit[isimhit]/1.e+3;
+        }
 
         // Loop over reconstructed hit by reconstructed hit
         for(int ireconhit=0; ireconhit < nreconhits; ireconhit++)
@@ -188,8 +182,6 @@ int rec_emcal_electrons_reader(double e_start = 1.0, double e_end = 1.0, const c
 		cluster_e = cluster_energy[icluster] / 1.e+3;
 		total_cluster_e += cluster_e;
 	}
-	// Writing into ASCII file
-	//fprintf(fp,"%d\t%d\t%2.4f\n",ievent,ncluster,total_cluster_e);
 
 	// Select events with one cluster
         if(ncluster == 1)
@@ -212,8 +204,8 @@ int rec_emcal_electrons_reader(double e_start = 1.0, double e_end = 1.0, const c
 			h11->Fill(eng_res, 1.0);
 		}
 
-		//if(total_cluster_e > 0.9*total_thr_e)
-		//	h13->Fill(momentum, 1.0);
+		if(total_cluster_e > 0.9*total_thr_e)
+			h13->Fill(momentum, 1.0);
 	}
   } 
   // Drawing and Saving figures
@@ -366,6 +358,5 @@ int rec_emcal_electrons_reader(double e_start = 1.0, double e_end = 1.0, const c
   c14->SaveAs("results/hAcceptance.png");
   c14->SaveAs("results/hAcceptance.pdf");
   
-  //fclose(fp);
   return 0;
 }