add a script for analysis and plot

Analysis/plotAll.C

+{
+TFile *inputFile = new TFile("output.root", "READ"); //JZCaPA output root file
+//Get the central location to align the ZDC1 photons distribution and WC LG efficiency map
+TH1D *xx = new TH1D("xx", "xx", 1000, -50, 50);
+ZDC1tree->Draw("x>>xx");
+double Xmin=0, Xmax=0, Xmean=0;
+for (int j = 0; j < 1000; j++){
+    double x = xx->GetBinContent(j);
+    if (x > 0){
+        if (Xmin == 0){
+            if (x > 0){
+                Xmin = j;
+            }
+        }
+        if (Xmax < j){
+            Xmax = j;
+        }
+    }  
+}
+Xmin = -50+Xmin/10;
+Xmax = -50+Xmax/10;
+Xmean = (Xmin+Xmax)/2;
+
+TH1D *zz = new TH1D("zz", "zz", 2000, 3100, 3300);
+ZDC1tree->Draw("z>>zz");
+double Zmin=0, Zmax=0, Zmean=0;
+for (int j = 0; j < 2000; j++){
+    double z = zz->GetBinContent(j);
+    if (z > 0){
+        if (Zmin == 0){
+            if (z > 0){
+                Zmin = j;
+            }
+        }
+        if (Zmax < j){
+            Zmax = j;
+        }
+    }  
+}
+Zmin = 3100+Zmin/10;
+Zmax = 3100+Zmax/10;
+Zmean = (Zmin+Zmax)/2;
+
+//Set the size of bins and X&Z range, the unit is millimeter
+double Bin = 0.1, halfX = 22.5, halfZ = 73.5; 
+if (Zmin > 3150) {
+    halfZ = 24.5; // for scan1,2,10,11,12 (1x3 WC module)
+}
+//Get ZDC1 photons distribution
+TH2D *hhh = new TH2D("hhh", "ZDC1 Photons Distribution", 2*halfZ/Bin, Zmean-halfZ, Zmean+halfZ, 2*halfX/Bin, Xmean-halfX, Xmean+halfX); // Range size = LG(3x3) size: 147mm x 45mm
+ZDC1tree->Draw("x:z>>hhh");
+gStyle->SetOptStat(0);
+hhh->SetXTitle("Z(mm)");
+hhh->SetYTitle("X(mm)");
+hhh->Draw("COLZ");
+gPad->SaveAs("JZCaPA_Output.pdf");
+
+//Calculate photons number in each single WC module
+int Zsize = 3;
+if (Zmin > 3150) {
+    Zsize = 1; // for scan1,2,10,11,12 (1x3 WC module)
+}
+TH2D *k1 = new TH2D("k1", "Single WC Module", Zsize, Zmin, Zmax, 3, Xmin, Xmax);
+ZDC1tree->Draw("x:z>>k1");
+gStyle->SetOptStat(0);
+k1->SetXTitle("Z(mm)");
+k1->SetYTitle("X(mm)");
+k1->Draw("COLZ");
+gPad->SaveAs("x3.pdf");
+
+
+double PhotonsInput = 50000000; //Pure light guide simulation input photons number
+inputFile = new TFile("temp.root", "READ"); //Input pure LG simulation root file 
+TH2D *mmm = new TH2D("mmm","LG Efficiency Map", 2*halfZ/Bin, -halfZ, halfZ, 2*halfX/Bin, -halfX, halfX);
+lightGuide->Draw("Z:X>>mmm");
+//Inject in the simulation uniform light but just over the rods surface (Just get photons on the rods surface)
+double RodsPhotonsOut = 0, AreaCounts = 0;
+TH2D *k2 = new TH2D("k2", "k2", 2*halfZ/Bin, Zmean-halfZ, Zmean+halfZ, 2*halfX/Bin, Xmean-halfX, Xmean+halfX); 
+for(int i = 1; i <= mmm->GetNbinsX(); i++){ 
+    for(int j = 1; j <= mmm->GetNbinsY(); j++){
+        if (mmm->GetBinContent(i, j) * hhh->GetBinContent(i, j) != 0){
+            k2->SetBinContent(i, j, mmm->GetBinContent(i, j));
+            RodsPhotonsOut += k2->GetBinContent(i, j);
+            AreaCounts++;
+        }
+    }
+}
+gStyle->SetOptStat(0);
+k2->SetXTitle("Z(mm)");
+k2->SetYTitle("X(mm)");
+k2->Draw("COLZ");
+gPad->SaveAs("Rods_Surface.pdf");
+std::cout << "==========================================================" << std::endl;
+std::cout << "===== RodsPhotonsOut: " << RodsPhotonsOut << " =====" << std::endl;
+std::cout << "===== Rods_Abs_Eff: " << RodsPhotonsOut/(AreaCounts*Bin*Bin/147/45*PhotonsInput) << " =====" << std::endl; //Input photons ratio of rods surface area to total area
+mmm->Scale(1.0/mmm->GetMaximum());
+mmm->Draw("COLZ");
+mmm->SetXTitle("Z(mm)");
+mmm->SetYTitle("X(mm)");
+mmm->SetAxisRange(0., 1.0, "Z");
+gPad->SaveAs("LG_Efficiency_Map.pdf");
+
+//Plot the final JZCaPA output with applied efficiency map 
+TH2D *fff = new TH2D("fff", "Relative Efficiency", 2*halfZ/Bin, Zmean-halfZ, Zmean+halfZ, 2*halfX/Bin, Xmean-halfX, Xmean+halfX); 
+for(int i = 1; i <= mmm->GetNbinsX(); i++){ 
+    for(int j = 1; j <= mmm->GetNbinsY(); j++){
+        fff->SetBinContent(i, j, hhh->GetBinContent(i, j) * mmm->GetBinContent(i, j));
+    }
+}
+
+//Calculate absolute efficiency of full WC modules (3x3)
+double JZCaPA_Photon_Out = 0;
+for(int i = 1; i <= fff->GetNbinsX(); i++){ 
+  for(int j = 1; j <= fff->GetNbinsY(); j++){
+    JZCaPA_Photon_Out += fff->GetBinContent(i, j);
+  }
+}
+std::cout << "===== Absolute Efficiency of Full Modules : " << JZCaPA_Photon_Out/xx->GetEntries() << " =====" << std::endl;
+
+//Get photons number of each single WC after applied efficiency map
+double wc[3][3];
+int k=1;
+for(int i = 1; i <= 1470; i++){ 
+    int t=1;
+    if (i >= k*490){
+        k++;
+    }
+    for(int j = 1; j <= 450; j++){
+        if (j >= t*150){
+            t++;
+        }
+        wc[k-1][t-1] += fff->GetBinContent(i, j);
+    }
+}
+
+std::cout << "===== Absolute Efficiency of Single WC =====" << std::endl;
+std::cout << "X-axis" << std::endl;
+if (Zmin > 3150){ //for scan1,2,10,11,12 (1x3 WC module)
+    std::cout << "|" << wc[0][2]/k1->GetBinContent(1, 3) << std::endl;
+    std::cout << "|" << wc[0][1]/k1->GetBinContent(1, 2) << std::endl;
+    std::cout << "|" << wc[0][0]/k1->GetBinContent(1, 1) << std::endl;
+}else{ //for scan3,4,6,7,8,9,13 (3x3 WC module)
+    std::cout << "|" << wc[0][2]/k1->GetBinContent(1, 3) << ", " << wc[1][2]/k1->GetBinContent(2, 3) << ", " << wc[2][2]/k1->GetBinContent(3, 3) << std::endl;
+    std::cout << "|" << wc[0][1]/k1->GetBinContent(1, 2) << ", " << wc[1][1]/k1->GetBinContent(2, 2) << ", " << wc[2][1]/k1->GetBinContent(3, 2) << std::endl;
+    std::cout << "|" << wc[0][0]/k1->GetBinContent(1, 1) << ", " << wc[1][0]/k1->GetBinContent(2, 1) << ", " << wc[2][0]/k1->GetBinContent(3, 1) << std::endl;
+}
+std::cout << "--------------------------->Z-axis" << std::endl;
+std::cout << "==========================================================" << std::endl;
+
+fff->SetXTitle("Z(mm)");
+fff->SetYTitle("X(mm)");
+fff->Scale(1.0/fff->GetMaximum());
+fff->SetAxisRange(0., 1.0, "Z");
+fff->Draw("COLZ");
+gPad->SaveAs("Relative Efficiency.pdf");
+
+}