1 package org.lcsim.mc.fast.reconstructedparticle;
2
3 import hep.physics.particle.Particle;
4 import hep.physics.particle.properties.ParticlePropertyManager;
5 import hep.physics.particle.properties.ParticlePropertyProvider;
6 import hep.physics.particle.properties.ParticleType;
7 import java.util.ArrayList;
8 import java.util.Map;
9 import java.util.HashMap;
10 import java.util.List;
11 import java.util.Random;
12 import org.lcsim.event.Cluster;
13 import org.lcsim.event.EventHeader;
14 import org.lcsim.event.ReconstructedParticle;
15 import org.lcsim.event.Track;
16 import org.lcsim.mc.fast.cluster.ronan.ReconHADCluster;
17 import org.lcsim.mc.fast.tracking.ReconTrack;
18 import org.lcsim.util.Driver;
19 import org.lcsim.conditions.ConditionsEvent;
20 import org.lcsim.conditions.ConditionsListener;
21 import org.lcsim.conditions.ConditionsSet;
22
23 import static java.lang.Math.abs;
24 import org.lcsim.mc.fast.cluster.ronan.ReconEMCluster;
25 import org.lcsim.util.aida.AIDA;
26
27
28
29
30
31 public class MCFastReconstructedParticleDriver extends Driver implements ConditionsListener {
32 private boolean refPoint000;
33 private ParticlePropertyProvider ppp;
34 private IDResolutionTables IDEff;
35 private AIDA aida = AIDA.defaultInstance();
36 private ParticleType eminus;
37 private ParticleType eplus;
38 private ParticleType klong;
39 private ParticleType muminus;
40 private ParticleType muplus;
41 private ParticleType neutron;
42 private ParticleType photon;
43 private ParticleType pizero;
44 private ParticleType piplus;
45 private ParticleType piminus;
46 private ParticleType pplus;
47 private ParticleType pminus;
48 private ParticleType kplus;
49 private ParticleType kminus;
50
51
52 public MCFastReconstructedParticleDriver() {
53 this(true);
54 }
55
56 public MCFastReconstructedParticleDriver(boolean refPoint000) {
57 this(ParticlePropertyManager.getParticlePropertyProvider(), refPoint000);
58 }
59
60 public MCFastReconstructedParticleDriver(ParticlePropertyProvider ppp) {
61 this(ppp, true);
62 }
63
64 public MCFastReconstructedParticleDriver(ParticlePropertyProvider ppp, boolean refPoint000) {
65 this.refPoint000 = refPoint000;
66 this.ppp = ppp;
67
68 eminus = ppp.get(11);
69 eplus = ppp.get(-11);
70 klong = ppp.get(130);
71 muminus = ppp.get(13);
72 muplus = ppp.get(-13);
73 neutron = ppp.get(2112);
74 photon = ppp.get(22);
75 pizero = ppp.get(111);
76 piplus = ppp.get(211);
77 piminus = ppp.get(-211);
78 pplus = ppp.get(2212);
79 pminus = ppp.get(-2212);
80 kplus = ppp.get(321);
81 kminus = ppp.get(-321);
82 }
83
84 protected void process(EventHeader event) {
85
86 boolean hist = getHistogramLevel() > 0;
87
88 if (IDEff == null) {
89 ConditionsSet idconditions = getConditionsManager().getConditions("IDEfficiency");
90 idconditions.addConditionsListener(this);
91 IDEff = new IDResolutionTables(idconditions);
92 }
93
94 Random rand = getRandom();
95
96 List<Track> tracks = event.get(Track.class, "Tracks");
97 List<Cluster> clusters = event.get(Cluster.class, "Clusters");
98
99
100 Map<Particle, Track> m_pt = new HashMap<Particle, Track>();
101 Map<Particle, Cluster> m_pc = new HashMap<Particle, Cluster>();
102 Map<Cluster, Track> m_ct = new HashMap<Cluster, Track>();
103 Map<Track, Cluster> m_tc = new HashMap<Track, Cluster>();
104
105 for (Track t : tracks)
106 m_pt.put(((ReconTrack) t).getMCParticle(), t);
107 for (Cluster c : clusters)
108 m_pc.put((c instanceof ReconEMCluster ? ((ReconEMCluster) c).getMCParticle() : ((ReconHADCluster) c).getMCParticle()), c);
109 for (Track t : tracks)
110 m_tc.put(t, m_pc.get(((ReconTrack) t).getMCParticle()));
111 for (Cluster c : clusters)
112 m_ct.put(c, m_pt.get(c instanceof ReconEMCluster ? ((ReconEMCluster) c).getMCParticle() : ((ReconHADCluster) c).getMCParticle()));
113
114 List<ReconstructedParticle> rpList = new ArrayList<ReconstructedParticle>();
115
116 for (Track t : tracks) {
117 ParticleType type = null;
118 if (t instanceof ReconTrack) {
119 ReconTrack rt = (ReconTrack) t;
120 Particle p = rt.getMCParticle();
121 int pdgid = p.getPDGID();
122
123
124 if ((abs(pdgid) == 11) && (rand.nextDouble() < IDEff.getElectronEff())) {
125 type = rt.getCharge() > 0 ? eplus : eminus;
126 } else if ((abs(pdgid) == 13) && (rand.nextDouble() < IDEff.getMuonEff())) {
127 type = rt.getCharge() > 0 ? muplus : muminus;
128 } else if ((abs(pdgid) == 2212) && (rand.nextDouble() < IDEff.getProtonEff())) {
129 type = rt.getCharge() > 0 ? pplus : pminus;
130 } else if ((abs(pdgid) == 321) && (rand.nextDouble() < IDEff.getKaonEff())) {
131 type = rt.getCharge() > 0 ? kplus : kminus;
132 } else {
133 type = rt.getCharge() > 0 ? piplus : piminus;
134 }
135
136 if ((p.getEnergy() > 10) && (hist)) {
137 if (Math.abs(t.getTrackParameter(4)) < 1) {
138 aida.histogram1D("track-particle", 150, -0.0003, 0.0003).fill((Math.sqrt(t.getPX() * t.getPX() + t.getPY() * t.getPY() + t.getPZ() * t.getPZ()) - Math.sqrt(p.getPX() * p.getPX() + p.getPY() * p.getPY() + p.getPZ() * p.getPZ())) / (p.getPX() * p.getPX() + p.getPY() * p.getPY() + p.getPZ() * p.getPZ()));
139 }
140 if ((Math.abs(t.getTrackParameter(4)) < 2.16) && (Math.abs(t.getTrackParameter(4)) > 1.96) && (p.getMomentum().magnitude() > 80) && (p.getMomentum().magnitude() < 120)) {
141 aida.histogram1D("track-particle-cut", 150, -0.01, 0.01).fill((Math.sqrt(t.getPX() * t.getPX() + t.getPY() * t.getPY() + t.getPZ() * t.getPZ()) - Math.sqrt(p.getPX() * p.getPX() + p.getPY() * p.getPY() + p.getPZ() * p.getPZ())) / (Math.sqrt(p.getPX() * p.getPX() + p.getPY() * p.getPY() + p.getPZ() * p.getPZ())));
142 }
143 }
144
145
146 MCFastReconstructedParticle rp = new MCFastReconstructedParticle(t, type, p, m_tc.get(t), IDEff.getWtChgTrkCal(), refPoint000);
147 rpList.add(rp);
148
149 if (hist) {
150 aida.histogram1D("recon-particle", 150, -5, 5).fill((rp.getEnergy() - p.getEnergy()) / (Math.sqrt(p.getEnergy())));
151 }
152 }
153 }
154
155
156 for (Cluster c : clusters) {
157
158 Particle p = null;
159 ParticleType type = null;
160
161 if (c instanceof ReconEMCluster) {
162 ReconEMCluster emc = (ReconEMCluster) c;
163 p = emc.getMCParticle();
164 if (m_ct.get(c) != null)
165 continue;
166 type = photon;
167 if (hist) {
168 aida.histogram1D("photonCLS-particle", 150, -3, 3).fill((emc.getEnergy() - emc.getMCParticle().getEnergy()) / (Math.sqrt(emc.getMCParticle().getEnergy())));
169 }
170 }
171
172 else if (c instanceof ReconHADCluster) {
173 ReconHADCluster emc = (ReconHADCluster) c;
174 p = emc.getMCParticle();
175 int pdgid = p.getPDGID();
176 if ((abs(pdgid) == 2112) && (rand.nextDouble() < IDEff.getNeutronEff())) {
177 type = neutron;
178 } else {
179 type = klong;
180 }
181 if (m_ct.get(c) != null || c.getEnergy() < type.getMass())
182 continue;
183 if (hist) {
184 aida.histogram1D("hadronCLS-particle", 150, -3, 3).fill((emc.getEnergy() - emc.getMCParticle().getEnergy()) / (Math.sqrt(emc.getMCParticle().getEnergy())));
185 }
186
187 }
188 MCFastReconstructedParticle rp = new MCFastReconstructedParticle(c, type, p);
189 rpList.add(rp);
190 if (hist) {
191 aida.histogram1D("recon-particle", 150, -10, 10).fill(rp.getEnergy() - p.getEnergy());
192 }
193
194 }
195
196 event.put("MCFastReconstructedParticles", rpList, ReconstructedParticle.class, 0);
197
198 }
199
200 public void conditionsChanged(ConditionsEvent event) {
201 ConditionsSet idconditions = getConditionsManager().getConditions("IDEfficiency");
202 IDEff = new IDResolutionTables(idconditions);
203 }
204
205 }