View Javadoc

1   /*
2    * SamplingFractionAnalysisDriver.java
3    *
4    * Created on May 19, 2008, 11:54 AM
5    *
6    * $Id: SamplingFractionAnalysisPolyCalDriver.java,v 1.4 2012/02/10 15:26:26 grefe Exp $
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   * @author Norman Graf
38   */
39  public class SamplingFractionAnalysisPolyCalDriver extends Driver
40  {
41  
42      private ConditionsSet _cond;
43      private CollectionManager _collectionmanager = CollectionManager.defaultInstance();
44      // we will accumulate the raw energy values in three depths:
45      // 1. Layers (0)1 through (20)21 of the EM calorimeter (note that if layer 0 is massless, SF==1.)
46      // 2. last ten layers of the EM calorimeter
47      // 3. the hadron calorimeter
48      //
49      private double[][] _acc = new double[3][3];
50      private double[] _vec = new double[3];
51      // let's use a clusterer to remove effects of calorimeter cells hit far, far away.
52      // use the only cross-detector clusterer we have:
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      // TODO fix this dependence on EM calorimeter geometry
59      boolean skipFirstLayer = false;
60  //    int firstEmStartLayer = 0;
61  //    int secondEmStartLayer = 20;
62  
63      private double[] _ecalLayering;
64      boolean _useFirstLayer;
65  
66  //    double emCalInnerRadius = 0.;
67      double emCalInnerZ = 0.;
68  
69      boolean _isHcalDigital = false;
70  
71      /**
72       * Creates a new instance of SamplingFractionAnalysisDriver
73       */
74      public SamplingFractionAnalysisPolyCalDriver()
75      {
76          _tree = aida.tree();
77      }
78  
79      protected void process(EventHeader event)
80      {
81          super.process(event);
82          //System.out.println("processing SamplingFractionAnalysisDriver");
83          // TODO make these values runtime definable        
84          String[] det
85                  = {
86                      "EcalBarrel", "EcalEndcap"
87                  };
88  //        String[] collNames = {"EcalBarrHits", "EcalEndcapHits", "HcalBarrHits", "HcalEndcapHits"};
89  //        double[] mipHistMaxBinVal = {.0005, .0005, .005, .005};
90  //        double timeCut = 100.; // cut on energy depositions coming more than 100 ns late
91  //        double ECalMipCut = .0001/3.; // determined from silicon using muons at normal incidence
92  //        double HCalMipCut = .0008/3.; // determined from scintillator using muons at normal incidence
93  //        double[] mipCut = {ECalMipCut, ECalMipCut, HCalMipCut, HCalMipCut};
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.;//.1;
105             double minE = .05; //.25;
106             _fcc = new FixedConeClusterer(radius, seed, minE, FixedConeDistanceMetric.DPHIDTHETA);
107 
108             // detector geometries here...
109             // barrel
110             //CylindricalCalorimeter calsubBarrel = (CylindricalCalorimeter)event.getDetector().getSubdetectors().get(det[0]);
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             // TODO remove this hardcoded dependence on the first layer
125 //            if (calsubBarrel.getLayering().getLayer(0).getSlices().get(0).isSensitive()) {
126 //                skipFirstLayer = true;
127 //                firstEmStartLayer += 1;
128 //                secondEmStartLayer += 1;
129 //            }
130 //            Layering layering = calsubBarrel.getLayering();
131 //            for(int i=0; i<layering.size(); ++i)
132 //            {
133 //                Layer l = layering.getLayer(i);
134 //                System.out.println("layering "+i);
135 //                List<LayerSlice> slices = l.getSlices();
136 //                for(int j=0; j<slices.size(); ++j)
137 //                {
138 //                    LayerSlice slice = slices.get(j);
139 //                    System.out.println("Layer "+i+" slice "+j+" is "+ slice.getMaterial().getName() +" and "+(slice.isSensitive() ? " is sensitive" : ""));
140 //                }
141 //            }
142 //            emCalInnerRadius = calsubBarrel.getInnerRadius();
143             //endcap
144 //            AbstractPolyhedraCalorimeter calsubEndcap = (AbstractPolyhedraCalorimeter) event.getDetector().getSubdetectors().get(det[1]);
145 //            emCalInnerZ = abs(calsubEndcap.getZMin());
146             if (skipFirstLayer) {
147                 System.out.println("processing " + event.getDetectorName() + " with an em calorimeter with a massless first gap");
148             }
149 //            System.out.println("Calorimeter bounds: r= " + emCalInnerRadius + " z= " + emCalInnerZ);
150             System.out.println("initialized...");
151 
152             _initialized = true;
153         }
154 
155         // organize the histogram tree by species and energy
156         List<MCParticle> mcparts = event.getMCParticles();
157         //MCParticle mcpart = mcparts.get(mcparts.size() - 1); // this only works if particle was generated by stdhep, does not work with GEANT gps
158         MCParticle mcpart = null;
159         // Look for the most energetic final state particle
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             // this analysis is intended for single particle calorimeter response.
186             // let's make sure that the primary particle did not interact in the
187             // tracker...
188 //        Hep3Vector endpoint = mcpart.getEndPoint();
189 //        // this is just crap. Why not use SpacePoint?
190 //        double radius = sqrt(endpoint.x()*endpoint.x()+endpoint.y()*endpoint.y());
191 //        double z = endpoint.z();
192 ////        System.out.println("Input MCParticle endpoint: r="+radius+" z= "+z);
193 //
194 //        boolean doit = true;
195 //        if(radius<emCalInnerRadius && abs(z) < emCalInnerZ) doit = false;
196 //        if(doit)
197             if (mcpart.getSimulatorStatus().isDecayedInCalorimeter()) {
198 //            // now let's check the em calorimeters...
199 //            // get all of the calorimeter hits...
200 //            List<CalorimeterHit> allHits = new ArrayList<CalorimeterHit>();
201 //            // and the list after cuts.
202 //            List<CalorimeterHit> hitsToCluster = new ArrayList<CalorimeterHit>();
203 //            int i = 0;
204 //            for(String name : collNames)
205 //            {
206 ////                System.out.println("fetching "+name+" from the event");
207 //                List<CalorimeterHit> hits = event.get(CalorimeterHit.class, name);
208 ////                System.out.println(name+ " has "+hits.size()+" hits");
209 //                // let's look at the hits and see if we need to cut on energy or time...
210 //                for(CalorimeterHit hit: hits)
211 //                {
212 //                    aida.histogram1D(name+" raw calorimeter cell energy",100, 0., mipHistMaxBinVal[i]).fill(hit.getRawEnergy());
213 //                    aida.histogram1D(name+" raw calorimeter cell energy full range",100, 0., 0.2).fill(hit.getRawEnergy());
214 ////                    aida.cloud1D(name+" raw calorimeter cell energies").fill(hit.getRawEnergy());
215 //                    aida.histogram1D(name+" calorimeter cell time",100,0., 200.).fill(hit.getTime());
216 //                    if(hit.getTime()<timeCut)
217 //                    {
218 //                        if(hit.getRawEnergy()>mipCut[i])
219 //                        {
220 //                            hitsToCluster.add(hit);
221 //                        }
222 //                    }
223 //                }
224 //                allHits.addAll(hits);
225 //                ++i;
226 //            }
227 //            System.out.println("ready to cluster "+hitsToCluster.size()+ " hits");
228                 String processedHitsName = _cond.getString("ProcessedHitsCollectionName");
229                 List<CalorimeterHit> hitsToCluster = _collectionmanager.getList(processedHitsName);//event.get(CalorimeterHit.class, processedHitsName);
230 
231                 if (_debug) {
232                     System.out.println("clustering " + hitsToCluster.size() + " hits");
233                 }
234                 // quick check
235 //            for(CalorimeterHit hit : hitsToCluster)
236 //            {
237 //                System.out.println("hit ");
238 //                System.out.println(hit.getLCMetaData().getName());
239 //            }
240                 // cluster the hits
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 //                System.out.println(c);
248                     aida.cloud1D("cluster energy for all clusters").fill(c.getEnergy());
249                 }
250 
251                 // proceed only if we found a single cluster above threshold
252                 // too restrictive! simply take the highest energy cluster
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                     // subtract the mass to get kinetic energy...
263                     double expectedEnergy = mcEnergy;
264                     // In case of protons and neutrons we expect only the kinetic energy in the shower
265                     if (mcpart.getPDGID() == 2212 || mcpart.getPDGID() == 2112) {
266                         expectedEnergy = mcEnergy - mcMass;
267                     }
268 //                System.out.println(mcpart.getType().getName()+" "+expectedEnergy);
269                     aida.cloud1D("measured - predicted energy").fill(clusterEnergy - expectedEnergy);
270 
271                     // let's now break down the cluster by component.
272                     // this analysis uses:
273                     // 1.) first 20 EM layers
274                     // 2.) next 10 EM layers
275                     // 3.) Had layers
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                         // FIXME Hard-coded name.
289                         if (detectorName.toUpperCase().startsWith("ECAL")) {
290 //                            if (layer >= firstEmStartLayer && layer < secondEmStartLayer) {
291 //                                type = 0;
292 //                            } else {
293 //                                type = 1;
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                         // FIXME Hard-coded name.
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 //                        if (_isHcalDigital == true) {
314 //                                nHcalHits += nHCalHits;
315 //                            }
316                     } // end of loop over hits in cluster
317                     if (_isHcalDigital == true) {
318                                 vals[3] = nHCalHits;
319                             }
320                     // set up linear least squares:
321                     // expectedEnergy = a*E1 + b*E2 +c*E3
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                 } // end of single cluster cut
333 
334 //            event.put("All Calorimeter Hits",allHits);
335 //            event.put("Hits To Cluster",hitsToCluster);
336                 event.put("Found Clusters", clusters);
337             }// end of check on decays outside tracker volume
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         // calculate the sampling fractions...
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 }