lcsim/src/org/lcsim/contrib/Cassell/recon/Cheat
diff -N QCalibrationFromData.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ QCalibrationFromData.java 8 Sep 2008 23:11:30 -0000 1.1
@@ -0,0 +1,544 @@
+package org.lcsim.contrib.Cassell.recon.Cheat;
+import org.lcsim.recon.cluster.analysis.*;
+import java.util.*;
+import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.Subdetector;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.util.*;
+import org.lcsim.event.*;
+import org.lcsim.event.base.*;
+import org.lcsim.recon.cluster.util.*;
+import org.lcsim.geometry.IDDecoder;
+import org.lcsim.util.aida.AIDA;
+import hep.physics.vec.*;
+import Jama.Matrix;
+
+/**
+ * Calculate neutral sampling fractions for each calorimeter using
+ * physics data minimizing event dE/rootE squared, then a response
+ * correction
+ * Ron Cassell
+ */
+public class QCalibrationFromData extends Driver
+{
+ private AIDA aida = AIDA.defaultInstance();
+ IDDecoder _decoder;
+ int ievt;
+ int mx = 2000000;
+ double[][] Aph;
+ double[][] Bph;
+ double[][] Anh;
+ double[][] Bnh;
+ double ctcut = .9;
+ double Ecutnh = 2.;
+ double Ecutg = 1.;
+ double phEtot;
+ double nhEtot;
+ int hitcut = 10;
+ int nsumsph = 4;
+ int nsumsnh = 6;
+ double alphanh = -.23;
+ double alphaph = .036;
+ double phres = .17;
+ double nhres = .7;
+ int nphcbins = 20;
+ int nnhcbins = 15;
+ int nfrontecallayers = 21;
+ List<MyDoubleArray> phl;
+ List<MyDoubleArray> nhl;
+ String CheatReconRname = "ReconPerfectReconParticles";
+ String PPRPflowRname = "PPRReconParticles";
+//
+// Only allow 1 detector
+//
+ boolean dchanged;
+
+ public QCalibrationFromData()
+ {
+ Aph = new double[nsumsph][nsumsph];
+ Bph = new double[nsumsph][1];
+ Anh = new double[nsumsnh][nsumsnh];
+ Bnh = new double[nsumsnh][1];
+ phl = new ArrayList<MyDoubleArray>();
+ nhl = new ArrayList<MyDoubleArray>();
+ phEtot = 0.;
+ nhEtot = 0.;
+ dchanged = false;
+//
+// Do the cheating reconstruction (includes Digisim)
+//
+ CheatReconDriver crd = new CheatReconDriver();
+ crd.setCheatReconstructedParticleOutputName(CheatReconRname);
+ add(crd);
+ PPRDriver pprd = new PPRDriver(CheatReconRname,PPRPflowRname);
+ TrivialClusterEnergyCalculator t = new TrivialClusterEnergyCalculator();
+ pprd.setPhotonClusterEnergyCalculator(t);
+ pprd.setNeutralHadronClusterEnergyCalculator(t);
+ add(pprd);
+ ievt = 0;
+ }
+ protected void process(EventHeader event)
+ {
+ System.out.println(" Processing event "+ievt);
+ super.process(event);
+ ievt ++;
+ List<ReconstructedParticle> rlist = event.get(ReconstructedParticle.class,PPRPflowRname);
+//
+// Pick out the photons and neutral hadrons
+//
+ for(ReconstructedParticle p:rlist)
+ {
+ double E = ( (MCReconstructedParticle) p.getParticles().get(0)).getMCParticle().getEnergy();
+ Hep3Vector mom = ( (MCReconstructedParticle) p.getParticles().get(0)).getMCParticle().getMomentum();
+ double Mass = ( (MCReconstructedParticle) p.getParticles().get(0)).getMCParticle().getMass();
+ double Charge = ( (MCReconstructedParticle) p.getParticles().get(0)).getMCParticle().getCharge();
+ Hep3Vector Origin = ( (MCReconstructedParticle) p.getParticles().get(0)).getMCParticle().getOrigin();
+ if( (p.getMass() == 0.)&&(p.getCharge() == 0.)&&(Mass == 0) )
+ {
+ if( Math.abs(mom.z())/mom.magnitude() < ctcut)
+ {
+//
+// Have a photon to use
+//
+ double dE = p.getEnergy() - E;
+ double[] gsums = new double[nsumsph];
+ for(Cluster c:p.getClusters())
+ {
+ for(CalorimeterHit hit:c.getCalorimeterHits())
+ {
+ IDDecoder idc = hit.getIDDecoder();
+ idc.setID(hit.getCellID());
+ int detector_index = idc.getValue("system");
+ double[] pos = hit.getPosition();
+ double R = Math.sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]);
+ double ctheta = Math.abs(pos[2])/R;
+ double st = Math.sqrt(1.-ctheta*ctheta);
+ double cb = 1. + alphaph*(1./st - 1.);
+ double ce = 1. + alphaph*(1./ctheta - 1.);
+ if(detector_index == 2)
+ {
+ int layer = idc.getValue("layer");
+ if(layer < nfrontecallayers)
+ {
+ gsums[0] += hit.getRawEnergy()/cb;
+ }
+ else
+ {
+ gsums[1] += hit.getRawEnergy()/cb;
+ }
+ }
+ else if(detector_index == 6)
+ {
+ int layer = idc.getValue("layer");
+ if(layer < nfrontecallayers)
+ {
+ gsums[2] += hit.getRawEnergy()/ce;
+ }
+ else
+ {
+ gsums[3] += hit.getRawEnergy()/ce;
+ }
+ }
+ }
+ }
+ phl.add(new MyDoubleArray(gsums,E));
+ if(E > Ecutg)
+ {
+ phEtot += E;
+ for(int row=0;row<nsumsph;row++)
+ {
+ Bph[row][0] += gsums[row];
+ for(int col=0;col<nsumsph;col++)
+ {
+ Aph[row][col] += gsums[row]*gsums[col]/E;
+ }
+ }
+ }
+ }
+ }
+ else if( (Mass > 0.)&&(Mass < 1.)&&(p.getCharge() == 0.)&&(Charge == 0)&&(Origin.magnitude() < 10.) )
+ {
+ if( Math.abs(mom.z())/mom.magnitude() < ctcut)
+ {
+//
+// Have a neutral hadron to use
+//
+ double[] gsums = new double[nsumsnh];
+ for(Cluster c:p.getClusters())
+ {
+ for(CalorimeterHit hit:c.getCalorimeterHits())
+ {
+ IDDecoder idc = hit.getIDDecoder();
+ idc.setID(hit.getCellID());
+ int detector_index = idc.getValue("system");
+ if(detector_index == 2)
+ {
+ int layer = idc.getValue("layer");
+ if(layer < nfrontecallayers)
+ {
+ gsums[0] += hit.getRawEnergy();
+ }
+ else
+ {
+ gsums[1] += hit.getRawEnergy();
+ }
+ }
+ else if(detector_index == 6)
+ {
+ int layer = idc.getValue("layer");
+ if(layer < nfrontecallayers)
+ {
+ gsums[2] += hit.getRawEnergy();
+ }
+ else
+ {
+ gsums[3] += hit.getRawEnergy();
+ }
+ }
+ else if(detector_index == 3)
+ {
+ double[] pos = hit.getPosition();
+ double R = Math.sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]);
+ double ctheta = Math.abs(pos[2])/R;
+ double st = Math.sqrt(1.-ctheta*ctheta);
+ gsums[4] += 1./(1. + alphanh*(1./st - 1.));
+ }
+ else if(detector_index == 7)
+ {
+ double[] pos = hit.getPosition();
+ double R = Math.sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]);
+ double st = Math.abs(pos[2])/R;
+ gsums[5] += 1./(1. + alphanh*(1./st - 1.));
+ }
+ }
+ }
+ nhl.add(new MyDoubleArray(gsums,E));
+ if(E > Ecutnh)
+ {
+ nhEtot += E;
+ for(int row=0;row<nsumsnh;row++)
+ {
+ Bnh[row][0] += gsums[row];
+ for(int col=0;col<nsumsnh;col++)
+ {
+ Anh[row][col] += gsums[row]*gsums[col]/E;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ protected void suspend()
+ {
+ Matrix mA = new Matrix(Aph);
+ Matrix mB = new Matrix(Bph);
+ Matrix x = mA.solve(mB);
+ double[] phsf = new double[nsumsph];
+ double phEmtot = 0.;
+ for(int i=0;i<nsumsph;i++)
+ {
+ phEmtot += x.get(i,0)*Bph[i][0];
+ }
+ double rat = phEmtot/phEtot;
+ for(int i=0;i<nsumsph;i++)
+ {
+ phsf[i] = rat/x.get(i,0);
+ }
+ System.out.println("*** Photons ***;");
+ System.out.println("nFrontLayersEcal: "+nfrontecallayers);
+ System.out.println("alpha: "+alphaph);
+ String sfout = "SFs: ";
+ for(int i=0;i<nsumsph-1;i++)
+ {
+ sfout += phsf[i]+",";
+ }
+ sfout += phsf[nsumsph-1];
+ System.out.println(sfout);
+ double Em = 0.;
+ for(int i=0;i<nsumsph;i++)
+ {
+ Em += phl.get(0).getVal(i)/phsf[i];
+ }
+ List<Double> phEm = new ArrayList<Double>();
+ List<Double> phE = new ArrayList<Double>();
+ Double xx = new Double(Em);
+ phEm.add(0,xx);
+ phE.add(0,new Double(phl.get(0).getE()));
+ for(int i=1;i<phl.size();i++)
+ {
+ Em = 0.;
+ for(int j=0;j<nsumsph;j++)
+ {
+ Em += phl.get(i).getVal(j)/phsf[j];
+ }
+ int k = 0;
+ for(int j=0;j<phEm.size();j++)
+ {
+ if(Em > phEm.get(j).doubleValue())break;
+ k++;
+ }
+ phEm.add(k,new Double(Em));
+ phE.add(k,new Double(phl.get(i).getE()));
+ }
+ double totEm = 0.;
+ for(int i=0;i<phEm.size();i++)
+ {
+ totEm += phEm.get(i).doubleValue();
+ }
+ double[] sEm = new double[nphcbins];
+ double[] sE = new double[nphcbins];
+ int[] N = new int[nphcbins];
+ double Emt = 0.;
+ int bin = 0;
+ double[] pbmax = {0.2,0.4,0.7,1.5,3.,8.,15.,40.,65.,1000.};
+ for(int i=0;i<phEm.size();i++)
+ {
+ double Egen = phE.get(i).doubleValue();
+ double Emeas = phEm.get(i).doubleValue();
+ double dE = Emeas - Egen;
+ int pbin = 0;
+ for(int k=0;k<pbmax.length;k++)
+ {
+ if(Emeas < pbmax[k])break;
+ pbin++;
+ }
+ aida.cloud1D("Photon/E",mx).fill(Egen);
+ aida.cloud1D("Photon/dE",mx).fill(dE);
+ aida.cloud1D("Photon/dE over E",mx).fill(dE/Egen);
+ aida.cloud1D("Photon/dE over rootE",mx).fill(dE/Math.sqrt(Egen));
+ aida.cloud1D("Photon/wted E",mx).fill(Egen,Egen);
+ aida.cloud1D("Photon/wted dE",mx).fill(dE,Egen);
+ aida.cloud1D("Photon/wted dE over E",mx).fill(dE/Egen,Egen);
+ aida.cloud1D("Photon/wted dE over rootE",mx).fill(dE/Math.sqrt(Egen),Egen);
+ aida.cloud1D("Photon/bin "+pbin+" dE",mx).fill(dE);
+ aida.cloud1D("Photon/bin "+pbin+" E",mx).fill(Egen);
+ aida.cloud1D("Photon/bin "+pbin+" dE over E",mx).fill(dE/Egen);
+ aida.cloud1D("Photon/bin "+pbin+" dE over rootE",mx).fill(dE/Math.sqrt(Egen));
+ aida.cloud1D("Photon/bin "+pbin+" wted E",mx).fill(Egen,Egen);
+ aida.cloud1D("Photon/bin "+pbin+" wted dE",mx).fill(dE,Egen);
+ aida.cloud1D("Photon/bin "+pbin+" wted dE over E",mx).fill(dE/Egen,Egen);
+ aida.cloud1D("Photon/bin "+pbin+" wted dE over rootE",mx).fill(dE/Math.sqrt(Egen),Egen);
+ sEm[bin] += phEm.get(i).doubleValue();
+ sE[bin] += phE.get(i).doubleValue();
+ N[bin]++;
+ Emt += phEm.get(i).doubleValue();
+ if(Emt > totEm/nphcbins*(bin+1))bin++;
+ }
+ double[] mEm = new double[nphcbins];
+ double[] mE = new double[nphcbins];
+ for(int i=0;i<nphcbins;i++)
+ {
+ int j = nphcbins - 1 - i;
+ mEm[j] = sEm[i]/N[i];
+ mE[j] = sE[i]/N[i];
+ }
+ String Emstr = "Emeas: ";
+ String Ecstr = "Ecorr: ";
+ for(int i=0;i<nphcbins-1;i++)
+ {
+ Emstr += mEm[i]+",";
+ Ecstr += mE[i]+",";
+ }
+ Emstr += mEm[nphcbins-1];
+ Ecstr += mE[nphcbins-1];
+ System.out.println(Emstr);
+ System.out.println(Ecstr);
+ mA = new Matrix(Anh);
+ mB = new Matrix(Bnh);
+ x = mA.solve(mB);
+ double[] nhsf = new double[nsumsnh];
+ double nhEmtot = 0.;
+ for(int i=0;i<nsumsnh;i++)
+ {
+ nhEmtot += x.get(i,0)*Bnh[i][0];
+ }
+ rat = nhEmtot/nhEtot;
+ for(int i=0;i<nsumsnh;i++)
+ {
+ nhsf[i] = rat/x.get(i,0);
+ }
+ System.out.println("*** Neutral hadrons ***;");
+ System.out.println("nFrontLayersEcal: "+nfrontecallayers);
+ System.out.println("alpha: "+alphanh);
+ sfout = "SFs: ";
+ for(int i=0;i<nsumsnh-1;i++)
+ {
+ sfout += nhsf[i]+",";
+ }
+ sfout += nhsf[nsumsnh-1];
+ System.out.println(sfout);
+ Em = 0.;
+ for(int i=0;i<nsumsnh;i++)
+ {
+ Em += nhl.get(0).getVal(i)/nhsf[i];
+ }
+ List<Double> nhEm = new ArrayList<Double>();
+ List<Double> nhE = new ArrayList<Double>();
+ xx = new Double(Em);
+ nhEm.add(0,xx);
+ nhE.add(0,new Double(nhl.get(0).getE()));
+ for(int i=1;i<nhl.size();i++)
+ {
+ Em = 0.;
+ for(int j=0;j<nsumsnh;j++)
+ {
+ Em += nhl.get(i).getVal(j)/nhsf[j];
+ }
+ int k = 0;
+ for(int j=0;j<nhEm.size();j++)
+ {
+ if(Em > nhEm.get(j).doubleValue())break;
+ k++;
+ }
+ nhEm.add(k,new Double(Em));
+ nhE.add(k,new Double(nhl.get(i).getE()));
+ }
+ totEm = 0.;
+ for(int i=0;i<nhEm.size();i++)
+ {
+ totEm += nhEm.get(i).doubleValue();
+ }
+ sEm = new double[nnhcbins];
+ sE = new double[nnhcbins];
+ N = new int[nnhcbins];
+ Emt = 0.;
+ bin = 0;
+ for(int i=0;i<nhEm.size();i++)
+ {
+ sEm[bin] += nhEm.get(i).doubleValue();
+ sE[bin] += nhE.get(i).doubleValue();
+ N[bin]++;
+ Emt += nhEm.get(i).doubleValue();
+ if(Emt > totEm/nnhcbins*(bin+1))bin++;
+ }
+ mEm = new double[nnhcbins];
+ mE = new double[nnhcbins];
+ for(int i=0;i<nnhcbins;i++)
+ {
+ int j = nnhcbins-1-i;
+ mEm[j] = sEm[i]/N[i];
+ mE[j] = sE[i]/N[i];
+ }
+ Emstr = "Emeas: ";
+ Ecstr = "Ecorr: ";
+ for(int i=0;i<nnhcbins-1;i++)
+ {
+ Emstr += mEm[i] + ",";
+ Ecstr += mE[i] + ",";
+ }
+ Emstr += mEm[nnhcbins-1];
+ Ecstr += mE[nnhcbins-1];
+ System.out.println(Emstr);
+ System.out.println(Ecstr);
+ Emt = 0.;
+ bin = 0;
+ double[] nbmax = {3.,8.,15.,25.,40.,60.,100.,1000.};
+ for(int i=0;i<nhEm.size();i++)
+ {
+ double Ec;
+ int ib = 0;
+ if(nhEm.get(i).doubleValue() <= mEm[0])
+ {
+ ib = -1;
+ Ec = nhEm.get(i).doubleValue()*mE[0]/mEm[0];
+ }
+ else if(nhEm.get(i).doubleValue() >= mEm[nnhcbins-1])
+ {
+ ib = 19;
+ Ec = nhEm.get(i).doubleValue()*mE[nnhcbins-1]/mEm[nnhcbins-1];
+ }
+ else
+ {
+ for(int j=1;j<nnhcbins;j++)
+ {
+ if(nhEm.get(i).doubleValue() < mEm[j])break;
+ ib++;
+ }
+ Ec = mE[ib] + (nhEm.get(i).doubleValue() - mEm[ib])*(mE[ib+1] - mE[ib])/(mEm[ib+1] - mEm[ib]);
+ }
+ int jb = bin;
+ double Egen = nhE.get(i).doubleValue();
+ double dE = (Ec - Egen);
+ int pbin = 0;
+ for(int k=0;k<nbmax.length;k++)
+ {
+ if(Egen < nbmax[k])break;
+ pbin++;
+ }
+ aida.cloud1D("NH/E",mx).fill(Egen);
+ aida.cloud1D("NH/dE",mx).fill(dE);
+ aida.cloud1D("NH/dE over E",mx).fill(dE/Egen);
+ aida.cloud1D("NH/dE over rootE",mx).fill(dE/Math.sqrt(Egen));
+ aida.cloud1D("NH/wted E",mx).fill(Egen,Egen);
+ aida.cloud1D("NH/wted dE",mx).fill(dE,Egen);
+ aida.cloud1D("NH/wted dE over E",mx).fill(dE/Egen,Egen);
+ aida.cloud1D("NH/wted dE over rootE",mx).fill(dE/Math.sqrt(Egen),Egen);
+ aida.cloud1D("NH/bin "+pbin+" E",mx).fill(Egen);
+ aida.cloud1D("NH/bin "+pbin+" dE",mx).fill(dE);
+ aida.cloud1D("NH/bin "+pbin+" dE over E",mx).fill(dE/Egen);
+ aida.cloud1D("NH/bin "+pbin+" dE over rootE",mx).fill(dE/Math.sqrt(Egen));
+ aida.cloud1D("NH/bin "+pbin+" wted E",mx).fill(Egen,Egen);
+ aida.cloud1D("NH/bin "+pbin+" wted dE",mx).fill(dE,Egen);
+ aida.cloud1D("NH/bin "+pbin+" wted dE over E",mx).fill(dE/Egen,Egen);
+ aida.cloud1D("NH/bin "+pbin+" wted dE over rootE",mx).fill(dE/Math.sqrt(Egen),Egen);
+ Emt += nhEm.get(i).doubleValue();
+ if(Emt > totEm/nnhcbins*(bin+1))bin++;
+ }
+ }
+ public void detectorChanged(Detector detector)
+ {
+ if(dchanged)
+ {
+ System.out.println("Only 1 detector allowed for calibration!! Aborting");
+ System.exit(0);
+ }
+ dchanged = true;
+ Subdetector subdetector = detector.getSubdetector("EMBarrel");
+ List<Integer> changeLayers = new ArrayList<Integer>();
+ double thickness = -1.;
+ for (int i=0, n=subdetector.getLayering().getLayerCount(); i<n; i++)
+ {
+ double thisThickness = subdetector.getLayering().getLayer(i).getThickness();
+ if (i==0)
+ {
+ thickness = thisThickness;
+ }
+
+ if (thisThickness != thickness)
+ {
+ changeLayers.add(i);
+ thickness = thisThickness;
+ }
+ }
+ int maxch = 1;
+ if(changeLayers.get(0) == 1)maxch = 2;
+ if(changeLayers.size() > maxch)
+ {
+ System.out.println("Too many changes in Ecal! Aborting");
+ System.exit(0);
+ }
+ nfrontecallayers = changeLayers.get(maxch-1);
+ }
+}
+class MyDoubleArray
+{
+ double[] x;
+ double E;
+ public MyDoubleArray(double[] a,double e)
+ {
+ x = a;
+ E = e;
+ }
+ public double getVal(int i)
+ {
+ if(i < x.length)return x[i];
+ System.out.println("getVal("+i+"): out of range, returning 0.");
+ return 0.;
+ }
+ public double getE()
+ {
+ return E;
+ }
+}