PandoraPFA/src/java/org/lcsim/recon/pfa/pandora
diff -u -r1.1 -r1.2
--- PandoraPFAProcessor.java 7 Nov 2007 18:41:30 -0000 1.1
+++ PandoraPFAProcessor.java 14 Nov 2007 04:18:30 -0000 1.2
@@ -0,0 +1,7904 @@
+package org.lcsim.recon.pfa.pandora;
+
+import java.util.ArrayList;
+import java.util.Iterator;
+import org.lcsim.recon.pfa.pandora.dummy.*;
+import org.lcsim.util.cpp.std;
+import java.util.List;
+import java.util.Map;
+import org.lcsim.recon.pfa.pandora.MyCaloHitExtended.CalorimeterHitType;
+import org.lcsim.recon.pfa.pandora.MyCaloHitExtended.CalorimeterHitZone;
+import org.lcsim.recon.pfa.pandora.MyTrackExtended.MyTrackStatus;
+import org.lcsim.recon.pfa.pandora.ProtoCluster.ProtoClusterClassificationType;
+import static java.lang.Math.sqrt;
+import static java.lang.Math.cos;
+import static java.lang.Math.sin;
+import static java.lang.Math.abs;
+import static java.lang.Math.atan2;
+import static java.lang.Math.acos;
+import static java.lang.Math.min;
+
+// Title : PandoraPFA cluster and particle flow object reconstruction
+// Author : M.Thomson
+// Version : v01-00
+// Last Update : 20/6/07
+// Recent Changes : Fix to use 1 or 2 byte Mokka Cell Codings
+// : Works with digital ECAL
+// : Improvements to reclustering
+// : Add PhotonRecovery and FragmentRemoval
+// : Added extra track killing (helix fit feature)
+// : Loosened 1 peak photon recovery interesting cut
+// : Bug fixes from Valgrind checks
+// : Clone shared photon track hits
+// : new photon recovery structure
+// : removed some obsolete and unhelpful code
+// 26/4/2007 : Performance
+// 45 GeV Jets : 29.5+-0.3 % 89.12
+// 100 GeV Jets : 30.7+-0.4 % 197.92
+// 180 GeV Jets : 41.9+-0.6 % 358.47
+// 250 GeV Jets : 53.0+-1.3 % 498.49
+//
+// 20/6/2007 : Add PerfectPFA
+// : get B-field from global
+// : allow possibility of energy-based
+// tagging a hit as being mip-like
+
+class PandoraPFAProcessor extends Processor {
+
+private class vec3 {
+ vec3(){
+ }
+ vec3(double x, double y, double z){
+ this.x = x;
+ this.y = y;
+ this.z = z;
+ }
+ double x;
+ double y;
+ double z;
+}
+
+private class ReclusterResult_t {
+ double chi2;
+ double chi;
+ int ntrackmatch;
+ int cluster0;
+ double chi0;
+ int cluster1;
+ double chi1;
+}
+
+ //*********** Algorithm Process parameters *************
+
+ private boolean _perfectPFA;
+ private int _maxLayer;
+
+ /** Vector of strings specifying ECAL
+ * collections of CalorimeterHits
+ **/
+ private List<String> _ecalCollections;
+
+ private String _relationCaloString;
+ private String _relationTrackString;
+
+ /** Vector of strings specifying HCAL
+ * collections of CalorimeterHits
+ **/
+ private List<String> _hcalCollections;
+
+ /** Vector of strings specifying Reconstructed Track
+ * collections
+ **/
+ private List<String> _trackCollections;
+
+ /** String specifying coding of cell ID
+ ** Mokka (default)
+ ** MokkaNew
+ ** Jupiter
+ ** SLIC
+ **/
+ private String _cellIDCoding;
+
+ /** String specifying coding of name of output cluster collection
+ **/
+ private String _clusterCollectionName;
+
+ /** String specifying output Particle
+ * collection
+ **/
+ private String _particleCollectionName;
+
+
+ /** ECAL and HCAL calibrations (energy to MIP)
+ * and threshold to apply to hits in MIPs
+ **/
+
+ private int _digitalHCAL;
+ private int _digitalECAL;
+ private double _ecalToMIP;
+ private double _hcalToMIP;
+ private double _ecalMIPThreshold;
+ private double _hcalMIPThreshold;
+ private double _ecalEMMIPToGeV;
+ private double _hcalEMMIPToGeV;
+ private double _ecalHadMIPToGeV;
+ private double _hcalHadMIPToGeV;
+
+ private int _maxHCALLayer;
+ /** algorithm configuration
+ **/
+
+ private int _isolationType;
+ private int _layersForIsolationECAL;
+ private double _distanceForIsolationECAL;
+ private int _layersForIsolationHCAL;
+ private double _distanceForIsolationHCAL;
+ private int _nSmallestClusterSizeForIsolation;
+ private double _isolationDistanceForReCombination;
+ private int _layersForDensityWeighting;
+ private int _densityWeightingPower;
+ private double _densityWeightCutECAL;
+ private double _densityWeightCutHCAL;
+ private int _minimumHitsInCluster;
+ private int _maximumHitsInSoftCluster;
+ private int _typeOfOrdering;
+ private int _layersToStepBackECAL;
+ private int _layersToStepBackHCAL;
+ private int _seedClustersWithTracks;
+ private int _extraLooperMatching;
+ private int _allowPhotonFragmentation;
+ private int _photonRecovery;
+ private double _photonFragPhotonDeltaChi2Cut;
+ private double _photonFragNonPhotonDeltaChi2Cut;
+ private int _allowReclustering;
+ private double _mipLikeMipCut;
+
+ private int _associateLoopingTracks;
+ private int _associateTrackSegments;
+ private int _associateShowersToMips;
+ private int _associateShowersToMipsII;
+ private int _associateBackScatters;
+ private int _associateMipStubsToShowers;
+ private int _associateFromStart;
+ private int _associateByProximity;
+ private int _associateByShowerCone;
+ private int _associateSoftClusters;
+ private int _fragmentRemoval;
+ private double _fragmentRemovalContactCut;
+ private int _fragmentRemovalContactLayersCut;
+ private double _fragmentRemovalMaxChi2ForMerge;
+ private double _fragmentRemovalDeltaChi2ForMerge;
+
+
+ private double _chiToAttemptReclustering;
+ private double _chiToForceReclustering;
+ private int _outerHcalRegionLayer;
+ private double _looperClosestApproachCutECAL;
+ private double _looperClosestApproachCutHCAL;
+ private double _looperTrackSeparationCut;
+ private int _clusterFormationStrategy;
+ private double _clusterFormationConeTanECAL;
+ private double _clusterFormationConeTanHCAL;
+ private double _clusterFormationPadsECAL;
+ private double _clusterFormationPadsHCAL;
+ private double _sameLayerPadCutECAL;
+ private double _sameLayerPadCutHCAL;
+ private int _trackingWeight;
+ private int _trackingCutOff;
+ private double _trackingTube;
+ private double _minimumClusterEnergyEM;
+ private double _minimumClusterEnergyHAD;
+ private double _softClusterEnergyHAD;
+ private double _coneAssociationCosine;
+ private double _maximumDistanceForConeAssociationECAL;
+ private double _maximumDistanceForConeAssociationHCAL;
+
+ private double _trackMergeCutEcal;
+ private double _trackMergeCutHcal;
+ private double _trackMergePerpCutEcal;
+ private double _trackMergePerpCutHcal;
+ private double _z0TrackCut;
+ private double _d0TrackCut;
+
+
+ private double _energyDeltaChi2ForCrudeMerging;
+ private double _chiToKillTrack;
+
+ //*********** Internal Workspace
+
+ private LCCollection _tracks;
+ private List<MyCaloHitExtended>[] _allCalHitsByPseudoLayer = new List[Pandora.MAX_NUMBER_OF_LAYERS];
+ private List<MyCaloHitExtended>[] _calHitsByPseudoLayer = new List[Pandora.MAX_NUMBER_OF_LAYERS];
+ private List<MyCaloHitExtended>[][] _calHitsBySectorAndPseudoLayer = new List[Pandora.NUMBER_OF_SECTORS][Pandora.MAX_NUMBER_OF_LAYERS];
+ private List<MyCaloHitExtended> _hitsForRemoval;
+ private List<MyCaloHitExtended>[] _isolatedCalHitsByPseudoLayer = new List[Pandora.MAX_NUMBER_OF_LAYERS];
+ private MyTrackExtendedVec _tracksAR;
+ private ProtoClusterVec _protoClusters;
+ private ProtoClusterVec _finalProtoClusters;
+
+
+ //*********** Internal Algorithm Parameters *************
+ private MCTree _mcTree;
+ private MyMCTree _myMcTree;
+ private LCCollection _relationCaloCollection;
+ private LCCollection _relationTrackCollection;
+ private LCRelationNavigator _caloNavigator;
+ private LCRelationNavigator _trackNavigator;
+
+ private int _nRun ;
+ private int _nEvt ;
+ private int _protoClustersMade;
+ private String _detectorName;
+ private String _rootFile;
+
+ /** Geometry used in reconstruction **/
+
+ // number of physical layers in calorimeter
+ private int _nEcalLayers;
+ private int _nHcalLayers;
+ // number of pseudolayers in calorimeter
+ private int _nPseudoLayers;
+ private int _nBarrelPseudoLayers;
+ private int _nEndcapPseudoLayers;
+
+ // Azimuthal structure
+ private int _symmetry;
+ private int _ECALsymmetry;
+ private int _HCALsymmetry;
+ private double _rOfBarrel;
+ private double _phiOfBarrel;
+ private double _zOfEndcap;
+ private double _zOfHcalEndcap;
+ private double _rOfEndcap;
+ private double _rInnerEcalEndcap;
+ private double _zOfBarrel;
+ private double _bField;
+ private double _tpcInnerR;
+ private double _tpcOuterR;
+ private double _tpcZmax;
+ private double _cosTPC;
+ private List<vec3> _barrelStaveDir;
+ private List<vec3> _barrelStaveDirHCAL;
+
+ // Layer Structure
+ private List<Integer> _ecalStructure;
+ private List<Double> _ecalPadSizes;
+ private List<Integer> _hcalStructure;
+ private List<Double> _hcalPadSizes;
+ private double[] _positionBarrelLayer = new double[Pandora.MAX_NUMBER_OF_LAYERS];
+ private double[] _positionEndcapLayer = new double[Pandora.MAX_NUMBER_OF_LAYERS];
+ private double[] _thicknessBarrelLayer = new double[Pandora.MAX_NUMBER_OF_LAYERS];
+ private double[] _thicknessEndcapLayer = new double[Pandora.MAX_NUMBER_OF_LAYERS];
+ private double[] _absThicknessBarrelLayer = new double[Pandora.MAX_NUMBER_OF_LAYERS];
+ private double[] _absThicknessEndcapLayer = new double[Pandora.MAX_NUMBER_OF_LAYERS];
+ private double[] _padSizeEcal = new double[Pandora.MAX_NUMBER_OF_LAYERS];
+ private double[] _padSizeHcal = new double[Pandora.MAX_NUMBER_OF_LAYERS];
+ private int _mipCells;
+ private int _mipMaxCellsHit;
+
+ private int _display;
+ private int _showProtoClusterFormation;
+ private int _showProtoClusterAssociation;
+ private int _printing;
+ private int _analyseClusters;
+
+ // maximum layer in event
+ private int _layerMax;
+ // HCAL lookup
+ private int[] _hcalLayerForIntRadius = new int[10000];
+
+ // some monitoring quantities/histograms
+ private double _etotal;
+ private double _etotalnu;
+ private double _etotalnufwd;
+
+ private static double pi = Math.PI;
+ private static double twopi = 2.0*pi;
+
+ void init() {
+
+ // Initialisation stage
+
+ // Print algorithm parameters
+ printParameters() ;
+
+ // zero counters/pointers
+ _nRun = 0 ;
+ _nEvt = 0 ;
+ //_recoDisplay = null;
+ _mcTree= null;
+ _myMcTree= null;
+
+ // if requested - start up internal root-based event display
+ //if(_recoDisplay==null && _display != 0)_recoDisplay = new MyRecoDisplay(1);
+
+}
+
+void processRunHeader( LCRunHeader run) {
+
+ _nRun++ ;
+ _detectorName = run.getDetectorName();
+
+ _bField = Global.parameters.getDoubleVal("BField");
+
+ // determine detector (still some hard-coded parameters)
+ std.cout(" DETECTOR : ",_detectorName,std.endl);
+
+ // TPC Geometry from GEAR
+ gear.TPCParameters pTPC = Global.GEAR.getTPCParameters();
+ gear.PadRowLayout2D pTPCpads = pTPC.getPadLayout();
+ List<Double>tpcExtent = pTPCpads.getPlaneExtent();
+ _tpcInnerR = tpcExtent.get(0);
+ _tpcOuterR = tpcExtent.get(1);
+ _tpcZmax = pTPC.getMaxDriftLength();
+ _cosTPC = _tpcZmax/sqrt(_tpcZmax*_tpcZmax+_tpcInnerR*_tpcInnerR);
+
+ // Calorimeter geometry from GEAR
+
+ gear.CalorimeterParameters pEcalBarrel = Global.GEAR.getEcalBarrelParameters();
+ gear.CalorimeterParameters pEcalEndcap = Global.GEAR.getEcalEndcapParameters();
+ gear.CalorimeterParameters pHcalBarrel = Global.GEAR.getHcalBarrelParameters();
+ gear.CalorimeterParameters pHcalEndcap = Global.GEAR.getHcalEndcapParameters();
+ gear.LayerLayout ecalBarrelLayout = pEcalBarrel.getLayerLayout();
+ gear.LayerLayout ecalEndcapLayout = pEcalEndcap.getLayerLayout();
+ gear.LayerLayout hcalBarrelLayout = pHcalBarrel.getLayerLayout();
+ gear.LayerLayout hcalEndcapLayout = pHcalEndcap.getLayerLayout();
+
+ // determine pseudo-geometry of detector (determined by ECAL barrel)
+ // symmetry 0 =cylinder, 1=prorotype, >1 = polygon
+ _ECALsymmetry = pEcalBarrel.getSymmetryOrder();
+ _HCALsymmetry = pHcalBarrel.getSymmetryOrder();
+ _symmetry = _ECALsymmetry;
+ _rOfBarrel = (double)pEcalBarrel.getExtent()[0];
+ _zOfBarrel = (double)pEcalBarrel.getExtent()[3];
+
+ if(_printing>1){
+ std.cout(" GEAR: ECAL Barrel NLayers : ",ecalBarrelLayout.getNLayers(),std.endl);
+ std.cout(" GEAR: ECAL EndCap NLayers : ",ecalEndcapLayout.getNLayers(),std.endl);
+ std.cout(" GEAR: HCAL Barrel NLayers : ",hcalBarrelLayout.getNLayers(),std.endl);
+ std.cout(" GEAR: HCAL EndCap NLayers : ",hcalEndcapLayout.getNLayers(),std.endl);
+ std.cout("GEAR : ",(double)pEcalEndcap.getExtent()[0],std.endl);
+ std.cout("GEAR : ",(double)pEcalEndcap.getExtent()[1],std.endl);
+ std.cout("GEAR : ",(double)pEcalEndcap.getExtent()[2],std.endl);
+ std.cout("GEAR : ",(double)pEcalEndcap.getExtent()[3],std.endl);
+ }
+
+ _zOfEndcap = (double)pEcalEndcap.getExtent()[2];
+ _rInnerEcalEndcap = (double)pEcalEndcap.getExtent()[0];
+ _zOfHcalEndcap = (double)pHcalEndcap.getExtent()[2];
+ _rOfEndcap = (double)pEcalEndcap.getExtent()[1];
+ _phiOfBarrel = (double)pEcalBarrel.getPhi0();
+ _nEcalLayers = ecalBarrelLayout.getNLayers();
+ _nHcalLayers = hcalBarrelLayout.getNLayers();
+
+ // Determine ECAL polygon angles
+ // Store radial vectors perpendicular to stave layers in _barrelStaveDir
+ if(_symmetry>1){
+ double nFoldSymmetry = _symmetry;
+ double phi0 = pEcalBarrel.getPhi0();
+ for(int i=0;i<_symmetry;++i){
+ double phi = phi0 + i*twopi/nFoldSymmetry;
+ vec3 staveDir = new vec3(cos(phi),sin(phi),0.);
+ _barrelStaveDir.add(staveDir);
+ }
+ }
+
+ // determine HCAL polygon angles
+ // Store radial vectors perpendicular to stave layers in _barrelStaveDirHCAL
+ if(_HCALsymmetry>1){
+ double nFoldSymmetry = _HCALsymmetry;
+ double phi0 = pHcalBarrel.getPhi0();
+ for(int i=0;i<_HCALsymmetry;++i){
+ double phi = phi0 + i*twopi/nFoldSymmetry;
+ vec3 staveDir = new vec3(cos(phi),sin(phi),0.);
+ _barrelStaveDirHCAL.add(staveDir);
+ }
+ }
+
+
+ // Using GEAR information define PSEUDOLAYERS in BARREL
+ // Internally PandoraPFA defines pseudolayers starting from 1
+ // with layer zero is reserved for track projections
+ int layer=1;
+ for(int i=0;i<ecalBarrelLayout.getNLayers();++i){
+ _positionBarrelLayer[layer] = ecalBarrelLayout.getDistance(i);
+ _thicknessBarrelLayer[layer] = ecalBarrelLayout.getThickness(i);
+ _absThicknessBarrelLayer[layer] = ecalBarrelLayout.getAbsorberThickness(i);
+ _padSizeEcal[i] = ecalBarrelLayout.getCellSize0(i);
+ // now position the active layer
+ double delta =(_absThicknessBarrelLayer[layer]+_thicknessBarrelLayer[layer])/2.0;
+ _positionBarrelLayer[layer]+= delta;
+ layer++;
+ }
+ for(int i=0;i<hcalBarrelLayout.getNLayers();++i){
+ _positionBarrelLayer[layer] = hcalBarrelLayout.getDistance(i);
+ _thicknessBarrelLayer[layer] = hcalBarrelLayout.getThickness(i);
+ _absThicknessBarrelLayer[layer] = hcalBarrelLayout.getAbsorberThickness(i);
+ _padSizeHcal[i] = hcalBarrelLayout.getCellSize0(i);
+ double delta =(_absThicknessBarrelLayer[layer]+_thicknessBarrelLayer[layer])/2.0;
+ _positionBarrelLayer[layer]+= delta;
+ layer++;
+ }
+ // add some non-physical layers in HCAL (pseudo-layer > max layer)
+ // these are only necessary if the number/thickness of barrel and endcap layers are different
+ for(int i=0;i<10;++i){
+ _positionBarrelLayer[layer+i] = _positionBarrelLayer[layer-1]+_thicknessBarrelLayer[layer-1];
+ _thicknessBarrelLayer[layer+i] = _thicknessBarrelLayer[layer-1];
+ _absThicknessBarrelLayer[layer+1] = _absThicknessBarrelLayer[layer-1];
+ _padSizeHcal[layer+i] = _padSizeHcal[layer-1];
+ }
+ // define layer zero for track projections
+ _positionBarrelLayer[0] = _positionBarrelLayer[1]-_thicknessBarrelLayer[1];
+ _thicknessBarrelLayer[0] = _thicknessBarrelLayer[1];
+ _absThicknessBarrelLayer[0] = _absThicknessBarrelLayer[1];
+ _padSizeEcal[0] = _padSizeEcal[1];
+
+ int nBarrelLayers = layer;
+
+
+ // Using GEAR information define PSEUDOLAYERS in ENDCAP
+ // Internally PandoraPFA defines pseudolayers starting from 1
+ // with layer zero is reserved for track projections
+ layer=1;
+ for(int i=0;i<ecalEndcapLayout.getNLayers();++i){
+ _positionEndcapLayer[layer] = ecalEndcapLayout.getDistance(i);
+
+ _thicknessEndcapLayer[layer] = ecalEndcapLayout.getThickness(i);
+ _absThicknessEndcapLayer[layer] = ecalEndcapLayout.getAbsorberThickness(i);
+ double delta =(_absThicknessEndcapLayer[layer]+_thicknessEndcapLayer[layer])/2.0;
+ _positionEndcapLayer[layer]+= delta;
+ layer++;
+ }
+ for(int i=0;i<hcalEndcapLayout.getNLayers();++i){
+ _positionEndcapLayer[layer] = hcalEndcapLayout.getDistance(i);
+ _thicknessEndcapLayer[layer] = hcalEndcapLayout.getThickness(i);
+ _absThicknessEndcapLayer[layer] = hcalEndcapLayout.getAbsorberThickness(i);
+ double delta =(_absThicknessEndcapLayer[layer]+_thicknessEndcapLayer[layer])/2.0;
+ _positionEndcapLayer[layer]+= delta;
+ layer++;
+ }
+ // add some non-physical layers in HCAL (pseudo-layer > max layer)
+ // these are only necessary if the number/thickness of barrel and endcap layers are different
+ for(int i=0;i<10;++i){
+ _positionEndcapLayer[layer+i] = _positionEndcapLayer[layer-1]+_thicknessEndcapLayer[layer-1];
+ _thicknessEndcapLayer[layer+i] = _thicknessEndcapLayer[layer-1];
+ _absThicknessEndcapLayer[layer+1] = _absThicknessEndcapLayer[layer-1];
+ }
+
+
+ // define layer zero : used to stote track projections
+ _positionEndcapLayer[0] = _positionEndcapLayer[1]-_thicknessEndcapLayer[1];
+ _thicknessEndcapLayer[0] = _thicknessEndcapLayer[1];
+ _absThicknessEndcapLayer[0] = _absThicknessEndcapLayer[1];
+ _padSizeEcal[0] = _padSizeEcal[1];
+
+ int nEndcapLayers = layer;
+
+ if(_printing>2){
+ for(int i=0;i<nBarrelLayers;++i){
+ std.cout("Barrel Layer ",i," : ",_positionBarrelLayer[i],std.endl);
+ }
+ for(int i=0;i<nBarrelLayers;++i){
+ std.cout("EndCap Layer ",i," : ",_positionEndcapLayer[i],std.endl);
+ }
+ }
+
+
+ // Number of Pseudo Layers
+ _nPseudoLayers = nBarrelLayers;
+ if(nEndcapLayers>nBarrelLayers)_nPseudoLayers=nEndcapLayers;
+ _nBarrelPseudoLayers = nBarrelLayers;
+ _nEndcapPseudoLayers = nEndcapLayers;
+
+
+ // set up correspondance between radius and HCAL barrel layer
+ this.MakeHCALBarrelLookup();
+
+ // Ultimately the hits passed to the clustering are self-describing
+ // However, a small amount of detector information has to be passed to the clusters
+ // This is done via the factory which is used to create all clusters
+ (ProtoClusterFactory.Instance()).ZOfEndcap(_zOfEndcap);
+ (ProtoClusterFactory.Instance()).ZOfBarrel(_zOfBarrel);
+ // INCLUDING A REALLY UGLY HACK TO SET PHOTON ID CUTS FOR LDC DETECTOR
+ double cut = 0.94;
+ if(_detectorName=="LDC01Sc")cut = 0.90;
+ if(_detectorName=="LDC01a")cut = 0.90;
+ if(_detectorName=="LDC01b")cut = 0.90;
+ if(_detectorName=="LDC01c")cut = 0.90;
+ if(_detectorName=="LDC01d")cut = 0.90;
+ (ProtoClusterFactory.Instance()).SetPhotonIDCuts(cut);
+
+}
+
+void MakeHCALBarrelLookup() {
+
+ // Build lookup table for correspondance between perp distance
+ // form origin in HCAL to HCAL layer
+
+ for(int iradius=0;iradius<10000;iradius++){
+ double radius = (double)iradius;
+ double barrelLayerDist=999.;
+ int closestLayer=0;
+ for(int ilayer = 0; ilayer <= _nPseudoLayers; ++ilayer){
+ if(abs(radius-_positionBarrelLayer[ilayer])<barrelLayerDist){
+ barrelLayerDist = abs(radius-_positionBarrelLayer[ilayer]);
+ closestLayer = ilayer;
+ }
+ }
+ _hcalLayerForIntRadius[iradius] = closestLayer;
+ }
+}
+
+
+
+
+void processEvent(LCEvent evt ) {
+
+ // Steering of the cluster formation and PFO construction
+
+ _nEvt++ ;
+
+ // clear the old event
+ this.Clear();
+
+ // PandoraPFA does not yet work for a digital HCAL
+ // would not be a great deal of work, just hasn't been done yet
+ if(_digitalHCAL!=0){
+ std.cout(" Digital HCAL option not yet implemented ! ",std.endl);
+ std.cout(" to request this feature send email to [log in to unmask] ",std.endl);
+ return;
+ }
+
+
+ // Unpack the MC information for DEBUGGING purposes (ultimately this will be removed)
+ List<String> strVec = evt.getCollectionNames() ;
+ for(String name : strVec){
+ LCCollection col = evt.getCollection(name);
+ int nelem = col.getNumberOfElements();
+ // on the first event print out available collections
+ if(_nEvt==1)std.cout("PandoraPFA.processEvent Input Collections : ",name," elements : ",nelem,std.endl);
+ if (col.getTypeName() == LCIO.MCPARTICLE ) {
+ _mcTree = new MCTree(col);
+ _myMcTree = new MyMCTree(col);
+ if(_printing>0)_myMcTree.PerfectPFO(1);
+ if(_printing<1)_myMcTree.PerfectPFO(0);
+ }
+ }
+
+ // **************************************************************
+ // **************** PandoraPFA INITIALIZATION STAGE *************
+ // **************************************************************
+ // Read CALO hits and tracks and store in internal extended objects
+ // the extended hits include information about energy and pseudolayer
+ // the extended tracks include extrapolation information
+ this.initialiseEvent(evt);
+ // Apply isolation cuts the CALO hits
+ this.stripIsolatedHits();
+ // Order hits within each pseudolayer either by energy, local denisty,...
+ this.OrderHitsInEachLayer();
+
+ // Could draw hits etc...
+ //_recoDisplay.DrawByMIP(_calHitsByPseudoLayer);
+ //_recoDisplay.DrawByDensity(_calHitsByPseudoLayer);
+ //_recoDisplay.DrawByEnergy(_allCalHitsByPseudoLayer);
+ //_recoDisplay.DrawByEnergy(_calHitsByPseudoLayer);
+
+ // *****************************************************************
+ // **************** PandoraPFA CLUSTER FORMATION STAGE *************
+ // *****************************************************************
+ // First create the design to which clusters are going to be built
+ // i.e. the basic formation/growth rules
+ // Configure the ProtoCluster factory with the protoClusterDesign
+
+ protoClusterDesign_t design = new protoClusterDesign_t();
+ design.tanConeAngleECAL = _clusterFormationConeTanECAL;
+ design.tanConeAngleHCAL = _clusterFormationConeTanHCAL;
+ design.additionalPadsECAL = _clusterFormationPadsECAL;
+ design.additionalPadsHCAL = _clusterFormationPadsHCAL;
+ design.sameLayerPadCutECAL = _sameLayerPadCutECAL;
+ design.sameLayerPadCutHCAL = _sameLayerPadCutHCAL;
+ design.maximumDistanceForConeAssociationECAL=_maximumDistanceForConeAssociationECAL;
+ design.maximumDistanceForConeAssociationHCAL=_maximumDistanceForConeAssociationHCAL;
+ design.trackingWeight =_trackingWeight;
+ design.trackingCutOff =_trackingCutOff;
+ design.trackingTube =_trackingTube;
+ // Pass design to factory
+ (ProtoClusterFactory.Instance()).SetDesign(design);
+ // Form the protoclusters
+
+ if(_perfectPFA){
+ // Perform perfect Particle Flow
+ this.PerfectPFA(evt);
+ }else{
+ // Form the protoclusters
+ this.MakeProtoClusters(_tracksAR,_protoClusters);
+
+ // ********************************************************************
+ // **************** PandoraPFA CLUSTER ASSOCIATION STAGE **************
+ // ********************************************************************
+ // associate protoClusters together
+ // NOTE this is one of the main parts of the algorithm!
+ ProtoClusterVec nfinalProtoClusters = new ProtoClusterVec();
+ if(_protoClusters.size()>0)this.associateProtoClusters(_protoClusters,_tracksAR,nfinalProtoClusters,true);
+
+ // ********************************************************************
+ // ************************** PandoraPFA CLEANUP *********************
+ // ********************************************************************
+ // cleanup clusters
+ if(_fragmentRemoval>0){
+ this.fragmentRemoval(nfinalProtoClusters,_tracksAR,_finalProtoClusters);
+ }
+ }
+
+ // ********************************************************************
+ // **************** PandoraPFA LCIO Object Creation *******************
+ // ********************************************************************
+ // create clusters collection
+ this.CreateClusterCollection(evt);
+
+ // create PFO collection
+ this.CreatePFOCollection(evt);
+
+ // PFA analysis
+ if(_analyseClusters>0)this.AnalysePFAPerformance(evt);
+ //if(_recoDisplay!=null)_recoDisplay.DrawByProtoClusterEnergy(_finalProtoClusters,3,RECO_VIEW_XY);
+
+ //if(_recoDisplay!=null)_recoDisplay.DrawByProtoClusterEnergy(_finalProtoClusters,3,RECO_VIEW_XZ);
+
+}
+
+
+void Clear(){
+
+ // Clean up after reconstructing an event
+ // Delete all temporary objects
+ // including the protoclusters, the extended hits and tracks
+
+ if(_printing>3)std.cout("Clear",std.endl);
+
+ // delete navigators
+ //if(_caloNavigator!=null)delete _caloNavigator;
+ //if(_trackNavigator!=null)delete _trackNavigator;
+ _caloNavigator = null;
+ _trackNavigator = null;
+
+ // delete the extended hits
+ for(int i=0;i<Pandora.MAX_NUMBER_OF_LAYERS;++i){
+ //for(int j =0;j<_allCalHitsByPseudoLayer[i].size();++j)delete _allCalHitsByPseudoLayer[i][j];
+ _allCalHitsByPseudoLayer[i].clear();
+ _calHitsByPseudoLayer[i].clear();
+ _isolatedCalHitsByPseudoLayer[i].clear();
+ }
+
+ // delete the extended tracks
+ //for(int i =0;i<_tracksAR.size();++i)delete _tracksAR.get(i);
+ _tracksAR.clear();
+
+ // delete all created ProtoClusters
+ (ProtoClusterFactory.Instance()).Delete();
+ _protoClusters.clear();
+ _finalProtoClusters.clear();
+ // and tell the factory to zero its cluster counter
+ ProtoClusterFactory.Instance().ResetCounter();
+
+ // remove the MC information
+ //if(_mcTree!=null)delete _mcTree;
+ _mcTree = null;
+ //if(_myMcTree!=null)delete _myMcTree;
+ _myMcTree = null;
+
+ // clear the event display
+ //if(_recoDisplay!=null){
+ // _recoDisplay.Clear();
+ // _recoDisplay.CleanUp();
+ //}
+
+}
+
+
+void PerfectPFA(LCEvent evt){
+
+
+ if(_caloNavigator==null || _trackNavigator==null){
+ std.cout(" ERROR : cant perform PerfectPFA",std.endl);
+ if(_caloNavigator==null)std.cout(" : no calorimeter hit navigator",std.endl);
+ if(_trackNavigator==null)std.cout(" :no track navigator",std.endl);
+ return;
+ }
+
+ List<MyMCPFO> mcPFOs = (_myMcTree.getMCPFOs());
+ ProtoCluster[] protoClusters = new ProtoCluster[mcPFOs.size()];
+ int[] nTracks = new int[mcPFOs.size()];
+ MyTrackExtended[][] extendedTracks = new MyTrackExtended[mcPFOs.size()][10];
+
+ for(int i =0;i<mcPFOs.size();++i)protoClusters[i] = null;
+ for(int i =0;i<mcPFOs.size();++i){
+ nTracks[i]=0;
+ for(int j =0;j<10;++j)extendedTracks[i][j] = null;
+ }
+ // loop over hits
+ for(int ilayer=0;ilayer<_maxLayer;++ilayer){
+ // loop over hits in layer and try and associate with existing clusters
+ for(int ihit=0;ihit<_calHitsByPseudoLayer[ilayer].size();++ihit){
+ MyCaloHitExtended hiti = _calHitsByPseudoLayer[ilayer].get(ihit);
+ CalorimeterHit calhit = (CalorimeterHit)(hiti.getCalorimeterHit());
+ // use navigator to match hit to simhit
+ LCObjectVec objectVec = _caloNavigator.getRelatedToObjects(calhit);
+ if (objectVec.size() > 0) {
+ SimCalorimeterHit simhit =
+ (SimCalorimeterHit)(objectVec.get(0));
+ if (simhit.getNMCParticles() > 0) {
+ // choose first mc particle contributing to hit
+ MCParticle part = simhit.getParticleCont(0);
+ // use myMcTree to match hit to perfect particle flow object
+ int jmc = _myMcTree.GetMCPFOIndex(part);
+ if(jmc<mcPFOs.size()){
+ // add it to the appropriate protocluster
+ if(protoClusters[jmc]==null){
+ protoClusters[jmc] = (ProtoClusterFactory.Instance()).Manufacture(hiti);
+ }else{
+ protoClusters[jmc].AddHit(hiti);
+ }
+ }
+ }
+ }
+ }
+ }
+
+ for(int i =0;i<mcPFOs.size();++i){
+ if(protoClusters[i]!=null){
+ protoClusters[i].Update(_maxLayer);
+ protoClusters[i].ClassifyYourself();
+ _finalProtoClusters.add(protoClusters[i]);
+ }
+ }
+
+ // now deal with the tracks
+ for(int it=0;it<_tracksAR.size();++it){
+ Track track = _tracksAR.get(it).getTrack();
+ std.cout(track,std.endl);
+ LCObjectVec objectVec = _trackNavigator.getRelatedToObjects(track);
+ if (objectVec.size() > 0){
+ MCParticle part = (MCParticle)(objectVec.get(0));
+ int jmc = _myMcTree.GetMCPFOIndex(part);
+ if(jmc>=0&&jmc<mcPFOs.size()){
+ extendedTracks[jmc][nTracks[jmc]]=_tracksAR.get(it);
+ nTracks[jmc]++;
+ }
+ }
+ }
+
+ double totalEnergy = 0.;
+ for(int i =0; i<mcPFOs.size(); ++i){
+ if(nTracks[i]>0){
+ double trackEnergy = 0;
+ for(int it=0;it<nTracks[i];it++){
+ MyTrackExtended trackAR = extendedTracks[i][it];
+ trackEnergy += trackAR.getEnergy();
+ }
+ totalEnergy+=trackEnergy;
+ }
+ }
+
+ for(int i =0;i<mcPFOs.size();++i){
+ if(nTracks[i]==0){
+ if(protoClusters[i]!=null){
+ protoClusters[i].SetClusterClassification(ProtoClusterClassificationType.PC_UNKNOWN);
+ if(mcPFOs.get(i).isPhoton())protoClusters[i].SetClusterClassification(ProtoClusterClassificationType.PC_PRIMARY_PHOTON);
+ if(protoClusters[i].IsPhoton()){
+ totalEnergy+=protoClusters[i].EnergyEM();
+ }else{
+ }
+ }
+ }
+ }
+
+}
+
+
+
+void associateProtoClusters(ProtoClusterVec input, MyTrackExtendedVec tracks, ProtoClusterVec finalClusters, boolean firstPass){
+
+ // Have initial clustering
+ // Use topological rules to join clusters together
+ // In many ways this is the guts of the algorithm and hides a lot of complexity
+
+ //if(_recoDisplay!=null && _showProtoClusterAssociation>0)_recoDisplay.DrawByProtoCluster(input,3);
+
+ // First look for looping tracks (end<.end matching)
+ if(_printing>2)std.cout("Find and associated looping tracks",std.endl);
+ if(_associateLoopingTracks>0)this.associateLoopingTracks(input,0);
+
+ // Look for associated track segements (end<.start matching)
+ if(_printing>2)std.cout("Match clusters using track segments",std.endl);
+ if(_associateTrackSegments>0)this.associateTrackSegments(input);
+
+
+ // Look for associated track segements (mip end<. shower start matching)
+ if(_printing>2)std.cout("Associate showers to mips",std.endl);
+ if(_associateShowersToMips>0)this.associateShowersToMips(input);
+ if(_associateShowersToMipsII>0)this.associateShowersToMipsII(input);
+ //if(_recoDisplay!=null && _showProtoClusterAssociation>0)_recoDisplay.DrawByProtoCluster(input,3);
+
+ if(_printing>2)std.cout("Associate back scatters",std.endl);
+ if(_associateBackScatters>0)this.associateBackScatters(input);
+ //if(_recoDisplay!=null && _showProtoClusterAssociation>0)_recoDisplay.DrawByProtoCluster(input,3);
+
+ if(_printing>2)std.cout("Associate mip stubs",std.endl);
+ if(_associateMipStubsToShowers>0)this.associateMipStubsToShowers(input);
+ //if(_recoDisplay!=null && _showProtoClusterAssociation>0)_recoDisplay.DrawByProtoCluster(input,3);
+
+ if(_printing>2)std.cout("Associate from Start",std.endl);
+ if(_associateFromStart>0)this.associateFromStart(input);
+ //if(_recoDisplay!=null && _showProtoClusterAssociation>0)_recoDisplay.DrawByProtoCluster(input,3);
+
+
+
+ // now use the cruder methods
+ ProtoClusterVec _joinedProtoClusters = new ProtoClusterVec();
+
+ this.combineProtoClusters(input,_joinedProtoClusters);
+ if(_printing>2)std.cout("Associate clusters by proximity",std.endl);
+ //if(_recoDisplay!=null && _showProtoClusterAssociation>0)_recoDisplay.DrawByProtoCluster(_joinedProtoClusters,3);
+
+
+ //std.cout(" NOT RUNNING this.associateByProximity ",std.endl);
+ if(_associateByProximity>0)this.associateByProximity(_joinedProtoClusters);
+
+
+ //if(_recoDisplay!=null && _showProtoClusterAssociation>0)_recoDisplay.DrawByProtoCluster(_joinedProtoClusters,3);
+
+ ProtoClusterVec _newJoinedProtoClusters = new ProtoClusterVec();
+ this.combineProtoClusters(_joinedProtoClusters,_newJoinedProtoClusters);
+ if(_printing>2)std.cout("Associate by shower cone",std.endl);
+ if(_associateByShowerCone>0)this.associateByShowerCone(_newJoinedProtoClusters);
+ //if(_recoDisplay!=null && _showProtoClusterAssociation>0)_recoDisplay.DrawByProtoCluster(_newJoinedProtoClusters,3);
+
+ ProtoClusterVec _combinedProtoClusters = new ProtoClusterVec();
+
+ if(_printing>2)std.cout("Combined Clusters",std.endl);
+ this.combineProtoClusters(_newJoinedProtoClusters,_combinedProtoClusters);
+ // _recoDisplay.DrawByProtoCluster(_combinedProtoClusters,3);
+
+ if(_allowPhotonFragmentation==1)this.findPhotonsMergedWithMips(_combinedProtoClusters);
+
+ if(_printing>2)std.cout("Assoc soft",std.endl);
+ if(_associateSoftClusters>0)this.associateSoftClusters(_combinedProtoClusters);
+ this.associateIsolatedHits(_combinedProtoClusters);
+ //_recoDisplay.DrawByProtoCluster(_combinedProtoClusters,3);
+
+ ProtoClusterVec _nearFinalProtoClusters = new ProtoClusterVec();
+ this.makeFinalProtoClusters(_combinedProtoClusters,_nearFinalProtoClusters);
+ if(_associateLoopingTracks>1)this.associateLoopingTracks(_nearFinalProtoClusters,1);
+
+ // the following methods split/merge clusters using tracking information
+ // not used in reclustering
+ if(firstPass && _allowReclustering>0){
+ // use tracking information to split
+ if(_photonRecovery>=2)this.PhotonRecovery(_nearFinalProtoClusters,tracks,false );
+ this.trackDrivenSplitMergedClusters(_nearFinalProtoClusters);
+ this.trackDrivenAssociation(_nearFinalProtoClusters,0);
+ this.ResolveTwoTrackAssociations(_nearFinalProtoClusters,tracks,false);
+ // do it again -
+ this.trackDrivenAssociation(_nearFinalProtoClusters,1);
+ this.NewResolveSingleTrackAssociations(_nearFinalProtoClusters,tracks,false);
+ this.NewResolveTwoTrackAssociations(_nearFinalProtoClusters,tracks,false);
+ this.NewResolveThreeTrackAssociations(_nearFinalProtoClusters,tracks,false);
+ this.NewTrackDrivenAssociation(_nearFinalProtoClusters,tracks,false);
+ //this.NewMergedPhotonSearch(_nearFinalProtoClusters,tracks,false );
+ if(_photonRecovery==1||_photonRecovery>2)this.PhotonRecovery(_nearFinalProtoClusters,tracks,false );
+ }
+ this.combineProtoClusters(_nearFinalProtoClusters,finalClusters);
+
+ // Soft Photon ID
+ for(int icluster=0;icluster<finalClusters.size();++icluster){
+ finalClusters.get(icluster).SoftPhotonID();
+ }
+
+
+}
+
+
+void fragmentRemoval(ProtoClusterVec protoClusters, MyTrackExtendedVec extendedTracks, ProtoClusterVec outputProtoClusters){
+
+ // have finished clustering
+ // Look in detail at all hadron clusters and decide
+ // Merge dubious looking clusters
+
+ this.ClusterTrackAssociation(protoClusters,extendedTracks,false);
+
+ if(_printing>1){
+ std.cout(" FRAGMENT REMOVAL *****************" ,std.endl);
+ std.cout(" FRAGMENT REMOVAL *****************" ,std.endl);
+ std.cout(" FRAGMENT REMOVAL *****************" ,std.endl);
+ std.cout(" FRAGMENT REMOVAL *****************" ,std.endl);
+ std.cout(" FRAGMENT REMOVAL *****************" ,std.endl);
+ }
+
+ // The first criterion considered is "contact"
+ // Contact is defined as the number of layers in which a hadron cluster
+ // is in contact (has hits close to) another cluster
+
+ int nclusters = protoClusters.size();
+ for(int icluster=0;icluster<nclusters;++icluster){
+ ProtoCluster pCluster = protoClusters.get(icluster);
+ MyTrackExtendedVec tracks = pCluster.GetTrackVec();
+ double chi2best = 9999.;
+ ProtoCluster bestParent = null;
+ if(tracks.size()==0 && !pCluster.IsPrimaryPhoton()){
+ double clusterE = pCluster.EnergyHAD();
+ double contactEnergy = 0.;
+ double contactClusterEnergy = 0.;
+ double contactTrackEnergy = 0.;
+ ProtoClusterVec contacts = new ProtoClusterVec();
+ int totalContact = 0;
+ for(int jcluster=0;jcluster<nclusters;++jcluster){
+ MyTrackExtendedVec xtracks = protoClusters.get(jcluster).GetTrackVec();
+ if(icluster!=jcluster && xtracks.size()!=0){
+ double nearTrackE = 0.;
+ double nearClusterE = protoClusters.get(jcluster).EnergyHAD();
+ double fraction = protoClusters.get(jcluster).FractionInRadialCone(pCluster,0.90);
+ for(int it=0;it<xtracks.size();it++){
+ nearTrackE+=xtracks.get(it).getEnergy();
+ double fractiont = pCluster.FractionInTrackCone(xtracks.get(it),0.90);
+ if(fractiont>fraction)fraction=fractiont;
+ }
+ clusterContact_t contact = protoClusters.get(jcluster).ContactLayers(pCluster,_fragmentRemovalContactCut);
+ double sigE = 0.6*sqrt(nearTrackE);
+ double chio = (nearTrackE-nearClusterE)/sigE;
+ double chin = (nearTrackE-nearClusterE-clusterE)/sigE;
+ double chi2o = chio*chio;
+ double chi2n = chin*chin;
+
+ if(contact.contactLayers>5 || (contact.contactLayers>2 && chio>-1.)){
+ totalContact+=contact.contactLayers;
+ contactEnergy+= nearTrackE-nearClusterE;
+ contactTrackEnergy+=nearTrackE;
+ contactClusterEnergy+=nearClusterE;
+ contacts.add(protoClusters.get(jcluster));
+ std.cout(" Contact : ",contact.contactLayers," ",nearTrackE," ",nearClusterE,std.endl);
+ }
+
+ if(contact.contactLayers>_fragmentRemovalContactLayersCut){
+ if(abs(chi2n)<_fragmentRemovalMaxChi2ForMerge){
+ if(_printing>1 && pCluster.FirstLayer()>_nEcalLayers-20){
+ std.cout(" Dodgy cluster ",contact.contactLayers," : ",pCluster.EnergyHAD()," : ",chi2o," : ",chi2n,std.endl);
+ ProtoClusterVec pCs = new ProtoClusterVec();
+ pCs.add(protoClusters.get(jcluster));
+ pCs.add(protoClusters.get(icluster));
+ //if(_recoDisplay!=null)_recoDisplay.DrawSingleProtoCluster(pCs,RECO_VIEW_XY);
+
+
+ }
+
+ int layersFromEdge=999;
+ layersFromEdge = this.LayersFromEdge(protoClusters.get(jcluster));
+
+ if(_printing>1 && pCluster.FirstLayer()<=_nEcalLayers-20){
+ std.cout(" slightly dodgy cluster ",contact.contactLayers," : ",pCluster.EnergyHAD()," : ",chi2o," : ",chi2n,std.endl);
+
+ if(chio>_chiToAttemptReclustering && chi2n<4.0 && layersFromEdge >3){
+ if(_printing>1){
+ std.cout("*****WILL REMOVE CLUSTER ",std.endl);
+ std.cout("layers from edge : ",layersFromEdge,std.endl);
+ std.cout("PhotonID : ",pCluster.IsPrimaryPhoton(),std.endl);
+ }
+ if(chi2n<chi2best){
+ chi2best = chi2n;
+ bestParent = protoClusters.get(jcluster);
+ }
+ }
+ }
+ if(pCluster.FirstLayer()>_nEcalLayers-10){
+ if(chi2n-chi2o<_fragmentRemovalDeltaChi2ForMerge){
+ if(_printing>1)std.cout("WILL REMOVE CLUSTER ",std.endl);
+ if(chi2n<chi2best){
+ chi2best = chi2n;
[truncated at 1000 lines; 6908 more skipped]