File: | programs/Analysis/hdview2/MyProcessor.cc |
Location: | line 1037, column 34 |
Description: | Dereference of null pointer |
1 | // Author: David Lawrence June 25, 2004 | |||
2 | // | |||
3 | // | |||
4 | // MyProcessor.cc | |||
5 | // | |||
6 | ||||
7 | #include <iostream> | |||
8 | #include <vector> | |||
9 | #include <string> | |||
10 | using namespace std; | |||
11 | ||||
12 | #include <TLorentzVector.h> | |||
13 | #include <TLorentzRotation.h> | |||
14 | #include <TEllipse.h> | |||
15 | #include <TArc.h> | |||
16 | #include <TBox.h> | |||
17 | #include <TLine.h> | |||
18 | #include <TText.h> | |||
19 | #include <TVector3.h> | |||
20 | #include <TColor.h> | |||
21 | ||||
22 | #include <particleType.h> | |||
23 | #include "hdview2.h" | |||
24 | #include "hdv_mainframe.h" | |||
25 | #include "MyProcessor.h" | |||
26 | #include "TRACKING/DTrackHit.h" | |||
27 | #include "TRACKING/DQuickFit.h" | |||
28 | #include "TRACKING/DMagneticFieldStepper.h" | |||
29 | #include "TRACKING/DTrackCandidate_factory.h" | |||
30 | #include "TRACKING/DMCTrackHit.h" | |||
31 | #include "TRACKING/DMCThrown.h" | |||
32 | #include "TRACKING/DTrackWireBased.h" | |||
33 | #include "TRACKING/DTrackTimeBased.h" | |||
34 | #include "PID/DChargedTrack.h" | |||
35 | #include "TRACKING/DReferenceTrajectory.h" | |||
36 | #include "JANA/JGeometry.h" | |||
37 | #include "TRACKING/DMCTrajectoryPoint.h" | |||
38 | #include "FCAL/DFCALHit.h" | |||
39 | #include "FDC/DFDCGeometry.h" | |||
40 | #include "CDC/DCDCTrackHit.h" | |||
41 | #include "FDC/DFDCPseudo.h" | |||
42 | #include "FDC/DFDCIntersection.h" | |||
43 | #include "HDGEOMETRY/DGeometry.h" | |||
44 | #include "FCAL/DFCALGeometry.h" | |||
45 | #include "FCAL/DFCALHit.h" | |||
46 | #include <PID/DParticleSet.h> | |||
47 | #include <PID/DPhysicsEvent.h> | |||
48 | #include "PID/DNeutralParticle.h" | |||
49 | #include "PID/DNeutralShower.h" | |||
50 | #include "PID/DTwoGammaFit.h" | |||
51 | #include "BCAL/DBCALHit.h" | |||
52 | #include "DVector2.h" | |||
53 | ||||
54 | extern hdv_mainframe *hdvmf; | |||
55 | ||||
56 | // These are declared in hdv_mainframe.cc, but as static so we need to do it here as well (yechh!) | |||
57 | static float FCAL_Zmin = 622.8; | |||
58 | static float FCAL_Rmin = 6.0; | |||
59 | static float FCAL_Rmax = 212.0/2.0; | |||
60 | static float BCAL_Rmin = 65.0; | |||
61 | static float BCAL_Zlen = 390.0; | |||
62 | static float BCAL_Zmin = 212.0 - BCAL_Zlen/2.0; | |||
63 | ||||
64 | ||||
65 | static vector<vector <DFDCWire *> >fdcwires; | |||
66 | ||||
67 | ||||
68 | bool DMCTrajectoryPoint_track_cmp(const DMCTrajectoryPoint *a,const DMCTrajectoryPoint *b){ | |||
69 | // sort by track number and then by particle type, then by z-coordinate | |||
70 | // (yes, I saw the same track number have different particle types!) | |||
71 | if(a->track != b->track)return a->track < b->track; | |||
72 | if(a->part != b->part )return a->part < b->part; | |||
73 | return a->z < b->z; | |||
74 | } | |||
75 | ||||
76 | ||||
77 | MyProcessor *gMYPROC=NULL__null; | |||
78 | ||||
79 | //------------------------------------------------------------------ | |||
80 | // MyProcessor | |||
81 | //------------------------------------------------------------------ | |||
82 | MyProcessor::MyProcessor() | |||
83 | { | |||
84 | Bfield = NULL__null; | |||
85 | loop = NULL__null; | |||
86 | ||||
87 | // Tell factory to keep around a few density histos | |||
88 | //gPARMS->SetParameter("TRKFIND:MAX_DEBUG_BUFFERS", 16); | |||
89 | ||||
90 | gMYPROC = this; | |||
91 | } | |||
92 | ||||
93 | //------------------------------------------------------------------ | |||
94 | // ~MyProcessor | |||
95 | //------------------------------------------------------------------ | |||
96 | MyProcessor::~MyProcessor() | |||
97 | { | |||
98 | ||||
99 | } | |||
100 | ||||
101 | //------------------------------------------------------------------ | |||
102 | // init | |||
103 | //------------------------------------------------------------------ | |||
104 | jerror_t MyProcessor::init(void) | |||
105 | { | |||
106 | // Make sure detectors have been drawn | |||
107 | //if(!drew_detectors)DrawDetectors(); | |||
108 | ||||
109 | ||||
110 | ||||
111 | vector<JEventLoop*> loops = app->GetJEventLoops(); | |||
112 | if(loops.size()>0){ | |||
113 | vector<string> facnames; | |||
114 | loops[0]->GetFactoryNames(facnames); | |||
115 | ||||
116 | hdvmf = new hdv_mainframe(gClient->GetRoot(), 1400, 700); | |||
117 | hdvmf->SetCandidateFactories(facnames); | |||
118 | hdvmf->SetWireBasedTrackFactories(facnames); | |||
119 | hdvmf->SetTimeBasedTrackFactories(facnames); | |||
120 | hdvmf->SetReconstructedFactories(facnames); | |||
121 | hdvmf->SetChargedTrackFactories(facnames); | |||
122 | fulllistmf = hdvmf->GetFullListFrame(); | |||
123 | debugermf = hdvmf->GetDebugerFrame(); | |||
124 | } | |||
125 | ||||
126 | return NOERROR; | |||
127 | } | |||
128 | ||||
129 | //------------------------------------------------------------------ | |||
130 | // brun | |||
131 | //------------------------------------------------------------------ | |||
132 | jerror_t MyProcessor::brun(JEventLoop *eventloop, int runnumber) | |||
133 | { | |||
134 | ||||
135 | // Read in Magnetic field map | |||
136 | DApplication* dapp = dynamic_cast<DApplication*>(eventloop->GetJApplication()); | |||
137 | Bfield = dapp->GetBfield(); | |||
138 | const DGeometry *dgeom = dapp->GetDGeometry(runnumber); | |||
139 | dgeom->GetFDCWires(fdcwires); | |||
140 | ||||
141 | RootGeom = dapp->GetRootGeom(); | |||
142 | geom = dapp->GetDGeometry(runnumber); | |||
143 | ||||
144 | MATERIAL_MAP_MODEL="DGeometry"; | |||
145 | gPARMS->SetDefaultParameter("TRKFIT:MATERIAL_MAP_MODEL", MATERIAL_MAP_MODEL); | |||
146 | ||||
147 | eventloop->GetCalib("PID/photon_track_matching", photon_track_matching); | |||
148 | DELTA_R_FCAL = photon_track_matching["DELTA_R_FCAL"]; | |||
149 | ||||
150 | return NOERROR; | |||
151 | } | |||
152 | ||||
153 | //------------------------------------------------------------------ | |||
154 | // evnt | |||
155 | //------------------------------------------------------------------ | |||
156 | jerror_t MyProcessor::evnt(JEventLoop *eventLoop, int eventnumber) | |||
157 | { | |||
158 | if(!eventLoop)return NOERROR; | |||
159 | loop = eventLoop; | |||
160 | last_jevent.FreeEvent(); | |||
161 | last_jevent = loop->GetJEvent(); | |||
162 | ||||
163 | string source = "<no source>"; | |||
164 | if(last_jevent.GetJEventSource())source = last_jevent.GetJEventSource()->GetSourceName(); | |||
165 | ||||
166 | cout<<"----------- New Event "<<eventnumber<<" -------------"<<endl; | |||
167 | hdvmf->SetEvent(eventnumber); | |||
168 | hdvmf->SetSource(source.c_str()); | |||
169 | hdvmf->DoMyRedraw(); | |||
170 | ||||
171 | return NOERROR; | |||
172 | } | |||
173 | ||||
174 | //------------------------------------------------------------------ | |||
175 | // FillGraphics | |||
176 | //------------------------------------------------------------------ | |||
177 | void MyProcessor::FillGraphics(void) | |||
178 | { | |||
179 | /// Create "graphics" objects for this event given the current GUI settings. | |||
180 | /// | |||
181 | /// This method will create DGraphicSet objects that represent tracks, hits, | |||
182 | /// and showers for the event. It creates objects for both hits and | |||
183 | /// reconstructed entities. The "graphics" objects created here are | |||
184 | /// really just collections of 3D space points with flags indicating | |||
185 | /// whether they should be drawn as markers or lines and with what | |||
186 | /// color and size. The actual graphics objects are created for the | |||
187 | /// various views of the detector in hdv_mainframe. | |||
188 | ||||
189 | graphics.clear(); | |||
190 | graphics_xyA.clear(); // The objects placed in these will be deleted by hdv_mainframe | |||
191 | graphics_xyB.clear(); // The objects placed in these will be deleted by hdv_mainframe | |||
192 | graphics_xz.clear(); // The objects placed in these will be deleted by hdv_mainframe | |||
193 | graphics_yz.clear(); // The objects placed in these will be deleted by hdv_mainframe | |||
194 | ||||
195 | if(!loop)return; | |||
| ||||
196 | ||||
197 | vector<const DTrackCandidate*> trCand; | |||
198 | loop->Get(trCand); | |||
199 | vector<const DTrackTimeBased*> trTB; | |||
200 | loop->Get(trTB); | |||
201 | vector<const DTrackWireBased*> trWB; | |||
202 | loop->Get(trWB); | |||
203 | hdv_debugerframe *p = hdvmf->GetDebugerFrame(); | |||
204 | p->SetNTrCand(trCand.size()); | |||
205 | p->SetNTrWireBased(trWB.size()); | |||
206 | p->SetNTrTimeBased(trTB.size()); | |||
207 | ||||
208 | // BCAL hits | |||
209 | if(hdvmf->GetCheckButton("bcal")){ | |||
210 | vector<const DBCALHit*> bcalhits; | |||
211 | loop->Get(bcalhits); | |||
212 | ||||
213 | for(unsigned int i=0; i<bcalhits.size(); i++){ | |||
214 | const DBCALHit *hit = bcalhits[i]; | |||
215 | TPolyLine *poly = hdvmf->GetBCALPolyLine(hit->module, hit->layer, hit->sector); | |||
216 | ||||
217 | if(!poly)continue; | |||
218 | ||||
219 | double a = hit->E/0.02; | |||
220 | double f = sqrt(a>1.0 ? 1.0:a<0.0 ? 0.0:a); | |||
221 | //double grey = 0.8; | |||
222 | //double s = 1.0 - f; | |||
223 | ||||
224 | //float r = s*grey; | |||
225 | //float g = s*grey; | |||
226 | //float b = f*(1.0-grey) + grey; | |||
227 | //float b = 0.; | |||
228 | //float g = 0.; | |||
229 | //float r = f*(1.0-grey) + grey; | |||
230 | ||||
231 | float r = 1.; | |||
232 | float g = 1.-f; | |||
233 | float b = 0.2; | |||
234 | if(f<=0.0){ | |||
235 | r = 1.; | |||
236 | g = 1.; | |||
237 | b = 0.9; | |||
238 | } | |||
239 | ||||
240 | poly->SetFillColor(TColor::GetColor(r,g,b)); | |||
241 | poly->SetLineColor(TColor::GetColor(r,g,b)); | |||
242 | poly->SetLineWidth(1); | |||
243 | poly->SetFillStyle(3001); | |||
244 | } | |||
245 | } | |||
246 | ||||
247 | // FCAL hits | |||
248 | if(hdvmf->GetCheckButton("fcal")){ | |||
249 | vector<const DFCALHit*> fcalhits; | |||
250 | loop->Get(fcalhits); | |||
251 | ||||
252 | for(unsigned int i=0; i<fcalhits.size(); i++){ | |||
253 | const DFCALHit *hit = fcalhits[i]; | |||
254 | TPolyLine *poly = hdvmf->GetFCALPolyLine(hit->x, hit->y); | |||
255 | if(!poly)continue; | |||
256 | ||||
257 | #if 0 | |||
258 | double a = hit->E/0.005; | |||
259 | double f = sqrt(a>1.0 ? 1.0:a<0.0 ? 0.0:a); | |||
260 | double grey = 0.8; | |||
261 | double s = 1.0 - f; | |||
262 | ||||
263 | float r = s*grey; | |||
264 | float g = s*grey; | |||
265 | float b = f*(1.0-grey) + grey; | |||
266 | #endif | |||
267 | double s = log10(hit->E/0.005)/log10(1.0/0.005); // s=1 for 1GeV energy deposit | |||
268 | float r = 1.; | |||
269 | float g = 1.-s; | |||
270 | float b = 0.2; | |||
271 | if(s<0.0){ | |||
272 | r = 1.; | |||
273 | g = 1.; | |||
274 | b = 0.9; | |||
275 | } | |||
276 | ||||
277 | poly->SetFillColor(TColor::GetColor(r,g,b)); | |||
278 | } | |||
279 | } | |||
280 | ||||
281 | ||||
282 | // CDC hits | |||
283 | if(hdvmf->GetCheckButton("cdc")){ | |||
284 | vector<const DCDCTrackHit*> cdctrackhits; | |||
285 | loop->Get(cdctrackhits); | |||
286 | ||||
287 | for(unsigned int i=0; i<cdctrackhits.size(); i++){ | |||
288 | const DCDCWire *wire = cdctrackhits[i]->wire; | |||
289 | ||||
290 | int color = (cdctrackhits[i]->tdrift>-20 && cdctrackhits[i]->tdrift<400) ? kCyan:kYellow; | |||
291 | DGraphicSet gset(color, kLine, 1.0); | |||
292 | DVector3 dpoint=wire->origin-(wire->L/2.0)*wire->udir; | |||
293 | TVector3 tpoint(dpoint.X(),dpoint.Y(),dpoint.Z()); | |||
294 | gset.points.push_back(tpoint); | |||
295 | dpoint=wire->origin+(wire->L/2.0)*wire->udir; | |||
296 | tpoint.SetXYZ(dpoint.X(),dpoint.Y(),dpoint.Z()); | |||
297 | gset.points.push_back(tpoint); | |||
298 | graphics.push_back(gset); | |||
299 | ||||
300 | // Rings for drift times. | |||
301 | // NOTE: These are not perfect since they still have TOF in them | |||
302 | if(hdvmf->GetCheckButton("cdcdrift") && wire->stereo==0.0){ | |||
303 | double x = wire->origin.X(); | |||
304 | double y = wire->origin.Y(); | |||
305 | double dist1 = cdctrackhits[i]->dist; | |||
306 | TEllipse *e = new TEllipse(x, y, dist1, dist1); | |||
307 | e->SetLineColor(38); | |||
308 | e->SetFillStyle(0); | |||
309 | graphics_xyA.push_back(e); | |||
310 | ||||
311 | double dist2 = dist1 - 4.0*55.0E-4; // what if our TOF was 4ns? | |||
312 | e = new TEllipse(x, y, dist2, dist2); | |||
313 | e->SetLineColor(38); | |||
314 | e->SetFillStyle(0); | |||
315 | graphics_xyA.push_back(e); | |||
316 | } | |||
317 | } | |||
318 | } | |||
319 | ||||
320 | // FDC wire | |||
321 | if(hdvmf->GetCheckButton("fdcwire")){ | |||
322 | vector<const DFDCHit*> fdchits; | |||
323 | loop->Get(fdchits); | |||
324 | ||||
325 | for(unsigned int i=0; i<fdchits.size(); i++){ | |||
326 | const DFDCHit *fdchit = fdchits[i]; | |||
327 | if(fdchit->type!=0)continue; | |||
328 | const DFDCWire *wire =fdcwires[fdchit->gLayer-1][fdchit->element-1]; | |||
329 | if(!wire){ | |||
330 | _DBG_std::cerr<<"MyProcessor.cc"<<":"<<330<< " "<<"Couldn't find wire for gLayer="<<fdchit->gLayer<<" and element="<<fdchit->element<<endl; | |||
331 | continue; | |||
332 | } | |||
333 | ||||
334 | // Wire | |||
335 | int color = (fdchit->t>-50 && fdchit->t<400) ? kCyan:kYellow; | |||
336 | DGraphicSet gset(color, kLine, 1.0); | |||
337 | DVector3 dpoint=wire->origin-(wire->L/2.0)*wire->udir; | |||
338 | TVector3 tpoint(dpoint.X(),dpoint.Y(),dpoint.Z()); | |||
339 | gset.points.push_back(tpoint); | |||
340 | dpoint=wire->origin+(wire->L/2.0)*wire->udir; | |||
341 | tpoint.SetXYZ(dpoint.X(),dpoint.Y(),dpoint.Z()); | |||
342 | gset.points.push_back(tpoint); | |||
343 | graphics.push_back(gset); | |||
344 | } | |||
345 | } | |||
346 | ||||
347 | // FDC intersection hits | |||
348 | if(hdvmf->GetCheckButton("fdcintersection")){ | |||
349 | vector<const DFDCIntersection*> fdcints; | |||
350 | loop->Get(fdcints); | |||
351 | DGraphicSet gsetp(46, kMarker, 0.5); | |||
352 | ||||
353 | for(unsigned int i=0; i<fdcints.size(); i++){ | |||
354 | TVector3 tpos(fdcints[i]->pos.X(),fdcints[i]->pos.Y(), | |||
355 | fdcints[i]->pos.Z()); | |||
356 | gsetp.points.push_back(tpos); | |||
357 | } | |||
358 | graphics.push_back(gsetp); | |||
359 | } | |||
360 | ||||
361 | // FDC psuedo hits | |||
362 | if(hdvmf->GetCheckButton("fdcpseudo")){ | |||
363 | vector<const DFDCPseudo*> fdcpseudos; | |||
364 | loop->Get(fdcpseudos); | |||
365 | DGraphicSet gsetp(38, kMarker, 0.5); | |||
366 | ||||
367 | for(unsigned int i=0; i<fdcpseudos.size(); i++){ | |||
368 | const DFDCWire *wire = fdcpseudos[i]->wire; | |||
369 | ||||
370 | // Pseudo point | |||
371 | TVector3 pos(fdcpseudos[i]->xy.X(), | |||
372 | fdcpseudos[i]->xy.Y(), wire->origin.Z()); | |||
373 | gsetp.points.push_back(pos); | |||
374 | } | |||
375 | graphics.push_back(gsetp); | |||
376 | } | |||
377 | ||||
378 | // DMCThrown | |||
379 | if(hdvmf->GetCheckButton("thrown")){ | |||
380 | vector<const DMCThrown*> mcthrown; | |||
381 | loop->Get(mcthrown); | |||
382 | for(unsigned int i=0; i<mcthrown.size(); i++){ | |||
383 | int color=14; | |||
384 | double size=1.5; | |||
385 | if(mcthrown[i]->charge()==0.0) color = kGreen; | |||
386 | if(mcthrown[i]->charge() >0.0) color = kBlue; | |||
387 | if(mcthrown[i]->charge() <0.0) color = kRed; | |||
388 | switch(mcthrown[i]->type){ | |||
389 | case Gamma: | |||
390 | case Positron: | |||
391 | case Electron: | |||
392 | size = 1.0; | |||
393 | break; | |||
394 | case Pi0: | |||
395 | case PiPlus: | |||
396 | case PiMinus: | |||
397 | size = 2.0; | |||
398 | break; | |||
399 | case Neutron: | |||
400 | case Proton: | |||
401 | case AntiProton: | |||
402 | size = 3.0; | |||
403 | break; | |||
404 | } | |||
405 | AddKinematicDataTrack(mcthrown[i], color, size); | |||
406 | } | |||
407 | } | |||
408 | ||||
409 | // CDC Truth points | |||
410 | if(hdvmf->GetCheckButton("cdctruth")){ | |||
411 | vector<const DMCTrackHit*> mctrackhits; | |||
412 | loop->Get(mctrackhits); | |||
413 | DGraphicSet gset(12, kMarker, 0.5); | |||
414 | for(unsigned int i=0; i<mctrackhits.size(); i++){ | |||
415 | const DMCTrackHit *hit = mctrackhits[i]; | |||
416 | if(hit->system != SYS_CDC)continue; | |||
417 | ||||
418 | TVector3 pos(hit->r*cos(hit->phi), hit->r*sin(hit->phi), hit->z); | |||
419 | gset.points.push_back(pos); | |||
420 | } | |||
421 | graphics.push_back(gset); | |||
422 | } | |||
423 | ||||
424 | // FDC Truth points | |||
425 | if(hdvmf->GetCheckButton("fdctruth")){ | |||
426 | vector<const DMCTrackHit*> mctrackhits; | |||
427 | loop->Get(mctrackhits); | |||
428 | DGraphicSet gset(12, kMarker, 0.5); | |||
429 | for(unsigned int i=0; i<mctrackhits.size(); i++){ | |||
430 | const DMCTrackHit *hit = mctrackhits[i]; | |||
431 | if(hit->system != SYS_FDC)continue; | |||
432 | ||||
433 | TVector3 pos(hit->r*cos(hit->phi), hit->r*sin(hit->phi), hit->z); | |||
434 | gset.points.push_back(pos); | |||
435 | } | |||
436 | graphics.push_back(gset); | |||
437 | } | |||
438 | ||||
439 | // Track Hits for Track Candidates and Candidate trajectory in Debuger Window | |||
440 | for(unsigned int n=0; n<trCand.size(); n++){ | |||
441 | if (n>9) | |||
442 | break; | |||
443 | char str1[128]; | |||
444 | sprintf(str1,"Candidate%d",n+1); | |||
445 | ||||
446 | if(hdvmf->GetCheckButton(str1)){ | |||
447 | ||||
448 | int color = n+1; | |||
449 | if (color > 4) | |||
450 | color++; | |||
451 | if (color > 6) | |||
452 | color++; | |||
453 | ||||
454 | AddKinematicDataTrack(trCand[n], color, 1.5); | |||
455 | ||||
456 | vector<const DCDCTrackHit*> cdctrackhits; | |||
457 | trCand[n]->Get(cdctrackhits); | |||
458 | for(unsigned int i=0; i<cdctrackhits.size(); i++){ | |||
459 | const DCDCWire *wire = cdctrackhits[i]->wire; | |||
460 | DGraphicSet gset(color, kLine, 1.0); | |||
461 | DVector3 dpoint=wire->origin-(wire->L/2.0)*wire->udir; | |||
462 | TVector3 tpoint(dpoint.X(),dpoint.Y(),dpoint.Z()); | |||
463 | gset.points.push_back(tpoint); | |||
464 | dpoint=wire->origin+(wire->L/2.0)*wire->udir; | |||
465 | tpoint.SetXYZ(dpoint.X(),dpoint.Y(),dpoint.Z()); | |||
466 | gset.points.push_back(tpoint); | |||
467 | graphics.push_back(gset); | |||
468 | ||||
469 | } // end loop of cdc hits of track candidate | |||
470 | vector<const DFDCPseudo*> fdcpseudos; | |||
471 | trCand[n]->Get(fdcpseudos); | |||
472 | DGraphicSet gsetp(color, kMarker, 0.5); | |||
473 | ||||
474 | for(unsigned int i=0; i<fdcpseudos.size(); i++){ | |||
475 | const DFDCWire *wire = fdcpseudos[i]->wire; | |||
476 | ||||
477 | // Pseudo point | |||
478 | TVector3 pos(fdcpseudos[i]->xy.X(), | |||
479 | fdcpseudos[i]->xy.Y(), wire->origin.Z()); | |||
480 | gsetp.points.push_back(pos); | |||
481 | } | |||
482 | graphics.push_back(gsetp); | |||
483 | ||||
484 | } | |||
485 | } | |||
486 | ||||
487 | // Wire Based Track Hits and trajectory for Debuger Window | |||
488 | for(unsigned int n=0; n<trWB.size(); n++){ | |||
489 | if (n>9) | |||
490 | break; | |||
491 | char str1[128]; | |||
492 | sprintf(str1,"WireBased%d",n+1); | |||
493 | ||||
494 | if(hdvmf->GetCheckButton(str1)){ | |||
495 | ||||
496 | int color = trWB[n]->candidateid; | |||
497 | if (color > 4) | |||
498 | color++; | |||
499 | if (color > 6) | |||
500 | color++; | |||
501 | ||||
502 | AddKinematicDataTrack(trWB[n], color, 1.5); | |||
503 | ||||
504 | vector<const DCDCTrackHit*> cdctrackhits; | |||
505 | trWB[n]->Get(cdctrackhits); | |||
506 | for(unsigned int i=0; i<cdctrackhits.size(); i++){ | |||
507 | const DCDCWire *wire = cdctrackhits[i]->wire; | |||
508 | DGraphicSet gset(color, kLine, 1.0); | |||
509 | DVector3 dpoint=wire->origin-(wire->L/2.0)*wire->udir; | |||
510 | TVector3 tpoint(dpoint.X(),dpoint.Y(),dpoint.Z()); | |||
511 | gset.points.push_back(tpoint); | |||
512 | dpoint=wire->origin+(wire->L/2.0)*wire->udir; | |||
513 | tpoint.SetXYZ(dpoint.X(),dpoint.Y(),dpoint.Z()); | |||
514 | gset.points.push_back(tpoint); | |||
515 | graphics.push_back(gset); | |||
516 | ||||
517 | } // end loop of cdc hits of track candidate | |||
518 | vector<const DFDCPseudo*> fdcpseudos; | |||
519 | trWB[n]->Get(fdcpseudos); | |||
520 | DGraphicSet gsetp(color, kMarker, 0.5); | |||
521 | ||||
522 | for(unsigned int i=0; i<fdcpseudos.size(); i++){ | |||
523 | const DFDCWire *wire = fdcpseudos[i]->wire; | |||
524 | ||||
525 | // Pseudo point | |||
526 | TVector3 pos(fdcpseudos[i]->xy.X(), | |||
527 | fdcpseudos[i]->xy.Y(), wire->origin.Z()); | |||
528 | gsetp.points.push_back(pos); | |||
529 | } | |||
530 | graphics.push_back(gsetp); | |||
531 | ||||
532 | } | |||
533 | } | |||
534 | ||||
535 | // Time Based Track Hits and trajectory for Debuger Window | |||
536 | for(unsigned int n=0; n<trTB.size(); n++){ | |||
537 | if (n>9) | |||
538 | break; | |||
539 | char str1[128]; | |||
540 | sprintf(str1,"TimeBased%d",n+1); | |||
541 | ||||
542 | if(hdvmf->GetCheckButton(str1)){ | |||
543 | ||||
544 | int color = trTB[n]->candidateid; | |||
545 | if (color > 4) | |||
546 | color++; | |||
547 | if (color > 6) | |||
548 | color++; | |||
549 | ||||
550 | AddKinematicDataTrack(trTB[n], color, 1.5); | |||
551 | ||||
552 | vector<const DCDCTrackHit*> cdctrackhits; | |||
553 | trTB[n]->Get(cdctrackhits); | |||
554 | for(unsigned int i=0; i<cdctrackhits.size(); i++){ | |||
555 | const DCDCWire *wire = cdctrackhits[i]->wire; | |||
556 | DGraphicSet gset(color, kLine, 1.0); | |||
557 | DVector3 dpoint=wire->origin-(wire->L/2.0)*wire->udir; | |||
558 | TVector3 tpoint(dpoint.X(),dpoint.Y(),dpoint.Z()); | |||
559 | gset.points.push_back(tpoint); | |||
560 | dpoint=wire->origin+(wire->L/2.0)*wire->udir; | |||
561 | tpoint.SetXYZ(dpoint.X(),dpoint.Y(),dpoint.Z()); | |||
562 | gset.points.push_back(tpoint); | |||
563 | graphics.push_back(gset); | |||
564 | ||||
565 | } // end loop of cdc hits of track candidate | |||
566 | vector<const DFDCPseudo*> fdcpseudos; | |||
567 | trTB[n]->Get(fdcpseudos); | |||
568 | DGraphicSet gsetp(color, kMarker, 0.5); | |||
569 | ||||
570 | for(unsigned int i=0; i<fdcpseudos.size(); i++){ | |||
571 | const DFDCWire *wire = fdcpseudos[i]->wire; | |||
572 | ||||
573 | // Pseudo point | |||
574 | TVector3 pos(fdcpseudos[i]->xy.X(), | |||
575 | fdcpseudos[i]->xy.Y(), wire->origin.Z()); | |||
576 | gsetp.points.push_back(pos); | |||
577 | } | |||
578 | graphics.push_back(gsetp); | |||
579 | ||||
580 | } | |||
581 | } | |||
582 | ||||
583 | ||||
584 | ||||
585 | ||||
586 | // TOF Truth points | |||
587 | if(hdvmf->GetCheckButton("toftruth")){ | |||
588 | vector<const DMCTrackHit*> mctrackhits; | |||
589 | loop->Get(mctrackhits); | |||
590 | DGraphicSet gset(kBlack, kMarker, 0.5); | |||
591 | for(unsigned int i=0; i<mctrackhits.size(); i++){ | |||
592 | const DMCTrackHit *hit = mctrackhits[i]; | |||
593 | if(hit->system != SYS_TOF)continue; | |||
594 | ||||
595 | TVector3 pos(hit->r*cos(hit->phi), hit->r*sin(hit->phi), hit->z); | |||
596 | gset.points.push_back(pos); | |||
597 | } | |||
598 | graphics.push_back(gset); | |||
599 | } | |||
600 | ||||
601 | // BCAL Truth points | |||
602 | if(hdvmf->GetCheckButton("bcaltruth")){ | |||
603 | vector<const DMCTrackHit*> mctrackhits; | |||
604 | loop->Get(mctrackhits); | |||
605 | DGraphicSet gset(kBlack, kMarker, 1.0); | |||
606 | for(unsigned int i=0; i<mctrackhits.size(); i++){ | |||
607 | const DMCTrackHit *hit = mctrackhits[i]; | |||
608 | if(hit->system != SYS_BCAL)continue; | |||
609 | ||||
610 | TVector3 pos(hit->r*cos(hit->phi), hit->r*sin(hit->phi), hit->z); | |||
611 | gset.points.push_back(pos); | |||
612 | ||||
613 | TMarker *m = new TMarker(pos.X(), pos.Y(), 2); | |||
614 | graphics_xyA.push_back(m); | |||
615 | } | |||
616 | //graphics.push_back(gset); | |||
617 | } | |||
618 | ||||
619 | // FCAL Truth points | |||
620 | if(hdvmf->GetCheckButton("fcaltruth")){ | |||
621 | vector<const DFCALGeometry*> fcalgeometries; | |||
622 | vector<const DFCALHit*> mcfcalhits; | |||
623 | loop->Get(fcalgeometries); | |||
624 | loop->Get(mcfcalhits); | |||
625 | if(fcalgeometries.size()>0){ | |||
626 | const DFCALGeometry *fgeom = fcalgeometries[0]; | |||
627 | ||||
628 | DGraphicSet gset(kBlack, kMarker, 0.25); | |||
629 | for(unsigned int i=0; i<mcfcalhits.size(); i++){ | |||
630 | const DFCALHit *hit = mcfcalhits[i]; | |||
631 | ||||
632 | DVector2 pos_face = fgeom->positionOnFace(hit->row, hit->column); | |||
633 | TVector3 pos(pos_face.X(), pos_face.Y(), FCAL_Zmin); | |||
634 | gset.points.push_back(pos); | |||
635 | ||||
636 | TMarker *m = new TMarker(pos.X(), pos.Y(), 2); | |||
637 | //m->SetColor(kGreen); | |||
638 | //m->SetLineWidth(1); | |||
639 | graphics_xyB.push_back(m); | |||
640 | ||||
641 | TMarker *m1 = new TMarker(pos.Z(), pos.X(), 2); | |||
642 | graphics_xz.push_back(m1); | |||
643 | TMarker *m2 = new TMarker(pos.Z(), pos.Y(), 2); | |||
644 | graphics_yz.push_back(m2); | |||
645 | } | |||
646 | //graphics.push_back(gset); | |||
647 | } | |||
648 | } | |||
649 | ||||
650 | // BCAL reconstructed photons | |||
651 | if(hdvmf->GetCheckButton("recon_photons_bcal")){ | |||
652 | vector<const DNeutralParticle*> neutrals; | |||
653 | loop->Get(neutrals); | |||
654 | ||||
655 | DGraphicSet gset(kYellow+2, kMarker, 1.25); | |||
656 | gset.marker_style=21; | |||
657 | for(unsigned int i=0; i<neutrals.size(); i++){ | |||
658 | vector<const DNeutralShower*> locNeutralShowers; | |||
659 | neutrals[i]->GetT(locNeutralShowers); | |||
660 | DetectorSystem_t locDetectorSystem = locNeutralShowers[0]->dDetectorSystem; | |||
661 | if(locDetectorSystem == SYS_BCAL){ | |||
662 | TVector3 pos( locNeutralShowers[0]->dSpacetimeVertex.X(), | |||
663 | locNeutralShowers[0]->dSpacetimeVertex.Y(), | |||
664 | locNeutralShowers[0]->dSpacetimeVertex.Z()); | |||
665 | gset.points.push_back(pos); | |||
666 | ||||
667 | double dist2 = 2.0 + 5.0*locNeutralShowers[0]->dEnergy; | |||
668 | TEllipse *e = new TEllipse(pos.X(), pos.Y(), dist2, dist2); | |||
669 | e->SetLineColor(kGreen); | |||
670 | e->SetFillStyle(0); | |||
671 | e->SetLineWidth(2); | |||
672 | graphics_xyA.push_back(e); | |||
673 | } | |||
674 | } | |||
675 | //graphics.push_back(gset); | |||
676 | } | |||
677 | ||||
678 | // FCAL reconstructed photons | |||
679 | if(hdvmf->GetCheckButton("recon_photons_fcal")){ | |||
680 | vector<const DNeutralParticle*> neutrals; | |||
681 | loop->Get(neutrals); | |||
682 | DGraphicSet gset(kOrange, kMarker, 1.25); | |||
683 | gset.marker_style=2; | |||
684 | for(unsigned int i=0; i<neutrals.size(); i++){ | |||
685 | vector<const DNeutralShower*> locNeutralShowers; | |||
686 | neutrals[i]->GetT(locNeutralShowers); | |||
687 | DetectorSystem_t locDetectorSystem = locNeutralShowers[0]->dDetectorSystem; | |||
688 | if(locDetectorSystem == SYS_FCAL){ | |||
689 | ||||
690 | TVector3 pos( locNeutralShowers[0]->dSpacetimeVertex.X(), | |||
691 | locNeutralShowers[0]->dSpacetimeVertex.Y(), | |||
692 | locNeutralShowers[0]->dSpacetimeVertex.Z()); | |||
693 | gset.points.push_back(pos); | |||
694 | ||||
695 | double dist2 = 2.0 + 10.0*locNeutralShowers[0]->dEnergy; | |||
696 | TEllipse *e = new TEllipse(pos.X(), pos.Y(), dist2, dist2); | |||
697 | e->SetLineColor(kGreen); | |||
698 | e->SetFillStyle(0); | |||
699 | e->SetLineWidth(2); | |||
700 | graphics_xyB.push_back(e); | |||
701 | ||||
702 | TEllipse *e1 = new TEllipse(pos.Z(), pos.X(), dist2, dist2); | |||
703 | e1->SetLineColor(kGreen); | |||
704 | e1->SetFillStyle(0); | |||
705 | e1->SetLineWidth(2); | |||
706 | graphics_xz.push_back(e1); | |||
707 | TEllipse *e2 = new TEllipse(pos.Z(), pos.Y(), dist2, dist2); | |||
708 | e2->SetLineColor(kGreen); | |||
709 | e2->SetFillStyle(0); | |||
710 | e2->SetLineWidth(2); | |||
711 | graphics_yz.push_back(e2); | |||
712 | ||||
713 | } | |||
714 | } | |||
715 | //graphics.push_back(gset); | |||
716 | } | |||
717 | ||||
718 | // Reconstructed photons matched with tracks | |||
719 | if(hdvmf->GetCheckButton("recon_photons_track_match")){ | |||
720 | vector<const DChargedTrack*> ctracks; | |||
721 | loop->Get(ctracks); | |||
722 | for(unsigned int i=0; i<ctracks.size(); i++){ | |||
723 | const DChargedTrack *locCTrack = ctracks[i]; | |||
724 | vector<const DNeutralShower*> locNeutralShowers; | |||
725 | locCTrack->GetT(locNeutralShowers); | |||
726 | ||||
727 | if (!locNeutralShowers.size()) continue; | |||
728 | ||||
729 | ||||
730 | // Decide if this hit BCAL of FCAL based on z of position on calorimeter | |||
731 | bool is_bcal = (locNeutralShowers[0]->dDetectorSystem == SYS_BCAL ); | |||
732 | ||||
733 | // Draw on all frames except FCAL frame | |||
734 | DGraphicSet gset(kRed, kMarker, 1.25); | |||
735 | gset.marker_style = is_bcal ? 22:3; | |||
736 | TVector3 tpos( locNeutralShowers[0]->dSpacetimeVertex.X(), | |||
737 | locNeutralShowers[0]->dSpacetimeVertex.Y(), | |||
738 | locNeutralShowers[0]->dSpacetimeVertex.Z()); | |||
739 | gset.points.push_back(tpos); | |||
740 | graphics.push_back(gset); | |||
741 | ||||
742 | // For BCAL hits, don't draw them on FCAL pane | |||
743 | if(is_bcal)continue; | |||
744 | ||||
745 | // Draw on FCAL pane | |||
746 | double dist2 = 2.0 + 2.0*locNeutralShowers[0]->dEnergy; | |||
747 | TEllipse *e = new TEllipse(tpos.X(), tpos.Y(), dist2, dist2); | |||
748 | e->SetLineColor(gset.color); | |||
749 | e->SetFillStyle(0); | |||
750 | e->SetLineWidth(1); | |||
751 | TMarker *m = new TMarker(tpos.X(), tpos.Y(), gset.marker_style); | |||
752 | m->SetMarkerColor(gset.color); | |||
753 | m->SetMarkerSize(1.75); | |||
754 | graphics_xyB.push_back(e); | |||
755 | graphics_xyB.push_back(m); | |||
756 | } | |||
757 | } | |||
758 | ||||
759 | // FCAL and BCAL thrown photon projections | |||
760 | if(hdvmf->GetCheckButton("thrown_photons_fcal") || hdvmf->GetCheckButton("thrown_photons_bcal")){ | |||
761 | vector<const DMCThrown*> throwns; | |||
762 | loop->Get(throwns); | |||
763 | DGraphicSet gset(kSpring, kMarker, 1.25); | |||
764 | for(unsigned int i=0; i<throwns.size(); i++){ | |||
765 | const DMCThrown *thrown = throwns[i]; | |||
766 | if(thrown->charge()!=0.0)continue; | |||
767 | ||||
768 | // This may seem a little funny, but it saves swimming the reference trajectory | |||
769 | // multiple times. The GetIntersectionWithCalorimeter() method will find the | |||
770 | // intersection point of the photon with whichever calorimeter it actually hits | |||
771 | // or if it doesn't hit either of them. Then, we decide to draw the marker based | |||
772 | // on whether the flag is set to draw the calorimeter it hit or not. | |||
773 | DVector3 pos; | |||
774 | DetectorSystem_t who; | |||
775 | GetIntersectionWithCalorimeter(thrown, pos, who); | |||
776 | ||||
777 | if(who!=SYS_FCAL && who!=SYS_BCAL)continue; | |||
778 | if(who==SYS_FCAL && !hdvmf->GetCheckButton("thrown_photons_fcal"))continue; | |||
779 | if(who==SYS_BCAL && !hdvmf->GetCheckButton("thrown_photons_bcal"))continue; | |||
780 | TVector3 tpos(pos.X(),pos.Y(),pos.Z()); | |||
781 | gset.points.push_back(tpos); | |||
782 | ||||
783 | // Only draw on FCAL pane if photon hits FCAL | |||
784 | if(who==SYS_BCAL)continue; | |||
785 | ||||
786 | double dist2 = 2.0 + 2.0*thrown->energy(); | |||
787 | TEllipse *e = new TEllipse(pos.X(), pos.Y(), dist2, dist2); | |||
788 | e->SetLineColor(kSpring); | |||
789 | e->SetFillStyle(0); | |||
790 | e->SetLineWidth(4); | |||
791 | graphics_xyB.push_back(e); | |||
792 | } | |||
793 | graphics.push_back(gset); | |||
794 | } | |||
795 | ||||
796 | // FCAL and BCAL thrown charged particle projections | |||
797 | if(hdvmf->GetCheckButton("thrown_charged_fcal") || hdvmf->GetCheckButton("thrown_charged_bcal")){ | |||
798 | vector<const DMCThrown*> throwns; | |||
799 | loop->Get(throwns); | |||
800 | ||||
801 | for(unsigned int i=0; i<throwns.size(); i++){ | |||
802 | const DMCThrown *thrown = throwns[i]; | |||
803 | if(thrown->charge()==0.0)continue; | |||
804 | ||||
805 | // This may seem a little funny, but it saves swimming the reference trajectory | |||
806 | // multiple times. The GetIntersectionWithCalorimeter() method will find the | |||
807 | // intersection point of the photon with whichever calorimeter it actually hits | |||
808 | // or if it doesn't hit either of them. Then, we decide to draw the marker based | |||
809 | // on whether the flag is set to draw the calorimeter it hit or not. | |||
810 | DVector3 pos; | |||
811 | DetectorSystem_t who; | |||
812 | GetIntersectionWithCalorimeter(thrown, pos, who); | |||
813 | ||||
814 | if(who!=SYS_FCAL && who!=SYS_BCAL)continue; | |||
815 | if(who==SYS_FCAL && !hdvmf->GetCheckButton("thrown_charged_fcal"))continue; | |||
816 | if(who==SYS_BCAL && !hdvmf->GetCheckButton("thrown_charged_bcal"))continue; | |||
817 | ||||
818 | DGraphicSet gset(thrown->charge()>0.0 ? kBlue:kRed, kMarker, 1.25); | |||
819 | TVector3 tpos(pos.X(),pos.Y(),pos.Z()); | |||
820 | gset.points.push_back(tpos); | |||
821 | graphics.push_back(gset); | |||
822 | ||||
823 | //double dist2 = 6.0 + 2.0*thrown->momentum().Mag(); | |||
824 | double dist2 = DELTA_R_FCAL; | |||
825 | TEllipse *e = new TEllipse(pos.X(), pos.Y(), dist2, dist2); | |||
826 | e->SetLineColor(thrown->charge()>0.0 ? kBlue:kRed); | |||
827 | e->SetFillStyle(0); | |||
828 | e->SetLineWidth(4); | |||
829 | graphics_xyB.push_back(e); | |||
830 | } | |||
831 | } | |||
832 | ||||
833 | // FCAL and BCAL reconstructed charged particle projections | |||
834 | if(hdvmf->GetCheckButton("recon_charged_bcal") || hdvmf->GetCheckButton("recon_charged_fcal")){ | |||
835 | ||||
836 | // Here we used the full time-based track list, even though it includes multiple | |||
837 | // hypotheses for each candidate. This is because currently, the photon/track | |||
838 | // matching code in PID/DPhoton_factory.cc uses the DTrackTimeBased objects and | |||
839 | // the current purpose of drawing these is to see matching of reconstructed | |||
840 | // charged tracks with calorimeter clusters. | |||
841 | vector<const DTrackTimeBased*> tracks; | |||
842 | loop->Get(tracks, hdvmf->GetFactoryTag("DTrackTimeBased")); | |||
843 | ||||
844 | for(unsigned int i=0; i<tracks.size(); i++){ | |||
845 | const DTrackTimeBased *track = tracks[i]; | |||
846 | ||||
847 | // See notes for above sections | |||
848 | DVector3 pos; | |||
849 | DetectorSystem_t who; | |||
850 | GetIntersectionWithCalorimeter(track, pos, who); | |||
851 | ||||
852 | if(who!=SYS_FCAL && who!=SYS_BCAL)continue; | |||
853 | if(who==SYS_FCAL && !hdvmf->GetCheckButton("recon_charged_fcal"))continue; | |||
854 | if(who==SYS_BCAL && !hdvmf->GetCheckButton("recon_charged_bcal"))continue; | |||
855 | ||||
856 | DGraphicSet gset(track->charge()>0.0 ? kBlue:kRed, kMarker, 1.25); | |||
857 | TVector3 tpos(pos.X(),pos.Y(),pos.Z()); | |||
858 | gset.points.push_back(tpos); | |||
859 | graphics.push_back(gset); | |||
860 | ||||
861 | if(who==SYS_BCAL)continue; // Don't draw tracks hitting BCAL on FCAL pane | |||
862 | ||||
863 | //double dist2 = 6.0 + 2.0*track->momentum().Mag(); | |||
864 | double dist2 = DELTA_R_FCAL; | |||
865 | TEllipse *e = new TEllipse(pos.X(), pos.Y(), dist2, dist2); | |||
866 | e->SetLineColor(track->charge()>0.0 ? kBlue:kRed); | |||
867 | e->SetFillStyle(0); | |||
868 | e->SetLineWidth(4); | |||
869 | graphics_xyB.push_back(e); | |||
870 | } | |||
871 | } | |||
872 | ||||
873 | // DMCTrajectoryPoints | |||
874 | if(hdvmf->GetCheckButton("trajectories")){ | |||
875 | vector<const DMCTrajectoryPoint*> mctrajectorypoints; | |||
876 | loop->Get(mctrajectorypoints); | |||
877 | //sort(mctrajectorypoints.begin(), mctrajectorypoints.end(), DMCTrajectoryPoint_track_cmp); | |||
878 | ||||
879 | poly_type drawtype = hdvmf->GetCheckButton("trajectories_lines") ? kLine:kMarker; | |||
880 | double drawsize = hdvmf->GetCheckButton("trajectories_lines") ? 1.0:0.3; | |||
881 | DGraphicSet gset(kBlack, drawtype, drawsize); | |||
882 | //gset.marker_style = 7; | |||
883 | TVector3 last_point; | |||
884 | int last_track=-1; | |||
885 | int last_part=-1; | |||
886 | for(unsigned int i=0; i<mctrajectorypoints.size(); i++){ | |||
887 | const DMCTrajectoryPoint *pt = mctrajectorypoints[i]; | |||
888 | ||||
889 | switch(pt->part){ | |||
890 | case Gamma: | |||
891 | if(!hdvmf->GetCheckButton("trajectories_photon"))continue; | |||
892 | break; | |||
893 | case Electron: | |||
894 | if(!hdvmf->GetCheckButton("trajectories_electron"))continue; | |||
895 | break; | |||
896 | case Positron: | |||
897 | if(!hdvmf->GetCheckButton("trajectories_positron"))continue; | |||
898 | break; | |||
899 | case Proton: | |||
900 | if(!hdvmf->GetCheckButton("trajectories_proton"))continue; | |||
901 | break; | |||
902 | case Neutron: | |||
903 | if(!hdvmf->GetCheckButton("trajectories_neutron"))continue; | |||
904 | break; | |||
905 | case PiPlus: | |||
906 | if(!hdvmf->GetCheckButton("trajectories_piplus"))continue; | |||
907 | break; | |||
908 | case PiMinus: | |||
909 | if(!hdvmf->GetCheckButton("trajectories_piminus"))continue; | |||
910 | break; | |||
911 | default: | |||
912 | if(!hdvmf->GetCheckButton("trajectories_other"))continue; | |||
913 | break; | |||
914 | } | |||
915 | ||||
916 | TVector3 v(pt->x, pt->y, pt->z); | |||
917 | ||||
918 | if(i>0){ | |||
919 | //if((v-last_point).Mag() > 10.0){ | |||
920 | if(pt->track!=last_track || pt->part!=last_part){ | |||
921 | if(hdvmf->GetCheckButton("trajectories_colors")){ | |||
922 | switch(last_part){ | |||
923 | case Gamma: | |||
924 | gset.color = kOrange; | |||
925 | break; | |||
926 | case Electron: | |||
927 | case PiMinus: | |||
928 | gset.color = kRed+2; | |||
929 | break; | |||
930 | case Positron: | |||
931 | case Proton: | |||
932 | case PiPlus: | |||
933 | gset.color = kBlue+1; | |||
934 | break; | |||
935 | case Neutron: | |||
936 | gset.color = kGreen+2; | |||
937 | break; | |||
938 | default: | |||
939 | gset.color = kBlack; | |||
940 | break; | |||
941 | } | |||
942 | }else{ | |||
943 | gset.color = kBlack; | |||
944 | } | |||
945 | graphics.push_back(gset); | |||
946 | gset.points.clear(); | |||
947 | } | |||
948 | } | |||
949 | ||||
950 | gset.points.push_back(v); | |||
951 | last_point = v; | |||
952 | last_track = pt->track; | |||
953 | last_part = pt->part; | |||
954 | } | |||
955 | ||||
956 | if(hdvmf->GetCheckButton("trajectories_colors")){ | |||
957 | switch(last_part){ | |||
958 | case Gamma: | |||
959 | gset.color = kOrange; | |||
960 | break; | |||
961 | case Electron: | |||
962 | case PiMinus: | |||
963 | gset.color = kRed+2; | |||
964 | break; | |||
965 | case Positron: | |||
966 | case Proton: | |||
967 | case PiPlus: | |||
968 | gset.color = kBlue+1; | |||
969 | break; | |||
970 | case Neutron: | |||
971 | gset.color = kGreen+2; | |||
972 | break; | |||
973 | default: | |||
974 | gset.color = kBlack; | |||
975 | break; | |||
976 | } | |||
977 | }else{ | |||
978 | gset.color = kBlack; | |||
979 | } | |||
980 | graphics.push_back(gset); | |||
981 | } | |||
982 | ||||
983 | // DTrackCandidate | |||
984 | if(hdvmf->GetCheckButton("candidates")){ | |||
985 | vector<const DTrackCandidate*> trackcandidates; | |||
986 | loop->Get(trackcandidates, hdvmf->GetFactoryTag("DTrackCandidate")); | |||
987 | for(unsigned int i=0; i<trackcandidates.size(); i++){ | |||
988 | int color=i+1; | |||
989 | double size=2.0; | |||
990 | //if(trackcandidates[i]->charge() >0.0) color += 100; // lighter shade | |||
991 | //if(trackcandidates[i]->charge() <0.0) color += 150; // darker shade | |||
992 | AddKinematicDataTrack(trackcandidates[i], color, size); | |||
993 | } | |||
994 | } | |||
995 | ||||
996 | // DTrackWireBased | |||
997 | if(hdvmf->GetCheckButton("wiretracks")){ | |||
998 | vector<const DTrackWireBased*> wiretracks; | |||
999 | loop->Get(wiretracks, hdvmf->GetFactoryTag("DTrackWireBased")); | |||
1000 | for(unsigned int i=0; i<wiretracks.size(); i++){ | |||
1001 | AddKinematicDataTrack(wiretracks[i], (wiretracks[i]->charge()>0.0 ? kBlue:kRed)+2, 1.25); | |||
1002 | } | |||
1003 | } | |||
1004 | ||||
1005 | // DTrackTimeBased | |||
1006 | if(hdvmf->GetCheckButton("timetracks")){ | |||
1007 | vector<const DTrackTimeBased*> timetracks; | |||
1008 | loop->Get(timetracks, hdvmf->GetFactoryTag("DTrackTimeBased")); | |||
1009 | for(unsigned int i=0; i<timetracks.size(); i++){ | |||
1010 | AddKinematicDataTrack(timetracks[i], (timetracks[i]->charge()>0.0 ? kBlue:kRed)+0, 1.00); | |||
1011 | } | |||
1012 | } | |||
1013 | ||||
1014 | // DChargedTrack | |||
1015 | if(hdvmf->GetCheckButton("chargedtracks")){ | |||
1016 | vector<const DChargedTrack*> chargedtracks; | |||
1017 | loop->Get(chargedtracks, hdvmf->GetFactoryTag("DChargedTrack")); | |||
1018 | for(unsigned int i=0; i<chargedtracks.size(); i++){ | |||
1019 | int color=kViolet-3; | |||
1020 | double size=1.25; | |||
1021 | if (chargedtracks[i]->Get_Charge() > 0) color=kMagenta; | |||
1022 | ||||
1023 | if (chargedtracks[i]->Get_BestFOM()->mass() > 0.9) size=2.5; | |||
1024 | AddKinematicDataTrack(chargedtracks[i]->Get_BestFOM(),color,size); | |||
1025 | } | |||
1026 | } | |||
1027 | ||||
1028 | // DNeutralParticles | |||
1029 | if(hdvmf->GetCheckButton("neutrals")){ | |||
1030 | vector<const DNeutralParticle*> photons; | |||
1031 | loop->Get(photons, hdvmf->GetFactoryTag("DNeutralParticle")); | |||
1032 | ||||
1033 | for(unsigned int i=0; i<photons.size(); i++){ | |||
1034 | int color = kBlack; | |||
1035 | vector<const DNeutralShower*> locNeutralShowers; | |||
1036 | photons[i]->GetT(locNeutralShowers); | |||
1037 | DetectorSystem_t locDetSys = locNeutralShowers[0]->dDetectorSystem; | |||
| ||||
1038 | if(locDetSys==SYS_FCAL)color = kOrange; | |||
1039 | if(locDetSys==SYS_BCAL)color = kYellow+2; | |||
1040 | //if(locDetSys==DPhoton::kCharge)color = kRed; | |||
1041 | AddKinematicDataTrack(photons[i]->Get_BestFOM(), color, 1.00); | |||
1042 | } | |||
1043 | } | |||
1044 | } | |||
1045 | ||||
1046 | //------------------------------------------------------------------ | |||
1047 | // UpdateTrackLabels | |||
1048 | //------------------------------------------------------------------ | |||
1049 | void MyProcessor::UpdateTrackLabels(void) | |||
1050 | { | |||
1051 | // Get the label pointers | |||
1052 | string name, tag; | |||
1053 | map<string, vector<TGLabel*> > &thrownlabs = hdvmf->GetThrownLabels(); | |||
1054 | map<string, vector<TGLabel*> > &reconlabs = hdvmf->GetReconstructedLabels(); | |||
1055 | hdvmf->GetReconFactory(name, tag); | |||
1056 | ||||
1057 | // Get Thrown particles | |||
1058 | vector<const DMCThrown*> throwns; | |||
1059 | if(loop)loop->Get(throwns); | |||
1060 | ||||
1061 | // Get the track info as DKinematicData objects | |||
1062 | vector<const DKinematicData*> trks; | |||
1063 | vector<const DKinematicData*> TrksCand; | |||
1064 | vector<const DTrackWireBased*> TrksWireBased; | |||
1065 | vector<const DTrackTimeBased*> TrksTimeBased; | |||
1066 | vector<const DTrackCandidate*> cand; | |||
1067 | if(loop)loop->Get(cand); | |||
1068 | for(unsigned int i=0; i<cand.size(); i++)TrksCand.push_back(cand[i]); | |||
1069 | ||||
1070 | if(loop)loop->Get(TrksWireBased); | |||
1071 | if(loop)loop->Get(TrksTimeBased); | |||
1072 | ||||
1073 | if(name=="DChargedTrack"){ | |||
1074 | vector<const DChargedTrack*> chargedtracks; | |||
1075 | if(loop)loop->Get(chargedtracks, tag.c_str()); | |||
1076 | for(unsigned int i=0; i<chargedtracks.size(); i++){ | |||
1077 | trks.push_back(chargedtracks[i]->Get_BestFOM()); | |||
1078 | } | |||
1079 | } | |||
1080 | if(name=="DTrackTimeBased"){ | |||
1081 | vector<const DTrackTimeBased*> timetracks; | |||
1082 | if(loop)loop->Get(timetracks, tag.c_str()); | |||
1083 | for(unsigned int i=0; i<timetracks.size(); i++)trks.push_back(timetracks[i]); | |||
1084 | } | |||
1085 | if(name=="DTrackWireBased"){ | |||
1086 | vector<const DTrackWireBased*> wiretracks; | |||
1087 | if(loop)loop->Get(wiretracks, tag.c_str()); | |||
1088 | for(unsigned int i=0; i<wiretracks.size(); i++)trks.push_back(wiretracks[i]); | |||
1089 | } | |||
1090 | if(name=="DTrackCandidate"){ | |||
1091 | vector<const DTrackCandidate*> candidates; | |||
1092 | if(loop)loop->Get(candidates, tag.c_str()); | |||
1093 | for(unsigned int i=0; i<candidates.size(); i++)trks.push_back(candidates[i]); | |||
1094 | } | |||
1095 | if(name=="DNeutralParticle"){ | |||
1096 | vector<const DNeutralParticle*> photons; | |||
1097 | if(loop)loop->Get(photons, tag.c_str()); | |||
1098 | for(unsigned int i=0; i<photons.size(); i++) { | |||
1099 | trks.push_back(photons[i]->Get_BestFOM()); | |||
1100 | } | |||
1101 | } | |||
1102 | if(name=="DTwoGammaFit"){ | |||
1103 | vector<const DTwoGammaFit*> twogammafits; | |||
1104 | if(loop)loop->Get(twogammafits, tag.c_str()); | |||
1105 | for(unsigned int i=0; i<twogammafits.size(); i++)trks.push_back(twogammafits[i]); | |||
1106 | } | |||
1107 | ||||
1108 | // Clear all labels (i.e. draw ---- in them) | |||
1109 | map<string, vector<TGLabel*> >::iterator iter; | |||
1110 | for(iter=reconlabs.begin(); iter!=reconlabs.end(); iter++){ | |||
1111 | vector<TGLabel*> &labs = iter->second; | |||
1112 | for(unsigned int i=1; i<labs.size(); i++){ | |||
1113 | labs[i]->SetText("--------"); | |||
1114 | } | |||
1115 | } | |||
1116 | for(iter=thrownlabs.begin(); iter!=thrownlabs.end(); iter++){ | |||
1117 | vector<TGLabel*> &labs = iter->second; | |||
1118 | for(unsigned int i=1; i<labs.size(); i++){ | |||
1119 | labs[i]->SetText("--------"); | |||
1120 | } | |||
1121 | } | |||
1122 | ||||
1123 | // Loop over thrown particles and fill in labels | |||
1124 | int ii=0; | |||
1125 | for(unsigned int i=0; i<throwns.size(); i++){ | |||
1126 | const DMCThrown *trk = throwns[i]; | |||
1127 | if(trk->type==0)continue; | |||
1128 | int row = thrownlabs["trk"].size()-(ii++)-1; | |||
1129 | if(row<1)break; | |||
1130 | ||||
1131 | stringstream trkno, type, p, theta, phi, z; | |||
1132 | trkno<<setprecision(4)<<i+1; | |||
1133 | thrownlabs["trk"][row]->SetText(trkno.str().c_str()); | |||
1134 | ||||
1135 | thrownlabs["type"][row]->SetText(ParticleType((Particle_t)trk->type)); | |||
1136 | ||||
1137 | p<<setprecision(4)<<trk->momentum().Mag(); | |||
1138 | thrownlabs["p"][row]->SetText(p.str().c_str()); | |||
1139 | ||||
1140 | theta<<setprecision(4)<<trk->momentum().Theta()*TMath::RadToDeg(); | |||
1141 | thrownlabs["theta"][row]->SetText(theta.str().c_str()); | |||
1142 | ||||
1143 | double myphi = trk->momentum().Phi(); | |||
1144 | if(myphi<0.0)myphi+=2.0*M_PI3.14159265358979323846; | |||
1145 | phi<<setprecision(4)<<myphi; | |||
1146 | thrownlabs["phi"][row]->SetText(phi.str().c_str()); | |||
1147 | ||||
1148 | z<<setprecision(4)<<trk->position().Z(); | |||
1149 | thrownlabs["z"][row]->SetText(z.str().c_str()); | |||
1150 | } | |||
1151 | ||||
1152 | // Loop over tracks and fill in labels | |||
1153 | for(unsigned int i=0; i<trks.size(); i++){ | |||
1154 | const DKinematicData *trk = trks[i]; | |||
1155 | int row = reconlabs["trk"].size()-i-1; | |||
1156 | if(row<1)break; | |||
1157 | ||||
1158 | stringstream trkno, type, p, theta, phi, z, chisq_per_dof, Ndof,cand; | |||
1159 | stringstream fom; | |||
1160 | trkno<<setprecision(4)<<i+1; | |||
1161 | reconlabs["trk"][row]->SetText(trkno.str().c_str()); | |||
1162 | ||||
1163 | double mass = trk->mass(); | |||
1164 | if(fabs(mass-0.13957)<1.0E-4)type<<"pi"; | |||
1165 | else if(fabs(mass-0.93827)<1.0E-4)type<<"proton"; | |||
1166 | else if(fabs(mass-0.493677)<1.0E-4)type<<"K"; | |||
1167 | else if(fabs(mass-0.000511)<1.0E-4)type<<"e"; | |||
1168 | else if (fabs(mass)<1.e-4 && fabs(trk->charge())<1.e-4) type << "gamma"; | |||
1169 | else type<<"q="; | |||
1170 | if (fabs(trk->charge())>1.e-4){ | |||
1171 | type<<(trk->charge()>0 ? "+":"-"); | |||
1172 | } | |||
1173 | reconlabs["type"][row]->SetText(type.str().c_str()); | |||
1174 | ||||
1175 | p<<setprecision(4)<<trk->momentum().Mag(); | |||
1176 | reconlabs["p"][row]->SetText(p.str().c_str()); | |||
1177 | ||||
1178 | theta<<setprecision(4)<<trk->momentum().Theta()*TMath::RadToDeg(); | |||
1179 | reconlabs["theta"][row]->SetText(theta.str().c_str()); | |||
1180 | ||||
1181 | double myphi = trk->momentum().Phi(); | |||
1182 | if(myphi<0.0)myphi+=2.0*M_PI3.14159265358979323846; | |||
1183 | phi<<setprecision(4)<<myphi; | |||
1184 | reconlabs["phi"][row]->SetText(phi.str().c_str()); | |||
1185 | ||||
1186 | z<<setprecision(4)<<trk->position().Z(); | |||
1187 | reconlabs["z"][row]->SetText(z.str().c_str()); | |||
1188 | ||||
1189 | // Get chisq and Ndof for DTrackTimeBased or DTrackWireBased objects | |||
1190 | const DTrackTimeBased *timetrack=dynamic_cast<const DTrackTimeBased*>(trk); | |||
1191 | const DTrackWireBased *track=dynamic_cast<const DTrackWireBased*>(trk); | |||
1192 | const DTwoGammaFit *twogammafit=dynamic_cast<const DTwoGammaFit*>(trk); | |||
1193 | ||||
1194 | if(timetrack){ | |||
1195 | chisq_per_dof<<setprecision(4)<<timetrack->chisq/timetrack->Ndof; | |||
1196 | Ndof<<timetrack->Ndof; | |||
1197 | fom << timetrack->FOM; | |||
1198 | }else if(track){ | |||
1199 | chisq_per_dof<<setprecision(4)<<track->chisq/track->Ndof; | |||
1200 | Ndof<<track->Ndof; | |||
1201 | fom << "N/A"; | |||
1202 | }else if(twogammafit){ | |||
1203 | chisq_per_dof<<setprecision(4)<<twogammafit->getChi2(); | |||
1204 | Ndof<<twogammafit->getNdf(); | |||
1205 | fom << twogammafit->getProb(); | |||
1206 | }else{ | |||
1207 | chisq_per_dof<<"N/A"; | |||
1208 | Ndof<<"N/A"; | |||
1209 | fom << "N/A"; | |||
1210 | } | |||
1211 | reconlabs["chisq/Ndof"][row]->SetText(chisq_per_dof.str().c_str()); | |||
1212 | reconlabs["Ndof"][row]->SetText(Ndof.str().c_str()); | |||
1213 | reconlabs["FOM"][row]->SetText(fom.str().c_str()); | |||
1214 | ||||
1215 | if (timetrack){ | |||
1216 | cand << timetrack->candidateid; | |||
1217 | } | |||
1218 | else if (track){ | |||
1219 | cand << track->candidateid; | |||
1220 | } | |||
1221 | else { | |||
1222 | cand << "--------"; | |||
1223 | } | |||
1224 | reconlabs["cand"][row]->SetText(cand.str().c_str()); | |||
1225 | } | |||
1226 | ||||
1227 | // Have the pop-up window with the full particle list update it's labels | |||
1228 | fulllistmf->UpdateTrackLabels(throwns, trks); | |||
1229 | debugermf->SetTrackCandidates(TrksCand); | |||
1230 | debugermf->SetTrackWireBased(TrksWireBased); | |||
1231 | debugermf->SetTrackTimeBased(TrksTimeBased); | |||
1232 | debugermf->UpdateTrackLabels(); | |||
1233 | } | |||
1234 | ||||
1235 | //------------------------------------------------------------------ | |||
1236 | // AddKinematicDataTrack | |||
1237 | //------------------------------------------------------------------ | |||
1238 | void MyProcessor::AddKinematicDataTrack(const DKinematicData* kd, int color, double size) | |||
1239 | { | |||
1240 | // Create a reference trajectory with the given kinematic data and swim | |||
1241 | // it through the detector. | |||
1242 | DReferenceTrajectory rt(Bfield); | |||
1243 | ||||
1244 | if(MATERIAL_MAP_MODEL=="DRootGeom"){ | |||
1245 | rt.SetDRootGeom(RootGeom); | |||
1246 | rt.SetDGeometry(NULL__null); | |||
1247 | }else if(MATERIAL_MAP_MODEL=="DGeometry"){ | |||
1248 | rt.SetDRootGeom(NULL__null); | |||
1249 | rt.SetDGeometry(geom); | |||
1250 | }else if(MATERIAL_MAP_MODEL!="NONE"){ | |||
1251 | _DBG_std::cerr<<"MyProcessor.cc"<<":"<<1251<< " "<<"WARNING: Invalid value for TRKFIT:MATERIAL_MAP_MODEL (=\""<<MATERIAL_MAP_MODEL<<"\")"<<endl; | |||
1252 | } | |||
1253 | ||||
1254 | rt.SetMass(kd->mass()); | |||
1255 | rt.Swim(kd->position(), kd->momentum(), kd->charge()); | |||
1256 | ||||
1257 | // Create a new graphics set and fill it with all of the trajectory points | |||
1258 | DGraphicSet gset(color, kLine, size); | |||
1259 | DReferenceTrajectory::swim_step_t *step = rt.swim_steps; | |||
1260 | for(int i=0; i<rt.Nswim_steps; i++, step++){ | |||
1261 | TVector3 tpoint(step->origin.X(),step->origin.Y(),step->origin.Z()); | |||
1262 | gset.points.push_back(tpoint); | |||
1263 | } | |||
1264 | ||||
1265 | // Push the graphics set onto the stack | |||
1266 | graphics.push_back(gset); | |||
1267 | } | |||
1268 | ||||
1269 | //------------------------------------------------------------------ | |||
1270 | // GetIntersectionWithCalorimeter | |||
1271 | //------------------------------------------------------------------ | |||
1272 | void MyProcessor::GetIntersectionWithCalorimeter(const DKinematicData* kd, DVector3 &pos, DetectorSystem_t &who) | |||
1273 | { | |||
1274 | // Create a reference trajectory with the given kinematic data and swim | |||
1275 | // it through the detector. | |||
1276 | DReferenceTrajectory rt(Bfield); | |||
1277 | ||||
1278 | if(MATERIAL_MAP_MODEL=="DRootGeom"){ | |||
1279 | rt.SetDRootGeom(RootGeom); | |||
1280 | rt.SetDGeometry(NULL__null); | |||
1281 | }else if(MATERIAL_MAP_MODEL=="DGeometry"){ | |||
1282 | rt.SetDRootGeom(NULL__null); | |||
1283 | rt.SetDGeometry(geom); | |||
1284 | }else if(MATERIAL_MAP_MODEL!="NONE"){ | |||
1285 | _DBG_std::cerr<<"MyProcessor.cc"<<":"<<1285<< " "<<"WARNING: Invalid value for TRKFIT:MATERIAL_MAP_MODEL (=\""<<MATERIAL_MAP_MODEL<<"\")"<<endl; | |||
1286 | } | |||
1287 | ||||
1288 | rt.SetMass(kd->mass()); | |||
1289 | rt.Swim(kd->position(), kd->momentum(), kd->charge()); | |||
1290 | ||||
1291 | // Find intersection with FCAL | |||
1292 | DVector3 pos_fcal; | |||
1293 | double s_fcal = 1.0E6; | |||
1294 | DVector3 origin(0.0, 0.0, FCAL_Zmin); | |||
1295 | DVector3 norm(0.0, 0.0, -1.0); | |||
1296 | rt.GetIntersectionWithPlane(origin, norm, pos_fcal, &s_fcal); | |||
1297 | if(pos_fcal.Perp()<FCAL_Rmin || pos_fcal.Perp()>FCAL_Rmax)s_fcal = 1.0E6; | |||
1298 | ||||
1299 | // Find intersection with BCAL | |||
1300 | DVector3 pos_bcal; | |||
1301 | double s_bcal = 1.0E6; | |||
1302 | rt.GetIntersectionWithRadius(BCAL_Rmin, pos_bcal, &s_bcal); | |||
1303 | if(pos_bcal.Z()<BCAL_Zmin || pos_bcal.Z()>(BCAL_Zmin+BCAL_Zlen))s_bcal = 1.0E6; | |||
1304 | ||||
1305 | if(s_fcal>1000.0 && s_bcal>1000.0){ | |||
1306 | // neither calorimeter hit | |||
1307 | who = SYS_NULL; | |||
1308 | pos.SetXYZ(0.0,0.0,0.0); | |||
1309 | }else if(s_fcal<s_bcal){ | |||
1310 | // FCAL hit | |||
1311 | who = SYS_FCAL; | |||
1312 | pos = pos_fcal; | |||
1313 | }else{ | |||
1314 | // BCAL hit | |||
1315 | who = SYS_BCAL; | |||
1316 | pos = pos_bcal; | |||
1317 | } | |||
1318 | } | |||
1319 | ||||
1320 | //------------------------------------------------------------------ | |||
1321 | // GetFactoryNames | |||
1322 | //------------------------------------------------------------------ | |||
1323 | void MyProcessor::GetFactoryNames(vector<string> &facnames) | |||
1324 | { | |||
1325 | vector<JEventLoop*> loops = app->GetJEventLoops(); | |||
1326 | if(loops.size()>0){ | |||
1327 | vector<string> facnames; | |||
1328 | loops[0]->GetFactoryNames(facnames); | |||
1329 | } | |||
1330 | } | |||
1331 | ||||
1332 | //------------------------------------------------------------------ | |||
1333 | // GetFactories | |||
1334 | //------------------------------------------------------------------ | |||
1335 | void MyProcessor::GetFactories(vector<JFactory_base*> &factories) | |||
1336 | { | |||
1337 | vector<JEventLoop*> loops = app->GetJEventLoops(); | |||
1338 | if(loops.size()>0){ | |||
1339 | factories = loops[0]->GetFactories(); | |||
1340 | } | |||
1341 | } | |||
1342 | ||||
1343 | //------------------------------------------------------------------ | |||
1344 | // GetNrows | |||
1345 | //------------------------------------------------------------------ | |||
1346 | unsigned int MyProcessor::GetNrows(const string &factory, string tag) | |||
1347 | { | |||
1348 | vector<JEventLoop*> loops = app->GetJEventLoops(); | |||
1349 | if(loops.size()>0){ | |||
1350 | // Here is something a little tricky. The GetFactory() method of JEventLoop | |||
1351 | // gets the factory of the specified data name and tag, but without trying | |||
1352 | // to substitute a user-specified tag (a'la -PDEFTAG:XXX=YYY) as is done | |||
1353 | // on normal Get() method calls. Therefore, we have to check for the default | |||
1354 | // tags ourselves and substitute it "by hand". | |||
1355 | if(tag==""){ | |||
1356 | map<string,string> default_tags = loops[0]->GetDefaultTags(); | |||
1357 | map<string, string>::const_iterator iter = default_tags.find(factory); | |||
1358 | if(iter!=default_tags.end())tag = iter->second.c_str(); | |||
1359 | } | |||
1360 | JFactory_base *fac = loops[0]->GetFactory(factory, tag.c_str()); | |||
1361 | ||||
1362 | // Since calling GetNrows will cause the program to quit if there is | |||
1363 | // not a valid event, then first check that there is one before calling it | |||
1364 | if(loops[0]->GetJEvent().GetJEventSource() == NULL__null)return 0; | |||
1365 | ||||
1366 | return fac==NULL__null ? 0:(unsigned int)fac->GetNrows(); | |||
1367 | } | |||
1368 | ||||
1369 | return 0; | |||
1370 | } | |||
1371 | ||||
1372 | //------------------------------------------------------------------ | |||
1373 | // GetDReferenceTrajectory | |||
1374 | //------------------------------------------------------------------ | |||
1375 | void MyProcessor::GetDReferenceTrajectory(string dataname, string tag, unsigned int index, DReferenceTrajectory* &rt, vector<const DCDCTrackHit*> &cdchits) | |||
1376 | { | |||
1377 | _DBG__std::cerr<<"MyProcessor.cc"<<":"<<1377<< std::endl; | |||
1378 | // initialize rt to NULL in case we don't find the one requested | |||
1379 | rt = NULL__null; | |||
1380 | cdchits.clear(); | |||
1381 | ||||
1382 | // Get pointer to the JEventLoop so we can get at the data | |||
1383 | vector<JEventLoop*> loops = app->GetJEventLoops(); | |||
1384 | if(loops.size()==0)return; | |||
1385 | JEventLoop* &loop = loops[0]; | |||
1386 | ||||
1387 | // Variables to hold track parameters | |||
1388 | DVector3 pos, mom(0,0,0); | |||
1389 | double q=0.0; | |||
1390 | double mass; | |||
1391 | ||||
1392 | // Find the specified track | |||
1393 | if(dataname=="DChargedTrack"){ | |||
1394 | vector<const DChargedTrack*> chargedtracks; | |||
1395 | vector<const DTrackTimeBased*> timebasedtracks; | |||
1396 | loop->Get(chargedtracks, tag.c_str()); | |||
1397 | if(index>=chargedtracks.size())return; | |||
1398 | q = chargedtracks[index]->Get_Charge(); | |||
1399 | pos = chargedtracks[index]->Get_BestFOM()->position(); | |||
1400 | mom = chargedtracks[index]->Get_BestFOM()->momentum(); | |||
1401 | chargedtracks[index]->Get_BestFOM()->GetT(timebasedtracks); | |||
1402 | timebasedtracks[0]->Get(cdchits); | |||
1403 | mass = chargedtracks[index]->Get_BestFOM()->mass(); | |||
1404 | } | |||
1405 | ||||
1406 | if(dataname=="DTrackTimeBased"){ | |||
1407 | vector<const DTrackTimeBased*> timetracks; | |||
1408 | loop->Get(timetracks, tag.c_str()); | |||
1409 | if(index>=timetracks.size())return; | |||
1410 | q = timetracks[index]->charge(); | |||
1411 | pos = timetracks[index]->position(); | |||
1412 | mom = timetracks[index]->momentum(); | |||
1413 | timetracks[index]->Get(cdchits); | |||
1414 | mass = timetracks[index]->mass(); | |||
1415 | } | |||
1416 | ||||
1417 | if(dataname=="DTrackWireBased"){ | |||
1418 | vector<const DTrackWireBased*> wiretracks; | |||
1419 | loop->Get(wiretracks, tag.c_str()); | |||
1420 | if(index>=wiretracks.size())return; | |||
1421 | q = wiretracks[index]->charge(); | |||
1422 | pos = wiretracks[index]->position(); | |||
1423 | mom = wiretracks[index]->momentum(); | |||
1424 | wiretracks[index]->Get(cdchits); | |||
1425 | mass = wiretracks[index]->mass(); | |||
1426 | } | |||
1427 | ||||
1428 | if(dataname=="DTrackCandidate"){ | |||
1429 | vector<const DTrackCandidate*> tracks; | |||
1430 | loop->Get(tracks, tag.c_str()); | |||
1431 | if(index>=tracks.size())return; | |||
1432 | q = tracks[index]->charge(); | |||
1433 | pos = tracks[index]->position(); | |||
1434 | mom = tracks[index]->momentum(); | |||
1435 | tracks[index]->Get(cdchits); | |||
1436 | mass = tracks[index]->mass(); | |||
1437 | } | |||
1438 | ||||
1439 | if(dataname=="DMCThrown"){ | |||
1440 | vector<const DMCThrown*> tracks; | |||
1441 | loop->Get(tracks, tag.c_str()); | |||
1442 | if(index>=tracks.size())return; | |||
1443 | const DMCThrown *t = tracks[index]; | |||
1444 | q = t->charge(); | |||
1445 | pos = t->position(); | |||
1446 | mom = t->momentum(); | |||
1447 | tracks[index]->Get(cdchits); | |||
1448 | mass = tracks[index]->mass(); | |||
1449 | _DBG_std::cerr<<"MyProcessor.cc"<<":"<<1449<< " "<<"mass="<<mass<<endl; | |||
1450 | } | |||
1451 | ||||
1452 | // Make sure we found a charged particle we can track | |||
1453 | if(q==0.0 || mom.Mag()<0.01)return; | |||
1454 | ||||
1455 | // Create a new DReference trajectory object. The caller takes | |||
1456 | // ownership of this and so they are responsible for deleting it. | |||
1457 | rt = new DReferenceTrajectory(Bfield); | |||
1458 | if(MATERIAL_MAP_MODEL=="DRootGeom"){ | |||
1459 | rt->SetDRootGeom(RootGeom); | |||
1460 | rt->SetDGeometry(NULL__null); | |||
1461 | }else if(MATERIAL_MAP_MODEL=="DGeometry"){ | |||
1462 | rt->SetDRootGeom(NULL__null); | |||
1463 | rt->SetDGeometry(geom); | |||
1464 | }else if(MATERIAL_MAP_MODEL!="NONE"){ | |||
1465 | _DBG_std::cerr<<"MyProcessor.cc"<<":"<<1465<< " "<<"WARNING: Invalid value for TRKFIT:MATERIAL_MAP_MODEL (=\""<<MATERIAL_MAP_MODEL<<"\")"<<endl; | |||
1466 | } | |||
1467 | rt->Swim(pos, mom, q); | |||
1468 | } | |||
1469 | ||||
1470 | //------------------------------------------------------------------ | |||
1471 | // GetAllWireHits | |||
1472 | //------------------------------------------------------------------ | |||
1473 | void MyProcessor::GetAllWireHits(vector<pair<const DCoordinateSystem*,double> > &allhits) | |||
1474 | { | |||
1475 | /// Argument is vector of pairs that contain a pointer to the | |||
1476 | /// DCoordinateSystem representing a wire and a double that | |||
1477 | /// represents the drift distance. To get info on the specific | |||
1478 | /// wire, one needs to attempt a dynamic_cast to both a DCDCWire | |||
1479 | /// and a DFDCWire and access the parameters of whichever one succeeds. | |||
1480 | ||||
1481 | // Get pointer to the JEventLoop so we can get at the data | |||
1482 | vector<JEventLoop*> loops = app->GetJEventLoops(); | |||
1483 | if(loops.size()==0)return; | |||
1484 | JEventLoop* &loop = loops[0]; | |||
1485 | ||||
1486 | // Get CDC wire hits | |||
1487 | vector<const DCDCTrackHit*> cdchits; | |||
1488 | loop->Get(cdchits); | |||
1489 | for(unsigned int i=0; i<cdchits.size(); i++){ | |||
1490 | pair<const DCoordinateSystem*,double> hit; | |||
1491 | hit.first = cdchits[i]->wire; | |||
1492 | hit.second = cdchits[i]->dist; | |||
1493 | allhits.push_back(hit); | |||
1494 | } | |||
1495 | ||||
1496 | // Get FDC wire hits | |||
1497 | vector<const DFDCPseudo*> fdchits; | |||
1498 | loop->Get(fdchits); | |||
1499 | for(unsigned int i=0; i<fdchits.size(); i++){ | |||
1500 | pair<const DCoordinateSystem*,double> hit; | |||
1501 | hit.first = fdchits[i]->wire; | |||
1502 | hit.second = 0.0055*fdchits[i]->time; | |||
1503 | allhits.push_back(hit); | |||
1504 | } | |||
1505 | } | |||
1506 | ||||
1507 | ||||
1508 |