lcsim/src/org/lcsim/recon/cheater
diff -u -r1.1 -r1.2
--- QSFCalibrationFromData.java 12 Dec 2008 21:39:01 -0000 1.1
+++ QSFCalibrationFromData.java 7 May 2010 19:00:56 -0000 1.2
@@ -1,444 +1 @@
-package org.lcsim.recon.cheater;
-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);
- }
-}
+package org.lcsim.recon.cheater;
import java.util.*;
import org.lcsim.geometry.Detector;
import org.lcsim.geometry.Subdetector;
import org.lcsim.util.Driver;
import org.lcsim.event.*;
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.*;
import org.lcsim.recon.util.CalorimeterInformation;
import org.lcsim.geometry.Calorimeter.CalorimeterType;
import org.lcsim.recon.util.CalInfoDriver;
/**
* 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
//
//
// 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;
CalorimeterInformation ci;
public QSFCalibrationFromData()
{
add(new CalInfoDriver());
//
// Always use digisim
//
add(new DigiPackageDriver());
//
// initialization
//
dchanged = false;
A = ne!
w 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);
if(ci == null)ci = CalorimeterInformation.instance();
String[] calcoll = {ci.getDigiCollectionName(CalorimeterType.EM_BARREL),
ci.getDigiCollectionName(CalorimeterType.EM_ENDCAP),
ci.getDigiCollectionName(CalorimeterType.HAD_BARREL),
ci.getDigiCollectionName(CalorimeterType.HAD_ENDCAP)};
//
// 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[ro!
w][col] += Esums[row] * Esums[col] / ZE;
i!
f ((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 (i!
nt ic = 0; ic < nct; ic++)
{
r9!
0 = 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 calibrati!
on!! Aborting");
System.exit(0);
}
!
dchange
d = 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);
}
}
\ No newline at end of file
