Commit in hps-java/src/main/java/org/lcsim/hps on MAIN
users/phansson/ecalGainAna.java+106-121.6 -> 1.7
              /ECalGainDriver.java+15-21.8 -> 1.9
recon/ecal/HPSEcalEdepToTriggerConverterDriver.java+5-41.3 -> 1.4
+126-18
3 modified files
ECal calibration with crystal weights
hps-java/src/main/java/org/lcsim/hps/users/phansson
ecalGainAna.java 1.6 -> 1.7
diff -u -r1.6 -r1.7
--- ecalGainAna.java	29 Aug 2012 20:49:23 -0000	1.6
+++ ecalGainAna.java	31 Aug 2012 01:26:33 -0000	1.7
@@ -33,6 +33,7 @@

         option.addOption("f", true, "Input file");
         option.addOption("s", false, "Save to file");
         option.addOption("n", true, "Name added to plots");

+        option.addOption("m", true, "Minimum weight for calibration");

         return option;
     }
 

@@ -83,6 +84,11 @@

             System.exit(1);
         }
 

+        double minCount = 20.0;
+        if (cmd.hasOption("m")) {
+            minCount = Double.valueOf(cmd.getOptionValue("m"));
+        }
+

         ITree sim = null;
         ITree real = null;
         try {

@@ -109,6 +115,11 @@

             }
         }
 

+        IHistogram2D weightPlotSim = (IHistogram2D) sim.find("Weights for correction factor");
+        IHistogram2D sumPlotSim = (IHistogram2D) sim.find("Sum for correction factor");
+        IHistogram2D weightPlotReal = (IHistogram2D) real.find("Weights for correction factor");
+        IHistogram2D sumPlotReal = (IHistogram2D) real.find("Sum for correction factor");
+

 //        IHistogram2D mpePlotSim = (IHistogram2D) sim.find("<E over p>");
 //        IHistogram2D spePlotSim = (IHistogram2D) sim.find("RMS(E over p)");
 //        IHistogram2D mpePlotReal = (IHistogram2D) real.find("<E over p>");

@@ -119,19 +130,90 @@

 
         IHistogramFactory hf = af.createHistogramFactory(null);
 

+        IHistogram2D corrPlotSim = hf.createHistogram2D("Correction factor - simulation", 47, -23.5, 23.5, 11, -5.5, 5.5);
+        IHistogram2D corrPlotReal = hf.createHistogram2D("Correction factor - data", 47, -23.5, 23.5, 11, -5.5, 5.5);
+
+        for (int x = -23; x <= 23; x++) {
+            for (int y = -5; y <= 5; y++) {
+//                if (weightPlotSim.binHeight(x + 23, y + 5) > minCount) {
+                corrPlotSim.fill(x, y, sumPlotSim.binHeight(x + 23, y + 5) / weightPlotSim.binHeight(x + 23, y + 5));
+//                }
+//                if (weightPlotReal.binHeight(x + 23, y + 5) > minCount) {
+                corrPlotReal.fill(x, y, sumPlotReal.binHeight(x + 23, y + 5) / weightPlotReal.binHeight(x + 23, y + 5));
+//                }
+            }
+        }
+
+

         IPlotter plotter = af.createPlotterFactory().create();
         plotter.setParameter("plotterWidth", "600");
         plotter.setParameter("plotterHeight", "400");
         plotter.createRegion();

-        IHistogram2D ecalPlot = hf.createHistogram2D("ECal map", 47, -23.5, 23.5, 11, -5.5, 5.5);
-        plotter.region(0).plot(ecalPlot);

         plotter.region(0).style().statisticsBoxStyle().setVisible(false);
         plotter.region(0).style().setParameter("hist2DStyle", "colorMap");
         plotter.region(0).style().dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
 
 

