1 /*
2 * SamplingFractionAnalysisDriver.java
3 *
4 * Created on May 19, 2008, 11:54 AM
5 *
6 * $Id: SamplingFractionAnalysisDriver.java,v 1.8 2010/05/28 18:34:00 ngraf Exp $
7 */
8
9 package org.lcsim.cal.calib;
10
11 import Jama.Matrix;
12 import hep.aida.ITree;
13 import hep.physics.vec.Hep3Vector;
14 import java.util.List;
15 import org.lcsim.event.CalorimeterHit;
16 import org.lcsim.event.Cluster;
17 import org.lcsim.event.EventHeader;
18 import org.lcsim.event.MCParticle;
19 import org.lcsim.geometry.IDDecoder;
20 import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
21 import org.lcsim.geometry.subdetector.CylindricalEndcapCalorimeter;
22 import org.lcsim.recon.cluster.fixedcone.FixedConeClusterer;
23 import org.lcsim.recon.cluster.fixedcone.FixedConeClusterer.FixedConeDistanceMetric;
24 import org.lcsim.util.Driver;
25 import org.lcsim.util.aida.AIDA;
26
27 import static java.lang.Math.abs;
28 import static java.lang.Math.sqrt;
29 import org.lcsim.conditions.ConditionsManager;
30 import org.lcsim.conditions.ConditionsManager.ConditionsSetNotFoundException;
31 import org.lcsim.conditions.ConditionsSet;
32
33 /**
34 *
35 * @author Norman Graf
36 */
37 public class SamplingFractionAnalysisDriver extends Driver
38 {
39 private ConditionsSet _cond;
40 private CollectionManager _collectionmanager = CollectionManager.defaultInstance();
41 // we will accumulate the raw energy values in three depths:
42 // 1. Layers (0)1 through (20)21 of the EM calorimeter (note that if layer 0 is massless, SF==1.)
43 // 2. last ten layers of the EM calorimeter
44 // 3. the hadron calorimeter
45 //
46 private double[][] _acc = new double[3][3];
47 private double[] _vec = new double[3];
48
49 // let's use a clusterer to remove effects of calorimeter cells hit far, far away.
50 // use the only cross-detector clusterer we have:
51 private FixedConeClusterer _fcc;
52
53 private AIDA aida = AIDA.defaultInstance();
54 private ITree _tree;
55
56 private boolean _initialized = false;
57 private boolean _debug = false;
58
59
60 // TODO fix this dependence on EM calorimeter geometry
61 boolean skipFirstLayer = false;
62 int firstEmStartLayer = 0;
63 int secondEmStartLayer = 20;
64
65 double emCalInnerRadius = 0.;
66 double emCalInnerZ = 0.;
67
68 /** Creates a new instance of SamplingFractionAnalysisDriver */
69 public SamplingFractionAnalysisDriver()
70 {
71 _tree = aida.tree();
72 }
73
74
75 protected void process(EventHeader event)
76 {
77 super.process(event);
78 //System.out.println("processing SamplingFractionAnalysisDriver");
79 // TODO make these values runtime definable
80 String[] det = {"EMBarrel","EMEndcap"};
81 // String[] collNames = {"EcalBarrHits", "EcalEndcapHits", "HcalBarrHits", "HcalEndcapHits"};
82 // double[] mipHistMaxBinVal = {.0005, .0005, .005, .005};
83 // double timeCut = 100.; // cut on energy depositions coming more than 100 ns late
84 // double ECalMipCut = .0001/3.; // determined from silicon using muons at normal incidence
85 // double HCalMipCut = .0008/3.; // determined from scintillator using muons at normal incidence
86 // double[] mipCut = {ECalMipCut, ECalMipCut, HCalMipCut, HCalMipCut};
87
88
89 if(!_initialized)
90 {
91 ConditionsManager mgr = ConditionsManager.defaultInstance();
92 try
93 {
94 _cond = mgr.getConditions("CalorimeterCalibration");
95 }
96 catch(ConditionsSetNotFoundException e)
97 {
98 System.out.println("ConditionSet CalorimeterCalibration not found for detector "+mgr.getDetector());
99 System.out.println("Please check that this properties file exists for this detector ");
100 }
101 double radius = .5;
102 double seed = 0.;//.1;
103 double minE = .05; //.25;
104 _fcc = new FixedConeClusterer(radius, seed, minE, FixedConeDistanceMetric.DPHIDTHETA);
105
106 // detector geometries here...
107 // barrel
108 CylindricalCalorimeter calsubBarrel = (CylindricalCalorimeter)event.getDetector().getSubdetectors().get(det[0]);
109 // TODO remove this hardcoded dependence on the first layer
110 if(calsubBarrel.getLayering().getLayer(0).getSlices().get(0).isSensitive())
111 {
112 skipFirstLayer = true;
113 firstEmStartLayer += 1;
114 secondEmStartLayer += 1;
115 }
116 // Layering layering = calsubBarrel.getLayering();
117 // for(int i=0; i<layering.size(); ++i)
118 // {
119 // Layer l = layering.getLayer(i);
120 // System.out.println("layering "+i);
121 // List<LayerSlice> slices = l.getSlices();
122 // for(int j=0; j<slices.size(); ++j)
123 // {
124 // LayerSlice slice = slices.get(j);
125 // System.out.println("Layer "+i+" slice "+j+" is "+ slice.getMaterial().getName() +" and "+(slice.isSensitive() ? " is sensitive" : ""));
126 // }
127 // }
128 emCalInnerRadius = calsubBarrel.getInnerRadius();
129 //endcap
130 CylindricalEndcapCalorimeter calsubEndcap = (CylindricalEndcapCalorimeter)event.getDetector().getSubdetectors().get(det[1]);
131 emCalInnerZ = abs(calsubEndcap.getZMin());
132 if(skipFirstLayer) System.out.println("processing "+event.getDetectorName()+" with an em calorimeter with a massless first gap");
133 System.out.println("Calorimeter bounds: r= "+emCalInnerRadius+ " z= "+emCalInnerZ);
134 System.out.println("initialized...");
135
136 _initialized = true;
137 }
138
139 // organize the histogram tree by species and energy
140 List<MCParticle> mcparts = event.getMCParticles();
141 MCParticle mcpart = mcparts.get(mcparts.size()-1);
142 String particleType = mcpart.getType().getName();
143 double mcEnergy = mcpart.getEnergy();
144 long mcIntegerEnergy = Math.round(mcEnergy);
145 boolean meV = false;
146 if(mcEnergy<.99)
147 {
148 mcIntegerEnergy = Math.round(mcEnergy*1000);
149 meV = true;
150 }
151
152 _tree.mkdirs(particleType);
153 _tree.cd(particleType);
154 _tree.mkdirs(mcIntegerEnergy+(meV ? "_MeV": "_GeV"));
155 _tree.cd(mcIntegerEnergy+(meV ? "_MeV": "_GeV"));
156
157 // this analysis is intended for single particle calorimeter response.
158 // let's make sure that the primary particle did not interact in the
159 // tracker...
160 Hep3Vector endpoint = mcpart.getEndPoint();
161 // this is just crap. Why not use SpacePoint?
162 double radius = sqrt(endpoint.x()*endpoint.x()+endpoint.y()*endpoint.y());
163 double z = endpoint.z();
164 // System.out.println("Input MCParticle endpoint: r="+radius+" z= "+z);
165
166 boolean doit = true;
167 if(radius<emCalInnerRadius && abs(z) < emCalInnerZ) doit = false;
168 if(doit)
169 {
170 // // now let's check the em calorimeters...
171 // // get all of the calorimeter hits...
172 // List<CalorimeterHit> allHits = new ArrayList<CalorimeterHit>();
173 // // and the list after cuts.
174 // List<CalorimeterHit> hitsToCluster = new ArrayList<CalorimeterHit>();
175 // int i = 0;
176 // for(String name : collNames)
177 // {
178 //// System.out.println("fetching "+name+" from the event");
179 // List<CalorimeterHit> hits = event.get(CalorimeterHit.class, name);
180 //// System.out.println(name+ " has "+hits.size()+" hits");
181 // // let's look at the hits and see if we need to cut on energy or time...
182 // for(CalorimeterHit hit: hits)
183 // {
184 // aida.histogram1D(name+" raw calorimeter cell energy",100, 0., mipHistMaxBinVal[i]).fill(hit.getRawEnergy());
185 // aida.histogram1D(name+" raw calorimeter cell energy full range",100, 0., 0.2).fill(hit.getRawEnergy());
186 //// aida.cloud1D(name+" raw calorimeter cell energies").fill(hit.getRawEnergy());
187 // aida.histogram1D(name+" calorimeter cell time",100,0., 200.).fill(hit.getTime());
188 // if(hit.getTime()<timeCut)
189 // {
190 // if(hit.getRawEnergy()>mipCut[i])
191 // {
192 // hitsToCluster.add(hit);
193 // }
194 // }
195 // }
196 // allHits.addAll(hits);
197 // ++i;
198 // }
199 // System.out.println("ready to cluster "+hitsToCluster.size()+ " hits");
200 String processedHitsName = _cond.getString("ProcessedHitsCollectionName");
201 List<CalorimeterHit> hitsToCluster = _collectionmanager.getList(processedHitsName);//event.get(CalorimeterHit.class, processedHitsName);
202
203 if(_debug) System.out.println("clustering "+hitsToCluster.size()+" hits");
204 // quick check
205 // for(CalorimeterHit hit : hitsToCluster)
206 // {
207 // System.out.println("hit ");
208 // System.out.println(hit.getLCMetaData().getName());
209 // }
210 // cluster the hits
211 List<Cluster> clusters = _fcc.createClusters(hitsToCluster);
212 if(_debug) System.out.println("found "+clusters.size()+" clusters");
213 aida.histogram1D("number of found clusters", 10, -0.5, 9.5).fill(clusters.size());
214 for(Cluster c : clusters)
215 {
216 // System.out.println(c);
217 aida.cloud1D("cluster energy for all clusters").fill(c.getEnergy());
218 }
219
220 // proceed only if we found a single cluster above threshold
221 // too restrictive! simply take the highest energy cluster
222 if(clusters.size()>0)
223 {
224 Cluster c = clusters.get(0);
225
226 aida.cloud1D("Highest energy cluster energy").fill(c.getEnergy());
227 aida.cloud1D("Highest energy cluster number of cells").fill(c.getCalorimeterHits().size());
228
229 double clusterEnergy = c.getEnergy();
230 double mcMass = mcpart.getType().getMass();
231 // subtract the mass to get kinetic energy...
232 double expectedEnergy = sqrt(mcEnergy*mcEnergy - mcMass*mcMass);
233 // System.out.println(mcpart.getType().getName()+" "+expectedEnergy);
234 aida.cloud1D("measured - predicted energy").fill(clusterEnergy - expectedEnergy);
235
236 // let's now break down the cluster by component.
237 // this analysis uses:
238 // 1.) first 20 EM layers
239 // 2.) next 10 EM layers
240 // 3.) Had layers
241 List<CalorimeterHit> hits = c.getCalorimeterHits();
242 double[] vals = new double[3];
243 double clusterRawEnergy = 0.;
244 for(CalorimeterHit hit : hits)
245 {
246 long id = hit.getCellID();
247 IDDecoder decoder = hit.getIDDecoder();
248 decoder.setID(id);
249 int layer = decoder.getLayer();
250 String detectorName = decoder.getSubdetector().getName();
251 int type = 0;
252 if(detectorName.startsWith("EM"))
253 {
254 if(layer>=firstEmStartLayer && layer<secondEmStartLayer)
255 {
256 type = 0;
257 }
258 else
259 {
260 type = 1;
261 }
262 }
263 if(detectorName.startsWith("HAD"))
264 {
265 type = 2;
266 }
267 clusterRawEnergy += hit.getRawEnergy();
268 vals[type] += hit.getRawEnergy();
269 } // end of loop over hits in cluster
270 // set up linear least squares:
271 // expectedEnergy = a*E1 + b*E2 +c*E3
272 for(int j=0; j<3; ++j)
273 {
274 _vec[j] += expectedEnergy*vals[j];
275 for(int k=0; k<3; ++k)
276 {
277 _acc[j][k] += vals[j]*vals[k];
278 }
279 }
280 } // end of single cluster cut
281
282 // event.put("All Calorimeter Hits",allHits);
283 // event.put("Hits To Cluster",hitsToCluster);
284 event.put("Found Clusters",clusters);
285 }// end of check on decays outside tracker volume
286 _tree.cd("/");
287 }
288
289 protected void endOfData()
290 {
291 System.out.println("done! endOfData.");
292 // calculate the sampling fractions...
293 Matrix A = new Matrix(_acc, 3, 3);
294 A.print(6,4);
295 Matrix b = new Matrix(3,1);
296 for(int i=0; i<3; ++i)
297 {
298 b.set(i,0,_vec[i]);
299 }
300 b.print(6,4);
301 try
302 {
303 Matrix x = A.solve(b);
304 x.print(6, 4);
305 System.out.println("SamplingFractions: "+ (skipFirstLayer ? "1., ":"")+1./x.get(0,0)+ ", " + 1./x.get(1,0)+", "+1./x.get(2,0));
306 }
307 catch(Exception e)
308 {
309 e.printStackTrace();
310 System.out.println("try reducing dimensionality...");
311 Matrix Ap = new Matrix(_acc, 2, 2);
312 Ap.print(6,4);
313 Matrix bp = new Matrix(2,1);
314 for(int i=0; i<2; ++i)
315 {
316 bp.set(i,0,_vec[i]);
317 }
318 bp.print(6,4);
319 try
320 {
321 Matrix x = Ap.solve(bp);
322 x.print(6, 4);
323 }
324 catch(Exception e2)
325 {
326 e2.printStackTrace();
327 }
328 }
329 }
330 }