lcsim/src/org/lcsim/contrib/Cassell/recon/Cheat
diff -u -r1.1 -r1.2
--- QCalibrationFromData.java 8 Sep 2008 23:11:30 -0000 1.1
+++ QCalibrationFromData.java 16 Nov 2008 02:44:06 -0000 1.2
@@ -1,544 +1,582 @@
-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;
- }
-}
+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.};
+ aida.cloud1D("Photon/E",mx).reset();
+ aida.cloud1D("Photon/dE",mx).reset();
+ aida.cloud1D("Photon/dE over E",mx).reset();
+ aida.cloud1D("Photon/dE over rootE",mx).reset();
+ aida.cloud1D("Photon/wted E",mx).reset();
+ aida.cloud1D("Photon/wted dE",mx).reset();
+ aida.cloud1D("Photon/wted dE over E",mx).reset();
+ aida.cloud1D("Photon/wted dE over rootE",mx).reset();
+ for(int pbin=0;pbin<pbmax.length+1;pbin++)
+ {
+ aida.cloud1D("Photon/bin "+pbin+" dE",mx).reset();
+ aida.cloud1D("Photon/bin "+pbin+" E",mx).reset();
+ aida.cloud1D("Photon/bin "+pbin+" dE over E",mx).reset();
+ aida.cloud1D("Photon/bin "+pbin+" dE over rootE",mx).reset();
+ aida.cloud1D("Photon/bin "+pbin+" wted E",mx).reset();
+ aida.cloud1D("Photon/bin "+pbin+" wted dE",mx).reset();
+ aida.cloud1D("Photon/bin "+pbin+" wted dE over E",mx).reset();
+ aida.cloud1D("Photon/bin "+pbin+" wted dE over rootE",mx).reset();
+ }
+ 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);
[truncated at 1000 lines; 130 more skipped]