+        plotter.region(0).plot(weightPlotSim);
+        plotter.region(0).setTitle(weightPlotSim.title() + " - simulation");
+        try {
+            plotter.writeToFile(outName + "_weights_sim.png", "png");
+        } catch (IOException e) {
+            throw new RuntimeException(e);
+        }
+        plotter.region(0).remove(weightPlotSim);
+        plotter.region(0).plot(weightPlotReal);
+        plotter.region(0).setTitle(weightPlotReal.title() + " - data");
+        plotter.region(0).refresh();
+        try {
+            plotter.writeToFile(outName + "_weights_data.png", "png");
+        } catch (IOException e) {
+            throw new RuntimeException(e);
+        }
+        plotter.region(0).remove(weightPlotReal);
+
+//        plotter.region(0).style().zAxisStyle().setParameter("scale", "log");
+        plotter.region(0).plot(corrPlotSim);
+//        plotter.region(0).setZLimits(0, 10.0);
+plotter.region(0).style().zAxisStyle().setParameter("upperLimit", "10.0");
+        plotter.region(0).setTitle(corrPlotSim.title());
+        plotter.region(0).refresh();
+  plotter.region(0).style().zAxisStyle().setParameter("upperLimit", "10.0");
+    try {
+            plotter.writeToFile(outName + "_corr_sim.png", "png");
+        } catch (IOException e) {
+            throw new RuntimeException(e);
+        }
+        plotter.region(0).remove(corrPlotSim);
+        plotter.region(0).plot(corrPlotReal);
+//        plotter.region(0).setZLimits(0, 10.0);
+//for (String opt :plotter.region(0).style().zAxisStyle().availableParameters()) {
+//            System.out.println(opt);
+//        }
+        plotter.region(0).setTitle(corrPlotReal.title());
+        plotter.region(0).refresh();
+        try {
+            plotter.writeToFile(outName + "_corr_data.png", "png");
+        } catch (IOException e) {
+            throw new RuntimeException(e);
+        }
+        plotter.region(0).remove(corrPlotReal);
+//        plotter.region(0).style().zAxisStyle().setParameter("scale", "lin");
+//        plotter.region(0).setZLimits();
+
+
+
+        IHistogram2D ecalPlot = hf.createHistogram2D("ECal map", 47, -23.5, 23.5, 11, -5.5, 5.5);
+        plotter.region(0).plot(ecalPlot);
+
+

         IHistogram1D pePlotsHVSim[][] = new IHistogram1D[4][12];
         IHistogram1D pePlotsHVReal[][] = new IHistogram1D[4][12];

+        double[][] weightsHVSim = new double[4][12];
+        double[][] sumsHVSim = new double[4][12];
+        double[][] weightsHVReal = new double[4][12];
+        double[][] sumsHVReal = new double[4][12];
+

         for (int quad = 1; quad <= 4; ++quad) {
             for (int module = 1; module <= 12; ++module) {
                 pePlotsHVSim[quad - 1][module - 1] = hf.createHistogram1D("E over p quadrant=" + quad + " HV group=" + module, 50, 0, 2);

@@ -142,6 +224,10 @@

             for (int y = -5; y <= 5; y++) {
                 pePlotsHVSim[HPSECalUtils.getQuadrant(x, y) - 1][HPSECalUtils.getHVGroup(x, y) - 1].add(pePlotsSim[x + 23][y + 5][0]);
                 pePlotsHVReal[HPSECalUtils.getQuadrant(x, y) - 1][HPSECalUtils.getHVGroup(x, y) - 1].add(pePlotsReal[x + 23][y + 5][0]);

+                weightsHVSim[HPSECalUtils.getQuadrant(x, y) - 1][HPSECalUtils.getHVGroup(x, y) - 1] += weightPlotSim.binHeight(x + 23, y + 5);
+                sumsHVSim[HPSECalUtils.getQuadrant(x, y) - 1][HPSECalUtils.getHVGroup(x, y) - 1] += sumPlotSim.binHeight(x + 23, y + 5);
+                weightsHVReal[HPSECalUtils.getQuadrant(x, y) - 1][HPSECalUtils.getHVGroup(x, y) - 1] += weightPlotReal.binHeight(x + 23, y + 5);
+                sumsHVReal[HPSECalUtils.getQuadrant(x, y) - 1][HPSECalUtils.getHVGroup(x, y) - 1] += sumPlotReal.binHeight(x + 23, y + 5);

             }
         }
 

@@ -158,7 +244,7 @@

         plotter.region(0).setTitle("Crystal gains before iterating");
         plotter.region(0).refresh();
         try {

-            plotter.writeToFile(outName + "_initial_gains.png", "png");

+            plotter.writeToFile(outName + "_gains_initial.png", "png");

         } catch (IOException e) {
             throw new RuntimeException(e);
         }

@@ -173,7 +259,7 @@

         plotter.region(0).setTitle("Matched cluster counts in simulation");
         plotter.region(0).refresh();
         try {

-            plotter.writeToFile(outName + "_sim_counts.png", "png");

+            plotter.writeToFile(outName + "_counts_sim.png", "png");

         } catch (IOException e) {
             throw new RuntimeException(e);
         }

@@ -188,7 +274,7 @@

         plotter.region(0).setTitle("Matched cluster counts in data");
         plotter.region(0).refresh();
         try {

-            plotter.writeToFile(outName + "_data_counts.png", "png");

+            plotter.writeToFile(outName + "_counts_data.png", "png");

         } catch (IOException e) {
             throw new RuntimeException(e);
         }

@@ -201,8 +287,7 @@

             throw new RuntimeException(e);
         }
 

