lcsim/src/org/lcsim/contrib/Cassell/recon/Cheat
diff -N QSFCalibrationFromData.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ QSFCalibrationFromData.java 8 Sep 2008 23:10:44 -0000 1.1
@@ -0,0 +1,444 @@
+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;
+import org.lcsim.digisim.*;
+import hep.aida.*;
+
+/**
+ * Calculate sampling fractions from fixed E
+ * physics data minimizing event dE/sqrtE
+ * Ron Cassell
+ */
+public class QSFCalibrationFromData extends Driver
+{
+ private AIDA aida = AIDA.defaultInstance();
+ IDDecoder _decoder;
+ int ievt;
+ int evtused;
+ int nevtmax = 40000;
+ int mx = 2000000;
+//
+// # of calorimeters for energy sums
+//
+ int nsums = 6;
+ int nsumsem = 7;
+ int nsumsm = 8;
+ double[][] A;
+ double[][] B;
+ double[][] Am;
+ double[][] Bm;
+ double[][] Aem;
+ double[][] Bem;
+//
+// Calorimeter collextions to use
+//
+ String[] calcoll = {"EcalBarrDigiHits","EcalEndcapDigiHits","HcalBarrDigiHits","HcalEndcapDigiHits","MuonEndcapDigiHits","MuonBarrDigiHits"};
+//
+// 2 sums per Ecal collection
+//
+ int[] cind = {0,2,4,5,6,7};
+//
+// digital flag per calorimeter collection
+//
+ boolean[] digital = {false,false,true,true,true,true};
+//
+// where the Ecal changes, overridden by detector changed method
+//
+ int nfrontecallayers = 21;
+//
+// Energy sums per event
+//
+ double[][] es;
+ RonRMS90Calculator calc;
+//
+// Generated Energy and cos theta per event
+//
+ double[] zE;
+ double[] ctheta;
+ boolean first;
+//
+// Keep track of the different cm Energies for plots
+//
+ int nE;
+ int[] Eevt;
+//
+// Only allow 1 detector
+//
+ boolean dchanged;
+ public QSFCalibrationFromData()
+ {
+//
+// Always use digisim
+//
+ add(new DigiPackageDriver());
+//
+// initialization
+//
+ dchanged = false;
+ A = new double[nsums][nsums];
+ B = new double[nsums][1];
+ Am = new double[nsumsm][nsumsm];
+ Bm = new double[nsumsm][1];
+ Aem = new double[nsumsem][nsumsem];
+ Bem = new double[nsumsem][1];
+ ievt = 0;
+ evtused = 0;
+ es = new double[nevtmax][nsumsm];
+ calc = new RonRMS90Calculator();
+ zE = new double[nevtmax];
+ ctheta = new double[nevtmax];
+ first = true;
+ nE = 0;
+ Eevt = new int[10];
+ }
+ protected void process(EventHeader event)
+ {
+//
+// Find the generated Energy and cos theta
+//
+ List<MCParticle> mcl = event.get(MCParticle.class,"MCParticle");
+ double Zmass = 0.;
+ double ZE = 0.;
+ double cost = 0.;
+ for(MCParticle p:mcl)
+ {
+ int id = Math.abs(p.getPDGID());
+ if( (id == 1 )||(id == 2 )||(id == 3 ) )
+ {
+ if(p.getParents().get(0).getPDGID() == 23 )
+ {
+ Zmass = p.getParents().get(0).getMass();
+ ZE = p.getParents().get(0).getEnergy();
+ Hep3Vector pp = p.getMomentum();
+ cost = Math.abs(pp.z()/pp.magnitude());
+ }
+ }
+ }
+//
+// Initialize cmE on first event
+//
+ if(first)
+ {
+ first = false;
+ int El = (int) (ZE + .5);
+ Eevt[0] = El;
+ nE = 1;
+ }
+ System.out.println(" Processing event "+ievt);
+ super.process(event);
+//
+// See if cmE changed
+//
+ int Ec = (int) (ZE + .5);
+ boolean match = false;
+ for(int i=0;i<nE;i++)
+ {
+ if(Ec == Eevt[i])
+ {
+ match = true;
+ break;
+ }
+ }
+ if(!match)
+ {
+ Eevt[nE] = Ec;
+ nE++;
+ }
+//
+// Cache generated energy and cos theta per event
+//
+ zE[ievt] = ZE;
+ ctheta[ievt] = cost;
+//
+// Calculate the energy sums for minimization
+//
+ double[] Esums = new double[nsumsm];
+ for(int i=0;i<calcoll.length;i++)
+ {
+ List<CalorimeterHit> hitlist = event.get(CalorimeterHit.class,calcoll[i]);
+ int ind = cind[i];
+ for(CalorimeterHit hit:hitlist)
+ {
+ int subind = 0;
+ if(i < 2)
+ {
+ IDDecoder idc = hit.getIDDecoder();
+ idc.setID(hit.getCellID());
+ int layer = idc.getValue("layer");
+ if(layer >= nfrontecallayers)subind = 1;
+ Esums[ind+subind] += hit.getRawEnergy();
+ }
+ else
+ {
+ if(digital[i])Esums[ind] += 1.;
+ else Esums[ind] += hit.getRawEnergy();
+ }
+ }
+ }
+//
+// Store the appropriate sums with a cos theta cut of 0.9
+//
+ for(int row=0;row<nsumsm;row++)
+ {
+ es[ievt][row] = Esums[row];
+ if(cost < .9)
+ {
+ Bm[row][0] += Esums[row];
+ if(row < nsums)B[row][0] += Esums[row];
+ if(row < nsumsem)Bem[row][0] += Esums[row];
+ for(int col=0;col<nsumsm;col++)
+ {
+ Am[row][col] += Esums[row]*Esums[col]/ZE;
+ if( (row < nsums)&&(col < nsums) ) A[row][col] += Esums[row]*Esums[col]/ZE;
+ if( (row < nsumsem)&&(col < nsumsem) ) Aem[row][col] += Esums[row]*Esums[col]/ZE;
+ }
+ evtused++;
+ }
+ }
+ ievt ++;
+ }
+ protected void suspend()
+ {
+//
+// Calculate the sampling fractions
+//
+ Matrix mA = new Matrix(A);
+ Matrix mB = new Matrix(B);
+ Matrix x = mA.solve(mB);
+ double[] sf = new double[nsums];
+ for(int i=0;i<nsums;i++)
+ {
+ sf[i] = 1./x.get(i,0);
+ }
+ mA = new Matrix(Am);
+ mB = new Matrix(Bm);
+ x = mA.solve(mB);
+ double[] sfm = new double[nsumsm];
+ for(int i=0;i<nsumsm;i++)
+ {
+ sfm[i] = 1./x.get(i,0);
+ }
+ mA = new Matrix(Aem);
+ mB = new Matrix(Bem);
+ x = mA.solve(mB);
+ double[] sfem = new double[nsumsem];
+ for(int i=0;i<nsumsem;i++)
+ {
+ sfem[i] = 1./x.get(i,0);
+ }
+//
+// Make some plots
+//
+ int nct = 12;
+ double[] ctlim = {.1,.2,.3,.4,.5,.6,.7,.8,.9,.95,.98,1.};
+ ICloud1D[][] edist = new ICloud1D[nE][3];
+ ICloud2D[][] rms90c = new ICloud2D[nE][3];
+ ICloud1D[][][] edistct = new ICloud1D[nE][3][nct];
+ ICloud1D[] mfrac = new ICloud1D[nE];
+ for(int iE=0;iE<nE;iE++)
+ {
+ String pre = "evtE="+Eevt[iE]+"GeV/";
+ edist[iE][0] = aida.cloud1D(pre+"no mu/Event dE",mx);
+ edist[iE][0].reset();
+ rms90c[iE][0] = aida.cloud2D(pre+"no mu/rms90 vs ct",mx);
+ rms90c[iE][0].reset();
+ edist[iE][1] = aida.cloud1D(pre+"both mu/Event dE",mx);
+ edist[iE][1].reset();
+ rms90c[iE][1] = aida.cloud2D(pre+"both mu/rms90 vs ct",mx);
+ rms90c[iE][1].reset();
+ edist[iE][2] = aida.cloud1D(pre+"ec mu/Event dE",mx);
+ edist[iE][2].reset();
+ rms90c[iE][2] = aida.cloud2D(pre+"ec mu/rms90 vs ct",mx);
+ rms90c[iE][2].reset();
+ mfrac[iE] = aida.cloud1D(pre+"Fraction of calE in mucal",mx);
+ mfrac[iE].reset();
+ double ll = 0.;
+ for(int ic=0;ic<nct;ic++)
+ {
+ if(ic > 0)ll = ctlim[ic-1];
+ String cst = "dE:"+ll+"<costheta<"+ctlim[ic];
+ edistct[iE][0][ic] = aida.cloud1D(pre+"no mu/"+cst,mx);
+ edistct[iE][0][ic].reset();
+ edistct[iE][1][ic] = aida.cloud1D(pre+"both mu/"+cst,mx);
+ edistct[iE][1][ic].reset();
+ edistct[iE][2][ic] = aida.cloud1D(pre+"ec mu/"+cst,mx);
+ edistct[iE][2][ic].reset();
+ }
+ }
+ for(int j=0;j<ievt;j++)
+ {
+ double evtE = 0;
+ double evtEm = 0;
+ double evtEem = 0;
+ double muE = 0.;
+ int iE = -1;
+ int Eint = (int) (zE[j] + .5);
+ for(int k=0;k<nE;k++)
+ {
+ if(Eint == Eevt[k])
+ {
+ iE = k;
+ break;
+ }
+ }
+ if(iE < 0)
+ {
+ System.out.println("Oops: index < 0, E = "+Eint);
+ iE = 0;
+ }
+ int ic = 0;
+ for(int k=0;k<nct;k++)
+ {
+ if(ctheta[j] <= ctlim[k])break;
+ ic++;
+ }
+ for(int i=0;i<nsumsm;i++)
+ {
+ evtEm += es[j][i]/sfm[i];
+ if(i < nsums)
+ {
+ evtE += es[j][i]/sf[i];
+ }
+ else
+ {
+ muE += es[j][i]/sfm[i];
+ }
+ if(i < nsumsem)
+ {
+ evtEem += es[j][i]/sfem[i];
+ }
+ }
+ if(ctheta[j] < .9)
+ {
+ edist[iE][0].fill(evtE-zE[j]);
+ edist[iE][1].fill(evtEm-zE[j]);
+ edist[iE][2].fill(evtEem-zE[j]);
+ }
+ edistct[iE][0][ic].fill(evtE-zE[j]);
+ edistct[iE][1][ic].fill(evtEm-zE[j]);
+ edistct[iE][2][ic].fill(evtEem-zE[j]);
+ mfrac[iE].fill(muE/evtEm);
+ }
+//
+// Write out the sampling fractions
+//
+ System.out.println("# front ecal layers = "+nfrontecallayers);
+ System.out.println(" SF order: EMB front, EMB back, EME front, EME back, HADB, HADE, MuonE, MuonB");
+ String sfout = "SFs: ";
+ for(int i=0;i<nsums-1;i++)
+ {
+ sfout += sf[i]+",";
+ }
+ sfout += sf[nsums-1];
+ System.out.println(sfout);
+ String sfoutm = "SFs with mucal: ";
+ for(int i=0;i<nsumsm-1;i++)
+ {
+ sfoutm += sfm[i]+",";
+ }
+ sfoutm += sfm[nsumsm-1];
+ System.out.println(sfoutm);
+ String sfoutem = "SFs with endcapmucal: ";
+ for(int i=0;i<nsumsem-1;i++)
+ {
+ sfoutem += sfem[i]+",";
+ }
+ sfoutem += sfem[nsumsem-1];
+ System.out.println(sfoutem);
+//
+// Print out some plot statistics using rms90
+//
+ for(int iE=0;iE<nE;iE++)
+ {
+ System.out.println("event E = "+Eevt[iE]);
+ System.out.println(" no mu");
+ double r90 = calc.calculateRMS90(edist[iE][0]);
+ double m90 = calc.getMEAN90();
+ System.out.println(" ct<.9: m90 = "+m90+": r90 = "+r90);
+ double ll = 0.;
+ for(int ic=0;ic<nct;ic++)
+ {
+ r90 = calc.calculateRMS90(edistct[iE][0][ic]);
+ m90 = calc.getMEAN90();
+ if(ic > 0)ll = ctlim[ic-1];
+ double mct = (ll + ctlim[ic])/2.;
+ rms90c[iE][0].fill(mct,r90);
+ System.out.println(" ct="+mct+": m90 = "+m90+": r90 = "+r90);
+ }
+ System.out.println(" both mu");
+ r90 = calc.calculateRMS90(edist[iE][1]);
+ m90 = calc.getMEAN90();
+ System.out.println(" ct<.9: m90 = "+m90+": r90 = "+r90);
+ ll = 0.;
+ for(int ic=0;ic<nct;ic++)
+ {
+ r90 = calc.calculateRMS90(edistct[iE][1][ic]);
+ m90 = calc.getMEAN90();
+ if(ic > 0)ll = ctlim[ic-1];
+ double mct = (ll + ctlim[ic])/2.;
+ rms90c[iE][1].fill(mct,r90);
+ System.out.println(" ct="+mct+": m90 = "+m90+": r90 = "+r90);
+ }
+ System.out.println(" ec mu");
+ r90 = calc.calculateRMS90(edist[iE][2]);
+ m90 = calc.getMEAN90();
+ System.out.println(" ct<.9: m90 = "+m90+": r90 = "+r90);
+ ll = 0.;
+ for(int ic=0;ic<nct;ic++)
+ {
+ r90 = calc.calculateRMS90(edistct[iE][2][ic]);
+ m90 = calc.getMEAN90();
+ if(ic > 0)ll = ctlim[ic-1];
+ double mct = (ll + ctlim[ic])/2.;
+ rms90c[iE][2].fill(mct,r90);
+ System.out.println(" ct="+mct+": m90 = "+m90+": r90 = "+r90);
+ }
+ }
+ }
+ 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);
+ }
+}