1
2
3
4
5
6
7
8 package org.lcsim.cal.calib;
9
10 import static java.lang.Math.abs;
11 import static java.lang.Math.sqrt;
12 import hep.aida.ITree;
13 import hep.physics.vec.Hep3Vector;
14
15 import java.util.List;
16
17 import org.lcsim.conditions.ConditionsManager;
18 import org.lcsim.conditions.ConditionsSet;
19 import org.lcsim.conditions.ConditionsManager.ConditionsSetNotFoundException;
20 import org.lcsim.event.CalorimeterHit;
21 import org.lcsim.event.Cluster;
22 import org.lcsim.event.EventHeader;
23 import org.lcsim.event.MCParticle;
24 import org.lcsim.geometry.IDDecoder;
25 import org.lcsim.geometry.Subdetector;
26 import org.lcsim.geometry.subdetector.AbstractPolyhedraCalorimeter;
27 import org.lcsim.geometry.subdetector.PolyhedraEndcapCalorimeter2;
28 import org.lcsim.recon.cluster.fixedcone.FixedConeClusterer;
29 import org.lcsim.recon.cluster.fixedcone.FixedConeClusterer.FixedConeDistanceMetric;
30 import org.lcsim.util.Driver;
31 import org.lcsim.util.aida.AIDA;
32
33 import Jama.Matrix;
34
35
36
37
38
39 public class SamplingFractionAnalysisPolyCalDriver extends Driver
40 {
41
42 private ConditionsSet _cond;
43 private CollectionManager _collectionmanager = CollectionManager.defaultInstance();
44
45
46
47
48
49 private double[][] _acc = new double[3][3];
50 private double[] _vec = new double[3];
51
52
53 private FixedConeClusterer _fcc;
54 private AIDA aida = AIDA.defaultInstance();
55 private ITree _tree;
56 private boolean _initialized = false;
57 private boolean _debug = false;
58
59 boolean skipFirstLayer = false;
60
61
62
63 private double[] _ecalLayering;
64 boolean _useFirstLayer;
65
66
67 double emCalInnerZ = 0.;
68
69 boolean _isHcalDigital = false;
70
71
72
73
74 public SamplingFractionAnalysisPolyCalDriver()
75 {
76 _tree = aida.tree();
77 }
78
79 protected void process(EventHeader event)
80 {
81 super.process(event);
82
83
84 String[] det
85 = {
86 "EcalBarrel", "EcalEndcap"
87 };
88
89
90
91
92
93
94
95 if (!_initialized) {
96 ConditionsManager mgr = ConditionsManager.defaultInstance();
97 try {
98 _cond = mgr.getConditions("CalorimeterCalibration");
99 } catch (ConditionsSetNotFoundException e) {
100 System.out.println("ConditionSet CalorimeterCalibration not found for detector " + mgr.getDetector());
101 System.out.println("Please check that this properties file exists for this detector ");
102 }
103 double radius = .5;
104 double seed = 0.;
105 double minE = .05;
106 _fcc = new FixedConeClusterer(radius, seed, minE, FixedConeDistanceMetric.DPHIDTHETA);
107
108
109
110
111 System.out.println("looking up subdet: " + det[0]);
112 Subdetector subdet = event.getDetector().getSubdetectors().get(det[0]);
113 System.out.println("proc subdet: " + subdet.getName());
114
115 _ecalLayering = _cond.getDoubleArray("ECalLayering");
116 _useFirstLayer = _cond.getDouble("IsFirstEmLayerSampling") == 1.;
117 skipFirstLayer = !_useFirstLayer;
118 ConditionsSet hcalProperties = mgr.getConditions("SamplingFractions/HcalBarrel");
119
120 _isHcalDigital = hcalProperties.getBoolean("digital");
121 System.out.println("HCal is " + (_isHcalDigital == true ? "" : "not") + " read out digitally");
122
123 AbstractPolyhedraCalorimeter calsubBarrel = (AbstractPolyhedraCalorimeter) event.getDetector().getSubdetectors().get(det[0]);
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146 if (skipFirstLayer) {
147 System.out.println("processing " + event.getDetectorName() + " with an em calorimeter with a massless first gap");
148 }
149
150 System.out.println("initialized...");
151
152 _initialized = true;
153 }
154
155
156 List<MCParticle> mcparts = event.getMCParticles();
157
158 MCParticle mcpart = null;
159
160 for (MCParticle myMCP : mcparts) {
161 if (myMCP.getGeneratorStatus() == MCParticle.FINAL_STATE) {
162 if (mcpart == null) {
163 mcpart = myMCP;
164 } else if (mcpart.getEnergy() < myMCP.getEnergy()) {
165 mcpart = myMCP;
166 }
167 }
168 }
169
170 if (mcpart != null) {
171 String particleType = mcpart.getType().getName();
172 double mcEnergy = mcpart.getEnergy();
173 long mcIntegerEnergy = Math.round(mcEnergy);
174 boolean meV = false;
175 if (mcEnergy < .99) {
176 mcIntegerEnergy = Math.round(mcEnergy * 1000);
177 meV = true;
178 }
179
180 _tree.mkdirs(particleType);
181 _tree.cd(particleType);
182 _tree.mkdirs(mcIntegerEnergy + (meV ? "_MeV" : "_GeV"));
183 _tree.cd(mcIntegerEnergy + (meV ? "_MeV" : "_GeV"));
184
185
186
187
188
189
190
191
192
193
194
195
196
197 if (mcpart.getSimulatorStatus().isDecayedInCalorimeter()) {
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228 String processedHitsName = _cond.getString("ProcessedHitsCollectionName");
229 List<CalorimeterHit> hitsToCluster = _collectionmanager.getList(processedHitsName);
230
231 if (_debug) {
232 System.out.println("clustering " + hitsToCluster.size() + " hits");
233 }
234
235
236
237
238
239
240
241 List<Cluster> clusters = _fcc.createClusters(hitsToCluster);
242 if (_debug) {
243 System.out.println("found " + clusters.size() + " clusters");
244 }
245 aida.histogram1D("number of found clusters", 10, -0.5, 9.5).fill(clusters.size());
246 for (Cluster c : clusters) {
247
248 aida.cloud1D("cluster energy for all clusters").fill(c.getEnergy());
249 }
250
251
252
253 if (clusters.size() > 0) {
254 int nHCalHits = 0;
255 Cluster c = clusters.get(0);
256
257 aida.cloud1D("Highest energy cluster energy").fill(c.getEnergy());
258 aida.cloud1D("Highest energy cluster number of cells").fill(c.getCalorimeterHits().size());
259
260 double clusterEnergy = c.getEnergy();
261 double mcMass = mcpart.getType().getMass();
262
263 double expectedEnergy = mcEnergy;
264
265 if (mcpart.getPDGID() == 2212 || mcpart.getPDGID() == 2112) {
266 expectedEnergy = mcEnergy - mcMass;
267 }
268
269 aida.cloud1D("measured - predicted energy").fill(clusterEnergy - expectedEnergy);
270
271
272
273
274
275
276 List<CalorimeterHit> hits = c.getCalorimeterHits();
277 double[] vals = new double[4];
278 int nHcalHits = 0;
279 double clusterRawEnergy = 0.;
280 for (CalorimeterHit hit : hits) {
281 long id = hit.getCellID();
282 IDDecoder decoder = hit.getIDDecoder();
283 decoder.setID(id);
284 int layer = decoder.getLayer();
285 String detectorName = decoder.getSubdetector().getName();
286 int type = 0;
287 int caltype = 0;
288
289 if (detectorName.toUpperCase().startsWith("ECAL")) {
290
291
292
293
294
295 for (int i = 1; i < _ecalLayering.length + 1; ++i) {
296 if (layer >= _ecalLayering[i - 1]) {
297 caltype = i - 1;
298 }
299 }
300 }
301
302 if (detectorName.toUpperCase().startsWith("HCAL")) {
303 type = 2;
304 caltype = 3;
305 nHCalHits += 1;
306 }
307 if(_debug)
308 {
309 System.out.println(detectorName+" layer: "+layer+" type: "+type+" caltype: "+caltype+" raw");
310 }
311 clusterRawEnergy += hit.getRawEnergy();
312 vals[caltype] += hit.getRawEnergy();
313
314
315
316 }
317 if (_isHcalDigital == true) {
318 vals[3] = nHCalHits;
319 }
320
321
322 for (int j = 0; j < 3; ++j) {
323 if(_debug)
324 {
325 System.out.println("clusterRawEnergy= "+clusterRawEnergy+" vals["+j+"]= "+vals[j]);
326 }
327 _vec[j] += expectedEnergy*vals[j+1];
328 for (int k = 0; k < 3; ++k) {
329 _acc[j][k] += vals[j+1] * vals[k+1] ;
330 }
331 }
332 }
333
334
335
336 event.put("Found Clusters", clusters);
337 }
338 } else {
339 System.out.println("null MCPointer at event " + event.getEventNumber());
340 }
341 _tree.cd("/");
342 }
343
344 protected void endOfData()
345 {
346 System.out.println("done! endOfData.");
347
348 Matrix A = new Matrix(_acc, 3, 3);
349 A.print(6, 4);
350 Matrix b = new Matrix(3, 1);
351 for (int i = 0; i < 3; ++i) {
352 b.set(i, 0, _vec[i]);
353 }
354 b.print(6, 4);
355 try {
356 Matrix x = A.solve(b);
357 x.print(6, 4);
358 System.out.println("SamplingFractions: " + (skipFirstLayer ? "1., " : "") + 1. / x.get(0, 0) + ", " + 1. / x.get(1, 0) + ", " + 1. / x.get(2, 0));
359 } catch (Exception e) {
360 e.printStackTrace();
361 System.out.println("try reducing dimensionality...");
362 Matrix Ap = new Matrix(_acc, 2, 2);
363 Ap.print(6, 4);
364 Matrix bp = new Matrix(2, 1);
365 for (int i = 0; i < 2; ++i) {
366 bp.set(i, 0, _vec[i]);
367 }
368 bp.print(6, 4);
369 try {
370 Matrix x = Ap.solve(bp);
371 x.print(6, 4);
372 } catch (Exception e2) {
373 e2.printStackTrace();
374 }
375 }
376 }
377
378 public void setDebug(boolean debug)
379 {
380 _debug = debug;
381 }
382 }