-        int minCount = 50;
-        

+

         gainStream.format("#x\ty\tgain\n");
         for (int x = -23; x <= 23; x++) {
             for (int y = -5; y <= 5; y++) {

@@ -211,12 +296,21 @@

                 IHistogram1D peReal = pePlotsReal[x + 23][y + 5][0];
                 IHistogram1D peHVSim = pePlotsHVSim[HPSECalUtils.getQuadrant(x, y) - 1][HPSECalUtils.getHVGroup(x, y) - 1];
                 IHistogram1D peHVReal = pePlotsHVReal[HPSECalUtils.getQuadrant(x, y) - 1][HPSECalUtils.getHVGroup(x, y) - 1];

+                double weightHVReal = weightsHVReal[HPSECalUtils.getQuadrant(x, y) - 1][HPSECalUtils.getHVGroup(x, y) - 1];
+                double sumHVReal = sumsHVReal[HPSECalUtils.getQuadrant(x, y) - 1][HPSECalUtils.getHVGroup(x, y) - 1];
+                double weightHVSim = weightsHVSim[HPSECalUtils.getQuadrant(x, y) - 1][HPSECalUtils.getHVGroup(x, y) - 1];
+                double sumHVSim = sumsHVSim[HPSECalUtils.getQuadrant(x, y) - 1][HPSECalUtils.getHVGroup(x, y) - 1];

                 if (gain != null) {

-                    if (peSim.allEntries() >= minCount && peReal.allEntries() >= minCount) {
-                        gain *= (peSim.mean() / peReal.mean());
-                    } else if (peHVSim.allEntries() >= minCount && peHVReal.allEntries() >= minCount) {
-                        gain *= (peHVSim.mean() / peHVReal.mean());

+                    if (weightPlotSim.binHeight(x + 23, y + 5) > minCount && weightPlotReal.binHeight(x + 23, y + 5) > minCount) {
+                        gain *= (corrPlotReal.binHeight(x + 23, y + 5) / corrPlotSim.binHeight(x + 23, y + 5));
+                    } else if (weightHVReal > minCount && weightHVSim > minCount) {
+//                        gain *= (sumHVReal / weightHVReal) / (sumHVSim / weightHVSim);

                     }

+//                    if (peSim.allEntries() >= minCount && peReal.allEntries() >= minCount) {
+//                        gain *= (peSim.mean() / peReal.mean());
+//                    } else if (peHVSim.allEntries() >= minCount && peHVReal.allEntries() >= minCount) {
+////                        gain *= (peHVSim.mean() / peHVReal.mean());
+//                    }

                     gainStream.format("%d\t%d\t%f\n", x, y, gain);
                     ecalPlot.fill(x, y, gain);
 //                    System.out.printf("%d\t%d\t%d\t%f\t%f\n", x, y, pePlotsSim[x + 23][y + 5][0].allEntries(), pePlotsSim[x + 23][y + 5][0].mean(), pePlotsSim[x + 23][y + 5][0].rms());

@@ -228,7 +322,7 @@

         plotter.region(0).setTitle("Crystal gains after iterating");
         plotter.region(0).refresh();
         try {

-            plotter.writeToFile(outName + "_final_gains.png", "png");

+            plotter.writeToFile(outName + "_gains_final.png", "png");

         } catch (IOException e) {
             throw new RuntimeException(e);
         }
hps-java/src/main/java/org/lcsim/hps/users/phansson
ECalGainDriver.java 1.8 -> 1.9
diff -u -r1.8 -r1.9
--- ECalGainDriver.java	29 Aug 2012 20:49:23 -0000	1.8
+++ ECalGainDriver.java	31 Aug 2012 01:26:33 -0000	1.9
@@ -62,6 +62,8 @@

     private IHistogram2D hitmap;
     private IHistogram1D[] h_PE_t = new IHistogram1D[5];
     private IHistogram1D[] h_PE_b = new IHistogram1D[5];

+    private IHistogram2D weightPlot;
+    private IHistogram2D sumPlot;

     int clTop = 0;
     int clBot = 0;
     int trTop = 0;

@@ -136,6 +138,9 @@

         plotter.region(1).style().setParameter("hist2DStyle", "colorMap");
         plotter.region(1).style().dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
 

+        weightPlot = aida.histogram2D("Weights for correction factor", 47, -23.5, 23.5, 11, -5.5, 5.5);
+        sumPlot = aida.histogram2D("Sum for correction factor", 47, -23.5, 23.5, 11, -5.5, 5.5);
+

         for (int iE = 0; iE <= 4; ++iE) {
             for (int irow = -5; irow <= 5; ++irow) {
                 for (int icol = -23; icol <= 23; ++icol) {

@@ -373,8 +378,16 @@

 
         double P = bestTrk.getPX() * 1000;
         double E = bestCl.getEnergy();

-        double Ep = E;
-        double Eoverp = Ep / P;

+        double Eoverp = E / P;
+        double PoverE = P / E;
+
+        for (CalorimeterHit hit : bestCl.getCalorimeterHits()) {
+            double weight = hit.getRawEnergy() / bestCl.getEnergy();
+            int ix = hit.getIdentifierFieldValue("ix");
+            int iy = hit.getIdentifierFieldValue("iy");
+            weightPlot.fill(ix, iy, weight);
+            sumPlot.fill(ix, iy, weight * PoverE);
+        }

 
         if (debug) {
             System.out.println("P " + P + " E " + E);
hps-java/src/main/java/org/lcsim/hps/recon/ecal
HPSEcalEdepToTriggerConverterDriver.java 1.3 -> 1.4
diff -u -r1.3 -r1.4
--- HPSEcalEdepToTriggerConverterDriver.java	27 Aug 2012 21:53:47 -0000	1.3
+++ HPSEcalEdepToTriggerConverterDriver.java	31 Aug 2012 01:26:33 -0000	1.4
@@ -27,6 +27,7 @@

     private double pulseIntegral = tp * Math.E / readoutPeriod;
     //normalization constant from cal gain (MeV/integral bit) to amplitude gain (amplitude bit/GeV)
     private double gainNorm = 1000.0 / pulseIntegral;

+    private double gainScale = 1.0;

 
     public HPSEcalEdepToTriggerConverterDriver() {
     }

@@ -99,7 +100,7 @@

 
     public CalorimeterHit makeTriggerHit(CalorimeterHit hit) {
         long id = hit.getCellID();

-        double amplitude = hit.getRawEnergy() * gainNorm / HPSEcalConditions.physicalToGain(id);

+        double amplitude = hit.getRawEnergy() * gainNorm / HPSEcalConditions.physicalToGain(id) * gainScale;

 
 //        double time = readoutPeriod * (Math.random() - 1);
         double time = 0 - hit.getTime();

@@ -134,12 +135,12 @@

 
     public CalorimeterHit makeReadoutHit(CalorimeterHit hit) {
         long id = hit.getCellID();

-        double amplitude = hit.getRawEnergy() * gainNorm / HPSEcalConditions.physicalToGain(id);

+        double amplitude = hit.getRawEnergy() * gainNorm / HPSEcalConditions.physicalToGain(id) * gainScale;

         if (amplitude < readoutThreshold) {
             return null;
         }

-        double integral = hit.getRawEnergy() * 1000.0;
-        

+        double integral = hit.getRawEnergy() * 1000.0 * gainScale;
+

 //        double thresholdCrossingTime = 0 - hit.getTime();
 //        while (true) {
 //            double currentValue = amplitude * pulseAmplitude(thresholdCrossingTime);
CVSspam 0.2.12

Use REPLY-ALL to reply to list

To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1