File: | plugins/monitoring/ST_online_tracking/JEventProcessor_ST_online_tracking.cc |
Location: | line 227, column 26 |
Description: | Value stored to 'phi_mom' is never read |
1 | // $Id$ |
2 | // |
3 | // File: JEventProcessor_ST_online_tracking.cc |
4 | // Created: Fri Jun 19 13:22:21 EDT 2015 |
5 | // Creator: mkamel (on Linux ifarm1401 2.6.32-431.el6.x86_64 x86_64) |
6 | // |
7 | // $Id$ |
8 | #include "JEventProcessor_ST_online_tracking.h" |
9 | using namespace jana; |
10 | using namespace std; |
11 | // Routine used to create our JEventProcessor |
12 | #include <JANA/JApplication.h> |
13 | #include <JANA/JFactory.h> |
14 | // ROOT header files |
15 | #include <TDirectory.h> |
16 | #include <TH1.h> |
17 | #include <TH2.h> |
18 | // ST header files |
19 | #include "START_COUNTER/DSCHit.h" |
20 | #include "START_COUNTER/DSCDigiHit.h" |
21 | // C++ header files |
22 | #include <stdint.h> |
23 | #include <vector> |
24 | #include <stdio.h> |
25 | // Include some JANA libraries |
26 | // PID libraries |
27 | #include "PID/DDetectorMatches.h" |
28 | #include "PID/DParticleID.h" |
29 | #include "PID/DChargedTrack.h" |
30 | // Tracking libraries |
31 | #include "TRACKING/DTrackTimeBased.h" |
32 | |
33 | // Define some constants |
34 | const double RAD2DEG = 57.29577951; // Convert radians to degrees |
35 | const uint32_t NCHANNELS = 30; // number of scintillator paddles |
36 | |
37 | // Declare 2D tracking histos |
38 | static TH2I *h2_r_vs_z; |
39 | static TH2I *h2_y_vs_x; |
40 | static TH2I *h2_phi_vs_sector; |
41 | static TH2I *h2_dphi_sector; |
42 | |
43 | static TH2I *h2_dedx_P_mag; |
44 | static TH2I *h2_dedx_P_mag_postv; |
45 | static TH2I *h2_dedx_P_mag_negtv; |
46 | |
47 | extern "C"{ |
48 | void InitPlugin(JApplication *app){ |
49 | InitJANAPlugin(app); |
50 | app->AddProcessor(new JEventProcessor_ST_online_tracking()); |
51 | } |
52 | } // "C" |
53 | |
54 | |
55 | //------------------ |
56 | // JEventProcessor_ST_online_tracking (Constructor) |
57 | //------------------ |
58 | JEventProcessor_ST_online_tracking::JEventProcessor_ST_online_tracking() |
59 | { |
60 | |
61 | } |
62 | |
63 | //------------------ |
64 | // ~JEventProcessor_ST_online_tracking (Destructor) |
65 | //------------------ |
66 | JEventProcessor_ST_online_tracking::~JEventProcessor_ST_online_tracking() |
67 | { |
68 | |
69 | } |
70 | |
71 | //------------------ |
72 | // init |
73 | //------------------ |
74 | jerror_t JEventProcessor_ST_online_tracking::init(void) |
75 | { |
76 | // This is called once at program startup. If you are creating |
77 | // and filling historgrams in this plugin, you should lock the |
78 | // ROOT mutex like this: |
79 | // |
80 | // japp->RootWriteLock(); |
81 | // ... fill historgrams or trees ... |
82 | // japp->RootUnLock(); |
83 | // |
84 | |
85 | // Create root folder for ST and cd to it, store main dir |
86 | TDirectory *main = gDirectory(TDirectory::CurrentDirectory()); |
87 | gDirectory(TDirectory::CurrentDirectory())->mkdir("st_tracking")->cd(); |
88 | // Book some 2D histos |
89 | h2_r_vs_z = new TH2I("r_vs_z", "Charged Track & ST Intersection; Z (cm); R (cm)", 280, 30, 100, 500, 0, 10); |
90 | h2_y_vs_x = new TH2I("y_vs_x", "Charged Track & ST Intersection; X (cm); Y (cm)", 2000, -20, 20, 1000, -10, 10); |
91 | |
92 | h2_phi_vs_sector = new TH2I("phi_vs_sector", "Charged Track & ST Intersection; Sector Number; #phi (deg)", NCHANNELS, 0.5, NCHANNELS + 0.5, 360, 0, 360); |
93 | h2_dphi_sector = new TH2I("h2_dphi_sector", "#delta #phi vs Sector; Sector Number; #delta #phi (deg)", NCHANNELS, 0.5, NCHANNELS + 0.5, 600, -15, 15); |
94 | h2_dedx_P_mag = new TH2I("h2_dedx_P_mag", "#frac{dE}{dx} vs Momentum; Momentum (Gev/c); dEdx (au)", 1000, 0, 2, 150, 0, .015); |
95 | h2_dedx_P_mag_postv= new TH2I("h2_dedx_P_mag_postv", "#frac{dE}{dx} vs Momentum (q > 0); Momentum (Gev/c); #frac{dE}{dx} (au)", 1000, 0, 2, 150, 0, .015); |
96 | h2_dedx_P_mag_negtv= new TH2I("h2_dedx_P_mag_negtv", "#frac{dE}{dx} vs Momentum (q < 0); Momentum (Gev/c); #frac{dE}{dx} (au)", 1000, 0, 2, 150, 0, .015); |
97 | |
98 | // cd back to main directory |
99 | gDirectory(TDirectory::CurrentDirectory())->cd("../"); |
100 | main->cd(); |
101 | |
102 | // for (uint32_t i = 0; i < NCHANNELS; i++) |
103 | // { |
104 | // phi_sec[i][0] = (6.0 + 12.0*double(i)) - 4.0; |
105 | // phi_sec[i][1] = (6.0 + 12.0*double(i)) + 4.0; |
106 | // } |
107 | return NOERROR; |
108 | } |
109 | //------------------ |
110 | // brun |
111 | //------------------ |
112 | jerror_t JEventProcessor_ST_online_tracking::brun(JEventLoop *eventLoop, int32_t runnumber) |
113 | { |
114 | // This is called whenever the run number changes |
115 | // Get the particleID object for each run |
116 | vector<const DParticleID *> dParticleID_algos; |
117 | eventLoop->Get(dParticleID_algos); |
118 | if(dParticleID_algos.size() < 1) |
119 | { |
120 | _DBG_std::cerr<<"plugins/monitoring/ST_online_tracking/JEventProcessor_ST_online_tracking.cc" <<":"<<120<<" "<<"Unable to get a DParticleID object! NO PID will be done!"<<endl; |
121 | return RESOURCE_UNAVAILABLE; |
122 | } |
123 | dParticleID = dParticleID_algos[0]; |
124 | return NOERROR; |
125 | } |
126 | |
127 | //------------------ |
128 | // evnt |
129 | //------------------ |
130 | jerror_t JEventProcessor_ST_online_tracking::evnt(JEventLoop *eventLoop, uint64_t eventnumber) |
131 | { |
132 | // This is called for every event. Use of common resources like writing |
133 | // to a file or filling a histogram should be mutex protected. Using |
134 | // loop->Get(...) to get reconstructed objects (and thereby activating the |
135 | // reconstruction algorithm) should be done outside of any mutex lock |
136 | // since multiple threads may call this method at the same time. |
137 | // Here's an example: |
138 | // |
139 | // vector<const MyDataClass*> mydataclasses; |
140 | // loop->Get(mydataclasses); |
141 | // |
142 | // japp->RootWriteLock(); |
143 | // ... fill historgrams or trees ... |
144 | // japp->RootUnLock(); |
145 | vector<const DSCDigiHit*> st_adc_digi_hits; |
146 | vector<const DParticleID*> pid_algorithm; |
147 | vector<const DSCHit*> st_hits; |
148 | eventLoop->Get(st_adc_digi_hits); |
149 | eventLoop->Get(pid_algorithm); |
150 | eventLoop->Get(st_hits); |
151 | vector<DVector3> sc_track_position; |
152 | vector<const DChargedTrack*> chargedTrackVector; |
153 | eventLoop->Get(chargedTrackVector); |
154 | // Grab the associated detector matches object |
155 | const DDetectorMatches* locDetectorMatches = NULL__null; |
156 | eventLoop->GetSingle(locDetectorMatches); |
157 | |
158 | // FILL HISTOGRAMS |
159 | // Since we are filling histograms local to this plugin, it will not interfere with other ROOT operations: can use plugin-wide ROOT fill lock |
160 | japp->RootFillLock(this); //ACQUIRE ROOT FILL LOCK |
161 | |
162 | sc_track_position.clear(); |
163 | |
164 | // Loop over charged tracks |
165 | for (uint32_t i = 0; i < chargedTrackVector.size(); i++) |
166 | { |
167 | // Grab the charged track |
168 | const DChargedTrack *thisChargedTrack = chargedTrackVector[i]; |
169 | |
170 | // Declare the time based track object |
171 | const DTrackTimeBased *timeBasedTrack; |
172 | |
173 | // Grab associated time based track object by selecting charged track with best FOM |
174 | thisChargedTrack->Get_BestTrackingFOM()->GetSingle(timeBasedTrack); |
175 | // Implement quality cuts for the time based tracks |
176 | float trackingFOMCut = 0.0027; // 3 sigma cut |
177 | if(timeBasedTrack->FOM < trackingFOMCut) continue; |
178 | // Get the charge of the track and cut on charged tracks |
179 | int q = timeBasedTrack->charge(); |
180 | // Grab the ST hit match params object and cut on only tracks matched to the ST |
181 | DSCHitMatchParams locSCHitMatchParams; |
182 | bool foundSC = dParticleID->Get_BestSCMatchParams(timeBasedTrack, locDetectorMatches, locSCHitMatchParams); |
183 | if (!foundSC) continue; |
184 | // Define vertex vector |
185 | DVector3 vertex; |
186 | // Vertex info |
187 | vertex = timeBasedTrack->position(); |
188 | // Cartesian Coordinates |
189 | double z_v = vertex.z(); |
190 | double r_v = vertex.Perp(); |
191 | // Liquid hydrogen target cuts |
192 | // Target length 30 cm, centered at z = 65.0 cm |
193 | // Upstream Diameter of target cell = 2.42 cm |
194 | // ID of scattering chamber = 7.49 cm |
195 | bool z_vertex_cut = 50.0 <= z_v && z_v <= 80.0; |
196 | bool r_vertex_cut = r_v < 0.5; |
197 | // applied vertex cut |
198 | if (!z_vertex_cut) continue; |
199 | if (!r_vertex_cut) continue; |
200 | // Declare a vector which quantizes the point of the intersection of a charged particle |
201 | // with a plane in the middle of the scintillator |
202 | DVector3 IntersectionPoint; |
203 | // Declare a vector which quantizes the unit vector of the charged particle track traversing |
204 | // through the scintillator with its origin at the intersection point |
205 | DVector3 IntersectionDir; |
206 | // Grab the paramteres associated to a track matched to the ST |
207 | vector<DSCHitMatchParams> st_params; |
208 | bool st_match = locDetectorMatches->Get_SCMatchParams(timeBasedTrack, st_params); |
209 | // If st_match = true, there is a match between this track and the ST |
210 | if (!st_match) continue; |
211 | |
212 | bool st_match_pid = dParticleID->MatchToSC(timeBasedTrack, |
213 | timeBasedTrack->rt, |
214 | st_params[0].dSCHit, |
215 | st_params[0].dSCHit->t, |
216 | locSCHitMatchParams, |
217 | true, |
218 | &IntersectionPoint, &IntersectionDir); |
219 | if(!st_match_pid) continue; |
220 | |
221 | DVector3 momentum_vec; |
222 | momentum_vec = timeBasedTrack->momentum(); |
223 | // applied vertex cut |
224 | DLorentzVector Lor_Mom = timeBasedTrack->lorentzMomentum(); |
225 | double P_mag = Lor_Mom.P(); |
226 | double phi_mom = momentum_vec.Phi()*RAD2DEG; |
227 | if (phi_mom < 0.0) phi_mom += 360.0; |
Value stored to 'phi_mom' is never read | |
228 | |
229 | if (st_match_pid) // Get the intersection point which can not be obtained from st_match |
230 | { |
231 | // Grab the sector |
232 | Int_t sector_m = st_params[0].dSCHit->sector; |
233 | //Acquire the energy loss per unit length in the ST (arbitrary units) |
234 | double dEdx = st_params[0].dEdx; |
235 | double dphi = st_params[0].dDeltaPhiToHit*RAD2DEG; |
236 | // Fill dEdx vs Momentum |
237 | h2_dedx_P_mag->Fill(P_mag,dEdx); |
238 | // Fill dEdx vs Momentum with cut on positive charges |
239 | if (q > 0) |
240 | { |
241 | h2_dedx_P_mag_postv->Fill(P_mag,dEdx); |
242 | } |
243 | // Fill dEdx vs Momentum with cut on negative charges |
244 | if (q < 0) |
245 | { |
246 | h2_dedx_P_mag_negtv->Fill(P_mag,dEdx); |
247 | } |
248 | // Obtain the intersection point with the ST |
249 | double phi_ip = IntersectionPoint.Phi()*RAD2DEG; |
250 | // Correct phi calculation |
251 | if (phi_ip < 0.0) phi_ip += 360.0; |
252 | // Acquire the intersection point |
253 | sc_track_position.push_back(IntersectionPoint); |
254 | //Fill 2D histos |
255 | h2_phi_vs_sector->Fill(sector_m,phi_ip); |
256 | h2_dphi_sector->Fill(sector_m,dphi); |
257 | } // PID match cut |
258 | } // Charged track loop |
259 | // Fill 2D histo |
260 | for (uint32_t i = 0; i < sc_track_position.size(); i++) |
261 | { |
262 | h2_r_vs_z->Fill(sc_track_position[i].z(),sc_track_position[i].Perp()); |
263 | h2_y_vs_x->Fill(sc_track_position[i].x(),sc_track_position[i].y()); |
264 | } |
265 | |
266 | japp->RootFillUnLock(this); //RELEASE ROOT FILL LOCK |
267 | |
268 | return NOERROR; |
269 | } |
270 | |
271 | //------------------ |
272 | // erun |
273 | //------------------ |
274 | jerror_t JEventProcessor_ST_online_tracking::erun(void) |
275 | { |
276 | // This is called whenever the run number changes, before it is |
277 | // changed to give you a chance to clean up before processing |
278 | // events from the next run number. |
279 | return NOERROR; |
280 | } |
281 | |
282 | //------------------ |
283 | // fini |
284 | //------------------ |
285 | jerror_t JEventProcessor_ST_online_tracking::fini(void) |
286 | { |
287 | // Called before program exit after event processing is finished. |
288 | return NOERROR; |
289 | } |
290 | |
291 | |
292 | |
293 | |
294 |