1 package org.lcsim.mc.fast.cluster.ronan;
2
3 import java.util.ArrayList;
4 import java.util.List;
5 import java.util.Random;
6
7 import org.lcsim.mc.fast.MCFast;
8 import org.lcsim.conditions.ConditionsEvent;
9 import org.lcsim.conditions.ConditionsListener;
10 import org.lcsim.conditions.ConditionsSet;
11 import org.lcsim.event.Cluster;
12 import org.lcsim.event.LCRelation;
13 import org.lcsim.event.EventHeader;
14 import org.lcsim.event.MCParticle;
15 import org.lcsim.util.Driver;
16 import org.lcsim.event.base.MyLCRelation;
17
18
19
20
21
22
23 public class MCFastRonan extends Driver implements ConditionsListener {
24 private final static int ElecID = 11;
25 private final static int NuEID = 12;
26 private final static int MuID = 13;
27 private final static int NuMuID = 14;
28 private final static int NuTauID = 16;
29 private final static int PhotonID = 22;
30 private final static int Neutralino1 = 1000022;
31 private final static int Neutralino2 = 1000023;
32 private final static int Neutralino3 = 1000025;
33 private final static int Neutralino4 = 1000035;
34 private boolean defaultMC = true;
35 private String fsname;
36
37 private ClusterResolutionTables clusterParm;
38
39 public void setFSList(String fslist) {
40 fsname = fslist;
41 defaultMC = false;
42 }
43
44 protected void process(EventHeader event) {
45 if (defaultMC) {
46 fsname = "MCParticle";
47 } else {
48 if (!event.hasCollection(MCParticle.class, fsname)) {
49 System.err.println("Collection " + fsname + " not found. Default Final State particles being used");
50 fsname = "MCParticle";
51 }
52 }
53 if (clusterParm == null) {
54 ConditionsSet conditions = getConditionsManager().getConditions("ClusterParameters");
55 conditions.addConditionsListener(this);
56 clusterParm = new ClusterResolutionTables(conditions);
57 }
58 List<Cluster> cl = new ArrayList<Cluster>();
59 List<LCRelation> lcrelationList = new ArrayList<LCRelation>();
60
61 boolean hist = getHistogramLevel() > 0;
62
63 List<MCParticle> particles = event.get(MCParticle.class, fsname);
64 for (MCParticle p : particles) {
65
66
67 if (defaultMC) {
68 if (p.getGeneratorStatus() != p.FINAL_STATE) {
69 continue;
70 }
71 }
72
73 int PDGID = p.getPDGID();
74 int absPDGID = Math.abs(PDGID);
75 double charge = p.getCharge();
76
77
78 boolean neutrino = absPDGID == NuEID || absPDGID == NuMuID || absPDGID == NuTauID || absPDGID == Neutralino1 || absPDGID == Neutralino2 || absPDGID == Neutralino3 || absPDGID == Neutralino4;
79 if (neutrino) {
80 continue;
81 }
82
83 double E = p.getEnergy();
84 if (Double.isNaN(E)) {
85 continue;
86 }
87
88 double pt2 = p.getMomentum().magnitudeSquared() - p.getPZ() * p.getPZ();
89 double pt = Math.sqrt(pt2);
90 double ptot = p.getMomentum().magnitude();
91 double cosTheta = p.getPZ() / p.getMomentum().magnitude();
92
93 Random rand = getRandom();
94
95
96 if (absPDGID == PhotonID || absPDGID == ElecID) {
97
98
99
100
101
102
103 if (E < clusterParm.getEMOnset()) {
104 continue;
105 }
106 if (Math.abs(cosTheta) > clusterParm.getPolarEMOuter()) {
107 continue;
108 }
109
110 cl.add(new ReconEMCluster(clusterParm, rand, p, hist));
111 lcrelationList.add(new MyLCRelation(cl.get(cl.size() - 1), (MCParticle) p));
112
113 }
114
115
116 else if (absPDGID != MuID) {
117
118
119
120
121
122
123
124 if (E < clusterParm.getHADOnset()) {
125 continue;
126 }
127 if (Math.abs(cosTheta) > clusterParm.getPolarHADOuter()) {
128 continue;
129 }
130
131 cl.add(new ReconHADCluster(clusterParm, rand, p, hist));
132 lcrelationList.add(new MyLCRelation(cl.get(cl.size() - 1), (MCParticle) p));
133 }
134 }
135 double neg_energy_total = 0.;
136 double pos_energy_weight_total = 0.;
137 for (Cluster rcl : cl) {
138 if (Math.abs(((ReconCluster) rcl).getMCParticle().getCharge()) > Double.MIN_VALUE)
139 continue;
140 neg_energy_total += ((ReconCluster) rcl).getNegEnergy();
141 pos_energy_weight_total += ((ReconCluster) rcl).getNegEnergy() < 0. ? 0. : Math.min(((ReconCluster) rcl).getSigma(), ((ReconCluster) rcl).getEnergy());
142 }
143 MCFast.log.info(" MCFast neg_energy_total= " + neg_energy_total + " pos_energy_weight_total= " + pos_energy_weight_total);
144 if (neg_energy_total < -Double.MIN_VALUE)
145 for (Cluster rcl : cl)
146 if (Math.abs(((ReconCluster) rcl).getMCParticle().getCharge()) < Double.MIN_VALUE && ((ReconCluster) rcl).getNegEnergy() >= 0.)
147 ((ReconCluster) rcl).adjustEnergy(neg_energy_total, pos_energy_weight_total);
148 event.put("Clusters", cl, Cluster.class, 0);
149 event.put("ClustersToMCP", lcrelationList, LCRelation.class, 0);
150 }
151
152 public void conditionsChanged(ConditionsEvent event) {
153 ConditionsSet conditions = getConditionsManager().getConditions("ClusterParameters");
154 clusterParm = new ClusterResolutionTables(conditions);
155 }
156 }