Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
hdview2/MyProcessor.cc
Go to the documentation of this file.
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 #include <TLegend.h>
22 
23 #include <particleType.h>
24 #include "hdview2.h"
25 #include "hdv_mainframe.h"
26 #include "hdv_debugerframe.h"
27 #include "MyProcessor.h"
28 #include "TRACKING/DTrackHit.h"
29 #include "TRACKING/DQuickFit.h"
32 #include "TRACKING/DMCTrackHit.h"
33 #include "TRACKING/DMCThrown.h"
36 #include "PID/DChargedTrack.h"
38 #include "JANA/JGeometry.h"
40 #include "FCAL/DFCALHit.h"
41 #include "CCAL/DCCALHit.h"
42 #include "CCAL/DCCALShower.h"
43 #include "TOF/DTOFGeometry.h"
44 #include "TOF/DTOFHit.h"
45 #include "TOF/DTOFTDCDigiHit.h"
46 #include "TOF/DTOFDigiHit.h"
47 #include "TOF/DTOFPaddleHit.h"
48 #include "TOF/DTOFPoint.h"
49 #include "FDC/DFDCGeometry.h"
50 #include "CDC/DCDCTrackHit.h"
51 #include "FDC/DFDCPseudo.h"
52 #include "FDC/DFDCIntersection.h"
53 #include "HDGEOMETRY/DGeometry.h"
54 #include "FCAL/DFCALGeometry.h"
55 #include "FCAL/DFCALHit.h"
56 #include "PID/DNeutralParticle.h"
57 #include "PID/DNeutralShower.h"
58 #include "BCAL/DBCALHit.h"
60 #include "TOF/DTOFPoint.h"
61 #include "START_COUNTER/DSCHit.h"
62 #include "DVector2.h"
63 #include "TRIGGER/DL1Trigger.h"
64 
65 extern hdv_mainframe *hdvmf;
66 
67 // These are declared in hdv_mainframe.cc, but as static so we need to do it here as well (yechh!)
68 static float FCAL_Zmin = 622.8;
69 static float FCAL_Rmin = 6.0;
70 static float FCAL_Rmax = 212.0/2.0;
71 static float BCAL_Rmin = 65.0;
72 static float BCAL_Zlen = 390.0;
73 static float BCAL_Zmin = 212.0 - BCAL_Zlen/2.0;
74 
75 
76 static vector<vector <DFDCWire *> >fdcwires;
77 
78 
80  // sort by track number and then by particle type, then by z-coordinate
81  // (yes, I saw the same track number have different particle types!)
82  if(a->track != b->track)return a->track < b->track;
83  if(a->part != b->part )return a->part < b->part;
84  return a->z < b->z;
85 }
86 
87 
89 
90 //------------------------------------------------------------------
91 // MyProcessor
92 //------------------------------------------------------------------
94 {
95  Bfield = NULL;
96  loop = NULL;
97 
98  // Tell factory to keep around a few density histos
99  //gPARMS->SetParameter("TRKFIND:MAX_DEBUG_BUFFERS", 16);
100 
101  RMAX_INTERIOR = 65.0;
102  RMAX_EXTERIOR = 88.0;
103  gPARMS->SetDefaultParameter("RT:RMAX_INTERIOR", RMAX_INTERIOR, "cm track drawing Rmax inside solenoid region");
104  gPARMS->SetDefaultParameter("RT:RMAX_EXTERIOR", RMAX_EXTERIOR, "cm track drawing Rmax outside solenoid region");
105 
106  BCALVERBOSE = 0;
107  gPARMS->SetDefaultParameter("BCALVERBOSE", BCALVERBOSE, "Verbosity level for BCAL objects and display");
108 
109  gMYPROC = this;
110 }
111 
112 //------------------------------------------------------------------
113 // ~MyProcessor
114 //------------------------------------------------------------------
116 {
117 
118 }
119 
120 //------------------------------------------------------------------
121 // init
122 //------------------------------------------------------------------
123 jerror_t MyProcessor::init(void)
124 {
125  // Make sure detectors have been drawn
126  //if(!drew_detectors)DrawDetectors();
127 
128 
129 
130  vector<JEventLoop*> loops = app->GetJEventLoops();
131  if(loops.size()>0){
132 
133  vector<string> facnames;
134  loops[0]->GetFactoryNames(facnames);
135 
136  hdvmf = new hdv_mainframe(gClient->GetRoot(), 1400, 700);
137  hdvmf->SetCandidateFactories(facnames);
140  hdvmf->SetReconstructedFactories(facnames);
141  hdvmf->SetChargedTrackFactories(facnames);
142  fulllistmf = hdvmf->GetFullListFrame();
143  debugermf = hdvmf->GetDebugerFrame();
144  BCALHitCanvas = hdvmf->GetBcalDispFrame();
145 
146  if (BCALHitCanvas){
147  BCALHitMatrixU = new TH2F("BCALHitMatrixU","BCAL Hits Upstream", 48*4+2, -1.5, 192.5, 10, 0., 10.);
148  BCALHitMatrixD = new TH2F("BCALHitMatrixD","BCAL Hits Downstream",48*4+2, -1.5, 192.5, 10, 0., 10.);
149  BCALParticles = new TH2F("BCALParticles","BCAL Hits Downstream",(48*4+2)*4, -1.87, 361.87, 1, 0., 1.);
150  BCALHitMatrixU->SetStats(0);
151  BCALHitMatrixD->SetStats(0);
152  BCALParticles->SetStats(0);
153  BCALHitMatrixU->GetXaxis()->SetTitle("Sector number");
154  BCALHitMatrixD->GetXaxis()->SetTitle("Sector number");
155  BCALParticles->GetXaxis()->SetTitle("Phi angle [deg]");
156  }
157  }
158 
159  return NOERROR;
160 }
161 
162 //------------------------------------------------------------------
163 // brun
164 //------------------------------------------------------------------
165 jerror_t MyProcessor::brun(JEventLoop *eventloop, int32_t runnumber)
166 {
167 
168  // Read in Magnetic field map
169  DApplication* dapp = dynamic_cast<DApplication*>(eventloop->GetJApplication());
170  Bfield = dapp->GetBfield(runnumber);
171  const DGeometry *dgeom = dapp->GetDGeometry(runnumber);
172  dgeom->GetFDCWires(fdcwires);
173 
174  RootGeom = dapp->GetRootGeom(runnumber);
175  geom = dapp->GetDGeometry(runnumber);
176 
177  MATERIAL_MAP_MODEL="DGeometry";
178  gPARMS->SetDefaultParameter("TRKFIT:MATERIAL_MAP_MODEL", MATERIAL_MAP_MODEL);
179 
180  eventloop->GetCalib("PID/photon_track_matching", photon_track_matching);
181  DELTA_R_FCAL = photon_track_matching["DELTA_R_FCAL"];
182 
183  return NOERROR;
184 }
185 
186 //------------------------------------------------------------------
187 // evnt
188 //------------------------------------------------------------------
189 jerror_t MyProcessor::evnt(JEventLoop *eventLoop, uint64_t eventnumber)
190 {
191  if(!eventLoop)return NOERROR;
192  loop = eventLoop;
193  last_jevent.FreeEvent();
194  last_jevent = loop->GetJEvent();
195 
196  // get the trigger bits
197  char trigstring[10];
198  const DL1Trigger *trig = NULL;
199  try {
200  loop->GetSingle(trig);
201  } catch (...) {}
202  if (trig) {
203  sprintf(trigstring,"0x%04x",trig->trig_mask);
204  } else {
205  sprintf(trigstring,"no bits");
206  }
207  hdvmf->SetTrig(trigstring);
208 
209  string source = "<no source>";
210  if(last_jevent.GetJEventSource())source = last_jevent.GetJEventSource()->GetSourceName();
211 
212  cout<<"----------- New Event "<<eventnumber<<" (run " << last_jevent.GetRunNumber()<<") -------------"<<endl;
213  hdvmf->SetEvent(eventnumber);
214  hdvmf->SetRun(last_jevent.GetRunNumber());
215  hdvmf->SetSource(source.c_str());
216  hdvmf->DoMyRedraw();
217 
218  japp->SetSequentialEventComplete();
219 
220  return NOERROR;
221 }
222 
223 //------------------------------------------------------------------
224 // FillGraphics
225 //------------------------------------------------------------------
227 {
228  /// Create "graphics" objects for this event given the current GUI settings.
229  ///
230  /// This method will create DGraphicSet objects that represent tracks, hits,
231  /// and showers for the event. It creates objects for both hits and
232  /// reconstructed entities. The "graphics" objects created here are
233  /// really just collections of 3D space points with flags indicating
234  /// whether they should be drawn as markers or lines and with what
235  /// color and size. The actual graphics objects are created for the
236  /// various views of the detector in hdv_mainframe.
237 
238  graphics.clear();
239  graphics_xyA.clear(); // The objects placed in these will be deleted by hdv_mainframe
240  graphics_xyB.clear(); // The objects placed in these will be deleted by hdv_mainframe
241  graphics_xz.clear(); // The objects placed in these will be deleted by hdv_mainframe
242  graphics_yz.clear(); // The objects placed in these will be deleted by hdv_mainframe
243  graphics_tof_hits.clear(); // The objects placed in these will be deleted by hdv_mainframe
244 
245  if(!loop)return;
246 
247  vector<const DSCHit *>schits;
248  loop->Get(schits);
249  vector<const DTrackCandidate*> trCand;
250  loop->Get(trCand);
251  vector<const DTrackTimeBased*> trTB;
252  loop->Get(trTB);
253  vector<const DTrackWireBased*> trWB;
254  loop->Get(trWB);
256  p->SetNTrCand(trCand.size());
257  p->SetNTrWireBased(trWB.size());
258  p->SetNTrTimeBased(trTB.size());
259 
260  if (BCALHitCanvas) {
261  vector<const DBCALHit*> locBcalHits;
262  loop->Get(locBcalHits);
263  BCALHitMatrixU->Reset();
264  BCALHitMatrixD->Reset();
265  for (unsigned int k=0;k<locBcalHits.size();k++){
266 
267  const DBCALHit* hit = locBcalHits[k];
268  float idxY = (float)hit->layer-1;
269  float idxX = (float) (hit->sector-1 + (hit->module-1)*4);
270  if (hit->end == DBCALGeometry::kUpstream){
271  if (hit->layer==1){
272  BCALHitMatrixU->Fill(idxX,idxY,hit->E);
273  } else if (hit->layer==2){
274  BCALHitMatrixU->Fill(idxX,idxY,hit->E);
275  BCALHitMatrixU->Fill(idxX,idxY+1.,hit->E);
276  } else if (hit->layer==3){
277  BCALHitMatrixU->Fill(idxX,idxY+1,hit->E);
278  BCALHitMatrixU->Fill(idxX,idxY+2.,hit->E);
279  BCALHitMatrixU->Fill(idxX,idxY+3.,hit->E);
280  } else if (hit->layer==4){
281  BCALHitMatrixU->Fill(idxX,idxY+3,hit->E);
282  BCALHitMatrixU->Fill(idxX,idxY+4.,hit->E);
283  BCALHitMatrixU->Fill(idxX,idxY+5.,hit->E);
284  BCALHitMatrixU->Fill(idxX,idxY+6.,hit->E);
285  }
286  } else {
287  if (hit->layer==1){
288  BCALHitMatrixD->Fill(idxX,idxY,hit->E);
289  } else if (hit->layer==2){
290  BCALHitMatrixD->Fill(idxX,idxY,hit->E);
291  BCALHitMatrixD->Fill(idxX,idxY+1.,hit->E);
292  } else if (hit->layer==3){
293  BCALHitMatrixD->Fill(idxX,idxY+1,hit->E);
294  BCALHitMatrixD->Fill(idxX,idxY+2.,hit->E);
295  BCALHitMatrixD->Fill(idxX,idxY+3.,hit->E);
296  } else if (hit->layer==4){
297  BCALHitMatrixD->Fill(idxX,idxY+3,hit->E);
298  BCALHitMatrixD->Fill(idxX,idxY+4.,hit->E);
299  BCALHitMatrixD->Fill(idxX,idxY+5.,hit->E);
300  BCALHitMatrixD->Fill(idxX,idxY+6.,hit->E);
301  }
302  }
303  }
304 
305  // Fill BCAL points into layer histograms
306  vector<const DBCALPoint*> locBcalPoint;
307  loop->Get(locBcalPoint);
308  for (int layer=0; layer<4; layer++) {
309  BCALPointZphiLayer[layer]->Reset();
310  BCALPointPhiTLayer[layer]->Reset();
311  }
312 
313  vector<double> BCALpoints_z;
314  vector<double> BCALpoints_phi;
315  vector<double>::const_iterator points_min_phi;
316  vector<double>::const_iterator points_max_phi;
317  vector<double>::const_iterator points_min_z;
318  vector<double>::const_iterator points_max_z;
319 
320  for(unsigned int k=0;k<locBcalPoint.size();k++){
321  BCALpoints_z.push_back(locBcalPoint[k]->z());
322  BCALpoints_phi.push_back(locBcalPoint[k]->phi()*TMath::RadToDeg());
323  }
324 
325  if(BCALpoints_phi.size() > 0){
326  points_min_phi = min_element(BCALpoints_phi.begin(),BCALpoints_phi.end());
327  points_max_phi = max_element(BCALpoints_phi.begin(),BCALpoints_phi.end());
328  }
329 
330  if(BCALpoints_z.size() > 0){
331  points_min_z = min_element(BCALpoints_z.begin(),BCALpoints_z.end());
332  points_max_z = max_element(BCALpoints_z.begin(),BCALpoints_z.end());
333  }
334 
335  BCALpoints_z.clear();
336  BCALpoints_phi.clear();
337 
338  for (unsigned int k=0;k<locBcalPoint.size();k++){
339  const DBCALPoint* point = locBcalPoint[k];
340  float pointphi = point->phi()*TMath::RadToDeg();
341  if (pointphi>360) pointphi-=360;
342  if (pointphi<0) pointphi+=360;
343  if (BCALVERBOSE>0) printf(" %3i (m,l,s) = (%2i,%i,%i) (z,phi) = (%6.1f,%4.0f) t=%7.2f E=%6.1f MeV\n",
344  k, point->module(), point->layer(), point->sector(),point->z(),
345  pointphi,point->t(),point->E()*1000);
346  float weight = log(point->E()*1000);
347  BCALPointZphiLayer[point->layer()-1]->Fill(pointphi,point->z(),weight);
348  BCALPointZphiLayer[point->layer()-1]->GetXaxis()->SetRangeUser((*points_min_phi-2),(*points_max_phi+2));
349  BCALPointZphiLayer[point->layer()-1]->GetYaxis()->SetRangeUser((*points_min_z-5),(*points_max_z+5));
350  float time = point->t();
351  if (point->t()>100) time=99; // put overflow in last bin
352  if (point->t()<-100) time=-99;
353  BCALPointPhiTLayer[point->layer()-1]->Fill(pointphi,time,weight);
354  }
355 
356  // Fill BCAL Clusters into histograms
357  vector<const DBCALCluster*> locBcalCluster;
358  loop->Get(locBcalCluster);
359 
360  unsigned int oldsize = BCALClusterZphiHistos.size();
361  for (unsigned int k=0;k<oldsize;k++){
362  delete BCALClusterZphiHistos[k];
363  }
364  BCALClusterZphiHistos.clear();
365 
366  delete ClusterLegend;
367  ClusterLegend = new TLegend(0.91,0.01,0.99,0.99);
368 
369  for(unsigned int k = 0 ; k < locBcalCluster.size() ; k++){
370  const DBCALCluster* cluster_range = locBcalCluster[k];
371 
372  vector<const DBCALPoint*> clusterpoints_range;
373  cluster_range->Get(clusterpoints_range);
374 
375  for(unsigned int l=0;l<clusterpoints_range.size();l++){
376  BCALpoints_z.push_back(clusterpoints_range[l]->z());
377  BCALpoints_phi.push_back(clusterpoints_range[l]->phi()*TMath::RadToDeg());
378  }
379 
380  if(BCALpoints_phi.size() > 0){
381  points_min_phi = min_element(BCALpoints_phi.begin(),BCALpoints_phi.end());
382  points_max_phi = max_element(BCALpoints_phi.begin(),BCALpoints_phi.end());
383  }
384 
385  if(BCALpoints_z.size() > 0){
386  points_min_z = min_element(BCALpoints_z.begin(),BCALpoints_z.end());
387  points_max_z = max_element(BCALpoints_z.begin(),BCALpoints_z.end());
388  }
389 
390  }
391 
392  BCALpoints_z.clear();
393  BCALpoints_phi.clear();
394 
395  for (unsigned int k=0;k<locBcalCluster.size();k++){
396  char name[255];
397  sprintf(name,"ClusterZphi%i",k);
398  BCALClusterZphiHistos.push_back(new TH2F(name,"BCAL Clusters;Phi angle [deg];Z position (cm)",48*4,0,360,48,-80,400));
399  int color=k+1;
400  if (k>3) color=k+2; // remove yellow
401  if (k>7) color=k+3; // remove white
402  FormatHistogram(BCALClusterZphiHistos[k],color);
403  BCALClusterZphiHistos[k]->SetLineWidth(2);
404  sprintf(name,"Cluster %i",k+1);
405  ClusterLegend->AddEntry(BCALClusterZphiHistos[k], name);
406 
407  const DBCALCluster* cluster = locBcalCluster[k];
408  if (BCALVERBOSE>0) printf("cluster %2i, (phi,theta,rho) = (%6.2f,%6.2f,%7.2f) t=%8.3f E=%8.2f MeV\n", k+1,cluster->phi(), cluster->theta(), cluster->rho(), cluster->t(), cluster->E()*1000);
409  vector<const DBCALPoint*> clusterpoints;
410  cluster->Get(clusterpoints);
411 
412  for (unsigned int l=0;l<clusterpoints.size();l++){
413  const DBCALPoint* clusterpoint = clusterpoints[l];
414  float weight = log(clusterpoint->E()*1000);
415  float pointphi = clusterpoint->phi()*TMath::RadToDeg();
416  if (pointphi>360) pointphi-=360;
417  if (pointphi<0) pointphi+=360;
418  if (BCALVERBOSE>1) printf(" (m,l,s) = (%2i,%i,%i) (z,phi) = (%6.1f,%4.0f) t=%7.2f E=%6.1f MeV\n",
419  clusterpoint->module(), clusterpoint->layer(), clusterpoint->sector(),clusterpoint->z(),
420  pointphi,clusterpoint->t(),clusterpoint->E()*1000);
421  BCALClusterZphiHistos[k]->Fill(pointphi,clusterpoint->z(),weight);
422  BCALClusterZphiHistos[k]->GetXaxis()->SetRangeUser((*points_min_phi - 2),(*points_max_phi + 2));
423  BCALClusterZphiHistos[k]->GetYaxis()->SetRangeUser((*points_min_z-10),(*points_max_z+10));
424  }
425  }
426 
427  vector<const DBCALIncidentParticle*> locBcalParticles;
428  loop->Get(locBcalParticles);
429  BCALParticles->Reset();
430  BCALPLables.clear();
431  for (unsigned int k=0;k<locBcalParticles.size();k++){
432  const DBCALIncidentParticle* part = locBcalParticles[k];
433 
434  float p = TMath::Sqrt(part->px*part->px + part->py*part->py + part->pz*part->pz);
435  float phi=999;
436  if (part->x!=0){
437  phi = TMath::ATan(TMath::Abs(part->y/part->x));
438  //cout<<phi<<" "<<part->y<<" / "<< part->x<<endl;
439  if (part->y>0){
440  if (part->x<0.){
441  phi = 3.1415926 - phi;
442  }
443  } else {
444  if (part->x<0){
445  phi += 3.1415926;
446  } else {
447  phi = 3.1415926*2. - phi;
448  }
449  }
450 
451  phi = phi*180./3.1415926;
452  }
453  //cout<<phi<<" "<<p<<endl;
454  BCALParticles->Fill(phi,0.5,p);
455  char l[20];
456  sprintf(l,"%d",part->ptype);
457  TText *t = new TText(phi,1.01,l);
458  t->SetTextSize(0.08);
459  t->SetTextFont(72);
460  t->SetTextAlign(21);
461  BCALPLables.push_back(t);
462  }
463 
464  // Get Maximum bin from cluster histograms
465  float clustermax=0;
466  if (BCALClusterZphiHistos.size()>0) {
467  float max = BCALClusterZphiHistos[0]->GetMaximum();
468  if (max>clustermax) clustermax=max;
469  }
470 
471  BCALHitCanvas->Clear();
472  BCALHitCanvas->Divide(1,3,0.001,0.001);
473  float leftmargin = 0.07;
474  BCALHitCanvas->cd(1);
475  gPad->SetGridx();
476  gPad->SetGridy();
477  gPad->SetLeftMargin(leftmargin);
478  // if (BCALHitMatrixU->GetMaximum() < 1) {
479  // BCALHitMatrixU->Scale(1000); // Scale histogram to MeV
480  // BCALHitMatrixU->GetZaxis()->SetTitle("Energy (MeV)");
481  // } else {
482  // BCALHitMatrixU->GetZaxis()->SetTitle("Energy (GeV)");
483  // }
484  // BCALHitMatrixU->Draw("colz");
485 
486  bool firstup=1;
487 
488  for (int layer=0; layer<4; layer++) {
489  if (BCALPointPhiTLayer[layer]->GetEntries()>0) {
490  BCALPointPhiTLayer[layer]->SetMaximum(clustermax);
491  if (firstup==1) {
492  BCALPointPhiTLayer[layer]->Draw("box");
493  firstup=0;
494  } else {
495  BCALPointPhiTLayer[layer]->Draw("box,same");
496  }
497  }
498  }
499  LayerLegend->Draw();
500 
501  BCALHitCanvas->cd(2);
502  gPad->SetGridx();
503  gPad->SetGridy();
504  gPad->SetLeftMargin(leftmargin);
505  // if (BCALHitMatrixD->GetMaximum() < 1) {
506  // BCALHitMatrixD->Scale(1000); // Scale histogram to MeV
507  // BCALHitMatrixD->GetZaxis()->SetTitle("Energy (MeV)");
508  // } else {
509  // BCALHitMatrixU->GetZaxis()->SetTitle("Energy (GeV)");
510  // }
511  // BCALHitMatrixD->Draw("colz");
512  for (int layer=0; layer<4; layer++) {
513  if (BCALPointZphiLayer[layer]->GetEntries()>0) {
514  BCALPointZphiLayer[layer]->SetMaximum(clustermax);
515  if (firstup==1) {
516  BCALPointZphiLayer[layer]->Draw("box");
517  firstup=0;
518  } else {
519  BCALPointZphiLayer[layer]->Draw("box,same");
520  }
521  }
522  }
523  LayerLegend->Draw();
524 
525  BCALHitCanvas->cd(3);
526  gPad->SetGridx();
527  gPad->SetGridy();
528  gPad->SetLeftMargin(leftmargin);
529  if (BCALClusterZphiHistos.size()>0) {
530  BCALClusterZphiHistos[0]->SetMaximum(clustermax);
531  BCALClusterZphiHistos[0]->Draw("box");
532  for (unsigned int k=1; k<BCALClusterZphiHistos.size(); k++) {
533  BCALClusterZphiHistos[k]->SetMaximum(clustermax);
534  BCALClusterZphiHistos[k]->Draw("box,same");
535  }
536  }
537  ClusterLegend->Draw();
538  // BCALParticles->Draw("colz");
539  // for (unsigned int n=0;n<BCALPLables.size();n++){
540  // BCALPLables[n]->Draw("same");
541  // }
542  BCALHitCanvas->Update();
543  }
544 
545 
546  // BCAL hits
547  if(hdvmf->GetCheckButton("bcal")){
548  vector<const DBCALHit*> bcalhits;
549  loop->Get(bcalhits);
550 
551  for(unsigned int i=0; i<bcalhits.size(); i++){
552  const DBCALHit *hit = bcalhits[i];
553  TPolyLine *poly = hdvmf->GetBCALPolyLine(hit->module, hit->layer, hit->sector);
554 
555  if(!poly)continue;
556 
557  // The aim is to have a log scale in energy
558  double E = 1000.0*hit->E; // Energy in MeV
559  double logE = log10(E);
560  // 3 = 1 GeV, 0 = 1 MeV, use range 0 through 4
561  // 0-1 White-Yellow
562  // 1-2 Yellow-Red
563  // 2-3 Red-Cyan
564  // 3-4 Cyan-Blue
565 
566  float r,g,b;
567  if (E<1){
568  r = 1.;
569  g = 1.;
570  b = 1.;
571  } else {
572  if (logE<1){
573  r = 1.;
574  g = 1.;
575  b = 1.-logE;
576  } else {
577  if (logE<2){
578  r = 1.;
579  g = 1.-(logE-1);
580  b = 0.;
581  } else {
582  if (logE<3){
583  r = 1.;
584  g = 0.;
585  b = 1.-(logE-2);
586  } else {
587  if (logE<4){
588  r = 1.-(logE-3);
589  g = 0.;
590  b = 1.;
591  } else {
592  r = 0;
593  g = 1.;
594  b = 0.;
595  //printf("High BCAL cell reconstructed energy: E=%.1f MeV\n",E);
596  }
597  }
598  }
599  }
600  }
601  if (r<0||g<0||b<0||r>1||g>1||b>1) printf("color error (r,g,b)=(%f,%f,%f)\n",r,g,b);
602 
603  poly->SetFillColor(TColor::GetColor(r,g,b));
604  poly->SetLineColor(TColor::GetColor(r,g,b));
605  poly->SetLineWidth(1);
606  poly->SetFillStyle(1001);
607  }
608  }
609 
610  // FCAL hits
611  if(hdvmf->GetCheckButton("fcal")){
612  vector<const DFCALHit*> fcalhits;
613  loop->Get(fcalhits);
614 
615  for(unsigned int i=0; i<fcalhits.size(); i++){
616  const DFCALHit *hit = fcalhits[i];
617  TPolyLine *poly = hdvmf->GetFCALPolyLine(hit->x, hit->y);
618  if(!poly)continue;
619 
620 #if 0
621  double a = hit->E/0.005;
622  double f = sqrt(a>1.0 ? 1.0:a<0.0 ? 0.0:a);
623  double grey = 0.8;
624  double s = 1.0 - f;
625 
626  float r = s*grey;
627  float g = s*grey;
628  float b = f*(1.0-grey) + grey;
629 #endif
630 
631  // The aim is to have a log scale in energy (see BCAL)
632  double E = 1000*hit->E; // Change Energy to MeV
633  if(E<0.0) continue;
634  double logE = log10(E);
635 
636  float r,g,b;
637  if (logE<0){
638  r = 1.;
639  g = 1.;
640  b = 1.;
641  } else {
642  if (logE<1){
643  r = 1.;
644  g = 1.;
645  b = 1.-logE;
646  } else {
647  if (logE<2){
648  r = 1.;
649  g = 1.-(logE-1);
650  b = 0.;
651  } else {
652  if (logE<3){
653  r = 1.;
654  g = 0.;
655  b = 1.-(logE-2);
656  } else {
657  if (logE<4){
658  r = 1.-(logE-3);
659  g = 0.;
660  b = 1.;
661  } else {
662  r = 0;
663  g = 0;
664  b = 0;
665  }
666  }
667  }
668  }
669  }
670  poly->SetFillColor(TColor::GetColor(r,g,b));
671  }
672  }
673 
674  // CCAL hits
675  if(hdvmf->GetCheckButton("ccal")){
676  vector<const DCCALHit*> ccalhits;
677  loop->Get(ccalhits);
678 
679  for(unsigned int i=0; i<ccalhits.size(); i++){
680  const DCCALHit *hit = ccalhits[i];
681  TPolyLine *poly = hdvmf->GetCCALPolyLine(hit->row, hit->column);
682  if(!poly)continue;
683 
684  // The aim is to have a log scale in energy (see BCAL)
685  double E = 1000*hit->E; // Change Energy to MeV
686  E /= 1000.0; // CCAL is currently not calibrated and appears to be at least 1E3 too large 2018-12-10 DL
687  if(E<0.0) continue;
688  double logE = log10(E);
689 
690  float r,g,b;
691  if (logE<0){
692  r = 1.;
693  g = 1.;
694  b = 1.;
695  } else {
696  if (logE<1){
697  r = 1.;
698  g = 1.;
699  b = 1.-logE;
700  } else {
701  if (logE<2){
702  r = 1.;
703  g = 1.-(logE-1);
704  b = 0.;
705  } else {
706  if (logE<3){
707  r = 1.;
708  g = 0.;
709  b = 1.-(logE-2);
710  } else {
711  if (logE<4){
712  r = 1.-(logE-3);
713  g = 0.;
714  b = 1.;
715  } else {
716  r = 0;
717  g = 0;
718  b = 0;
719  }
720  }
721  }
722  }
723  }
724  poly->SetFillColor(TColor::GetColor(r,g,b));
725  }
726  }
727 
728  // TOF hits
729  if(hdvmf->GetCheckButton("tof")){
730 
731  double hit_north[45];
732  double hit_south[45];
733  double hit_up[45];
734  double hit_down[45];
735 
736  memset(hit_north,0,sizeof(hit_north));
737  memset(hit_south,0,sizeof(hit_south));
738  memset(hit_up,0,sizeof(hit_up));
739  memset(hit_down,0,sizeof(hit_down));
740 
741  vector<const DTOFHit*> tofhits;
742  loop->Get(tofhits);
743 
744  for(unsigned int i=0; i<tofhits.size(); i++){
745  const DTOFHit *tof_hit = tofhits[i];
746 
747  int plane = tof_hit->plane;
748  int bar = tof_hit->bar;
749  int end = tof_hit->end;
750  float t = tof_hit->t;
751 
752  int translate_side;
753  TPolyLine *pmtPline;
754 
755 
756  // Flash the PMTs that do have fADC hits
757 
758  /*
759  double dE_padd = 0.2/5.2E5 * 40000;
760  double thold = 0.2/5.2E5 * 115;
761  if (tof_hit->has_fADC && (tof_hit->dE - dE_padd > thold)){
762  */
763 
764  if (tof_hit->has_fADC){
765  switch(plane)
766  {
767  case 0:
768  if(end == 1){
769  //cout << "Down : " << bar << endl;
770  translate_side = 0;
771  pmtPline = hdvmf->GetTOFPolyLine(translate_side, bar);
772  pmtPline->SetFillColor(2);
773  }
774  else if(end == 0){
775  //cout << "Up : " << bar << endl;
776  translate_side = 2;
777  pmtPline = hdvmf->GetTOFPolyLine(translate_side, bar);
778  pmtPline->SetFillColor(2);
779  }
780  else{
781  cerr << "Out of range TOF end" << endl;
782  }
783  break;
784  case 1:
785  if(end == 0){
786  //cout << "North : " << bar << endl;
787  translate_side = 3;
788  pmtPline = hdvmf->GetTOFPolyLine(translate_side, bar);
789  pmtPline->SetFillColor(2);
790  }
791  else if(end == 1){
792  //cout << "South : " << bar << endl;
793  translate_side = 1;
794  pmtPline = hdvmf->GetTOFPolyLine(translate_side, bar);
795  pmtPline->SetFillColor(2);
796  }
797  else{
798  cerr << "Out of range TOF end" << endl;
799  }
800  break;
801  default:
802  cerr << "Out of range TOF plane" << endl;
803  exit(0);
804  } // close the switch plane loop
805  } // if for the fADC
806 
807  // Draw the position from the events that do have tdc hits
808  // with the current status of the TOFHit object those hits appear with no match from the fADC
809  //Float_t hit_dist;
810 
811  if (tof_hit->has_TDC){
812  switch(plane)
813  {
814  case 0:
815  if(end == 1){
816  if (hit_down[bar]<=0 || (t < hit_down[bar]) ){
817  hit_down[bar] = t;
818  }
819  }
820  else if(end == 0){
821  if (hit_up[bar]<=0 || (t < hit_up[bar]) ){
822  hit_up[bar] = t;
823  }
824  }
825  else{
826  cerr << "Out of range TOF end" << endl;
827  }
828  break;
829  case 1:
830  if(end == 0){
831  if (hit_north[bar]<=0 || (t < hit_north[bar]) ){
832  hit_north[bar] = t;
833  }
834  }
835  else if(end == 1){
836  if (hit_south[bar]<=0 || (t < hit_south[bar]) ){
837  hit_south[bar] = t;
838  }
839  }
840  else{
841  cerr << "Out of range TOF end" << endl;
842  }
843  break;
844  default:
845  cerr << "Out of range TOF plane" << endl;
846  exit(0);
847  }
848 
849  } // close the switch if there is a TDC Hit
850 
851  } // close the for TOFHit object
852 
853  // Draw the TDC Points here
854 
855  Float_t hit_dist;
856  Float_t distY_Horz = -126; // Horizontal plane start counting from the Bottom to Top
857  Float_t distX_Vert = -126; // Vertical plane start counting from the South to North
858  int tdc_hits = 0;
859 
860  for(Int_t i_tdc = 1; i_tdc <= 44; i_tdc++){
861  if ( i_tdc == 20 || i_tdc == 21 || i_tdc == 24 || i_tdc == 25 ){
862  distY_Horz = distY_Horz + 1.5;
863  distX_Vert = distX_Vert + 1.5;
864  }
865  else{
866  distY_Horz = distY_Horz + 3.0;
867  distX_Vert = distX_Vert + 3.0;
868  }
869  if(hit_north[i_tdc] > 0 && hit_south[i_tdc] > 0){
870  hit_dist = (15.2*(Float_t(hit_south[i_tdc] - hit_north[i_tdc])/2));
871  TArc *tdc_cir = new TArc(hit_dist,distY_Horz,2);
872  tdc_cir->SetFillColor(kGreen);
873 
874  graphics_tof_hits.push_back(tdc_cir);
875  tdc_hits++;
876  }
877  if(hit_up[i_tdc] > 0 && hit_down[i_tdc] > 0){
878  hit_dist = (15.2*(Float_t(hit_down[i_tdc] - hit_up[i_tdc])/2) );
879  TArc *tdc_cir = new TArc(distX_Vert,hit_dist,2);
880  tdc_cir->SetFillColor(kBlue);
881 
882  graphics_tof_hits.push_back(tdc_cir);
883  tdc_hits++;
884  }
885  if ( i_tdc == 20 || i_tdc == 21 || i_tdc == 24 || i_tdc == 25 ){
886  distY_Horz = distY_Horz + 1.5;
887  distX_Vert = distX_Vert + 1.5;
888  }
889  else{
890  distY_Horz = distY_Horz + 3.0;
891  distX_Vert = distX_Vert + 3.0;
892  }
893  }
894 
895  } // close the if check button for the TOF
896 
897  // Start counter hits
898  for (unsigned int i=0;i<schits.size();i++){
899  DGraphicSet gset(6,kLine,2.0);
900  double r_start=7.7493;
901  double phi0=0.2094395*(schits[i]->sector-1); // span 12 deg in phi
902  double phi1=0.2094395*(schits[i]->sector);
903  TVector3 point1(r_start*cos(phi0),r_start*sin(phi0),38.75);
904  gset.points.push_back(point1); // upstream end of sctraight section of scint
905  TVector3 point2(r_start*cos(phi1),r_start*sin(phi1),38.75);
906  gset.points.push_back(point2);
907  TVector3 point3(r_start*cos(phi1),r_start*sin(phi1),78.215);
908  gset.points.push_back(point3); // downstream end of sctraight section of scint
909  TVector3 point4(r_start*cos(phi0),r_start*sin(phi0),78.215);
910  gset.points.push_back(point4);
911  TVector3 point5(r_start*cos(phi0),r_start*sin(phi0),38.75);
912  gset.points.push_back(point5);
913 
914  /*
915  // ST dimensions
916  Double_t dtr = 1.74532925e-02; // Conversion factor from degrees to radians
917  Double_t st_straight = 39.465; // Distance of the straight section along scintillator path
918  Double_t st_arc_angle = 18.5; // Angle of the bend arc
919  Double_t st_arc_radius = 12.0; // Radius of the bend arc
920  Double_t st_to_beam = 7.74926; // Distance from beam to bottom side of scintillator paddle
921  Double_t st_len_cone = 16.056; // Length of the cone section along scintillator path
922  Double_t st_thickness = 0.3; // Thickness of the scintillator paddles
923  Double_t st_beam_to_arc_radius = 4.25; // Distance from the beam line to the arc radius
924 
925  // Nose Arrays
926  Double_t tp_nose_z[5];
927  Double_t tp_nose_y[5];
928  Double_t bm_nose_z[5];
929  Double_t bm_nose_y[5];
930 
931  // Offsets for Hall coordinates
932  Double_t z_center = 65.0; // Target Center (0,0,65)
933  Double_t us_end_pt = -26.25; // Distance of upstream end relative to target
934  Double_t ds_end_pt = 97.4; // Distance of downstream end relative
935 
936  // Top start counter paddle straight section
937  TPolyLine *top_paddle_straight;
938  Double_t tp_straight_z[5] = {us_end_pt + z_center, us_end_pt + z_center, us_end_pt + st_straight + z_center, us_end_pt + st_straight + z_center, us_end_pt + z_center};
939  Double_t tp_straight_y[5] = {st_to_beam, st_to_beam + st_thickness, st_to_beam + st_thickness, st_to_beam, st_to_beam};
940  */
941  graphics.push_back(gset);
942  }
943 
944  // CDC hits
945  if(hdvmf->GetCheckButton("cdc")){
946  vector<const DCDCTrackHit*> cdctrackhits;
947  loop->Get(cdctrackhits);
948 
949  for(unsigned int i=0; i<cdctrackhits.size(); i++){
950  const DCDCWire *wire = cdctrackhits[i]->wire;
951 
952  int color = (cdctrackhits[i]->tdrift>-20 && cdctrackhits[i]->tdrift<800) ? kCyan:kYellow;
953  DGraphicSet gset(color, kLine, 1.0);
954  DVector3 dpoint=wire->origin-(wire->L/2.0)*wire->udir;
955  TVector3 tpoint(dpoint.X(),dpoint.Y(),dpoint.Z());
956  gset.points.push_back(tpoint);
957  dpoint=wire->origin+(wire->L/2.0)*wire->udir;
958  tpoint.SetXYZ(dpoint.X(),dpoint.Y(),dpoint.Z());
959  gset.points.push_back(tpoint);
960  graphics.push_back(gset);
961 
962  // Rings for drift times.
963  // NOTE: These are not perfect since they still have TOF in them
964  if(hdvmf->GetCheckButton("cdcdrift") && fabs(wire->stereo)<0.05){
965  double x = wire->origin.X();
966  double y = wire->origin.Y();
967  double dist1 = cdctrackhits[i]->dist;
968  TEllipse *e = new TEllipse(x, y, dist1, dist1);
969  e->SetLineColor(38);
970  e->SetFillStyle(0);
971  graphics_xyA.push_back(e);
972 
973  double dist2 = dist1 - 4.0*55.0E-4; // what if our TOF was 4ns?
974  e = new TEllipse(x, y, dist2, dist2);
975  e->SetLineColor(38);
976  e->SetFillStyle(0);
977  graphics_xyA.push_back(e);
978  }
979  }
980  }
981 
982  // FDC wire
983  if(hdvmf->GetCheckButton("fdcwire")){
984  vector<const DFDCHit*> fdchits;
985  loop->Get(fdchits);
986 
987  for(unsigned int i=0; i<fdchits.size(); i++){
988  const DFDCHit *fdchit = fdchits[i];
989  if(fdchit->type!=0)continue;
990  const DFDCWire *wire =fdcwires[fdchit->gLayer-1][fdchit->element-1];
991  if(!wire){
992  _DBG_<<"Couldn't find wire for gLayer="<<fdchit->gLayer<<" and element="<<fdchit->element<<endl;
993  continue;
994  }
995 
996  // Wire
997  int color = (fdchit->t>-50 && fdchit->t<2000) ? kCyan:kYellow;
998  DGraphicSet gset(color, kLine, 1.0);
999  DVector3 dpoint=wire->origin-(wire->L/2.0)*wire->udir;
1000  TVector3 tpoint(dpoint.X(),dpoint.Y(),dpoint.Z());
1001  gset.points.push_back(tpoint);
1002  dpoint=wire->origin+(wire->L/2.0)*wire->udir;
1003  tpoint.SetXYZ(dpoint.X(),dpoint.Y(),dpoint.Z());
1004  gset.points.push_back(tpoint);
1005  graphics.push_back(gset);
1006  }
1007  }
1008 
1009  // FDC intersection hits
1010  if(hdvmf->GetCheckButton("fdcintersection")){
1011  vector<const DFDCIntersection*> fdcints;
1012  loop->Get(fdcints);
1013  DGraphicSet gsetp(46, kMarker, 0.5);
1014 
1015  for(unsigned int i=0; i<fdcints.size(); i++){
1016  TVector3 tpos(fdcints[i]->pos.X(),fdcints[i]->pos.Y(),
1017  fdcints[i]->pos.Z());
1018  gsetp.points.push_back(tpos);
1019  }
1020  graphics.push_back(gsetp);
1021  }
1022 
1023  // FDC psuedo hits
1024  if(hdvmf->GetCheckButton("fdcpseudo")){
1025  vector<const DFDCPseudo*> fdcpseudos;
1026  loop->Get(fdcpseudos);
1027  DGraphicSet gsetp(38, kMarker, 0.5);
1028 
1029  for(unsigned int i=0; i<fdcpseudos.size(); i++){
1030  const DFDCWire *wire = fdcpseudos[i]->wire;
1031 
1032  // Pseudo point
1033  TVector3 pos(fdcpseudos[i]->xy.X(),
1034  fdcpseudos[i]->xy.Y(), wire->origin.Z());
1035  gsetp.points.push_back(pos);
1036  }
1037  graphics.push_back(gsetp);
1038  }
1039 
1040  // DMCThrown
1041  if(hdvmf->GetCheckButton("thrown")){
1042  vector<const DMCThrown*> mcthrown;
1043  loop->Get(mcthrown);
1044  for(unsigned int i=0; i<mcthrown.size(); i++){
1045  int color=14;
1046  double size=1.5;
1047  if(mcthrown[i]->charge()==0.0) color = kGreen;
1048  if(mcthrown[i]->charge() >0.0) color = kBlue;
1049  if(mcthrown[i]->charge() <0.0) color = kRed;
1050  switch(mcthrown[i]->type){
1051  case Gamma:
1052  case Positron:
1053  case Electron:
1054  size = 1.0;
1055  break;
1056  case Pi0:
1057  case PiPlus:
1058  case PiMinus:
1059  size = 2.0;
1060  break;
1061  case Neutron:
1062  case Proton:
1063  case AntiProton:
1064  size = 3.0;
1065  break;
1066  }
1067  AddKinematicDataTrack(mcthrown[i], color, size);
1068  }
1069  }
1070 
1071  // CDC Truth points
1072  if(hdvmf->GetCheckButton("cdctruth")){
1073  vector<const DMCTrackHit*> mctrackhits;
1074  loop->Get(mctrackhits);
1075  DGraphicSet gset(12, kMarker, 0.5);
1076  for(unsigned int i=0; i<mctrackhits.size(); i++){
1077  const DMCTrackHit *hit = mctrackhits[i];
1078  if(hit->system != SYS_CDC)continue;
1079 
1080  TVector3 pos(hit->r*cos(hit->phi), hit->r*sin(hit->phi), hit->z);
1081  gset.points.push_back(pos);
1082  }
1083  graphics.push_back(gset);
1084  }
1085 
1086  // FDC Truth points
1087  if(hdvmf->GetCheckButton("fdctruth")){
1088  vector<const DMCTrackHit*> mctrackhits;
1089  loop->Get(mctrackhits);
1090  DGraphicSet gset(12, kMarker, 0.5);
1091  for(unsigned int i=0; i<mctrackhits.size(); i++){
1092  const DMCTrackHit *hit = mctrackhits[i];
1093  if(hit->system != SYS_FDC)continue;
1094 
1095  TVector3 pos(hit->r*cos(hit->phi), hit->r*sin(hit->phi), hit->z);
1096  gset.points.push_back(pos);
1097  }
1098  graphics.push_back(gset);
1099  }
1100 
1101  // Track Hits for Track Candidates and Candidate trajectory in Debuger Window
1102  for(unsigned int n=0; n<trCand.size(); n++){
1103  if (n>9)
1104  break;
1105  char str1[128];
1106  sprintf(str1,"Candidate%d",n+1);
1107 
1108  if(hdvmf->GetCheckButton(str1)){
1109 
1110  int color = n+1;
1111  if (color > 4)
1112  color++;
1113  if (color > 6)
1114  color++;
1115 
1116  AddKinematicDataTrack(trCand[n], color, 1.5);
1117 
1118  vector<const DCDCTrackHit*> cdctrackhits;
1119  trCand[n]->Get(cdctrackhits,"",1);
1120  for(unsigned int i=0; i<cdctrackhits.size(); i++){
1121  const DCDCWire *wire = cdctrackhits[i]->wire;
1122  DGraphicSet gset(color, kLine, 1.0);
1123  DVector3 dpoint=wire->origin-(wire->L/2.0)*wire->udir;
1124  TVector3 tpoint(dpoint.X(),dpoint.Y(),dpoint.Z());
1125  gset.points.push_back(tpoint);
1126  dpoint=wire->origin+(wire->L/2.0)*wire->udir;
1127  tpoint.SetXYZ(dpoint.X(),dpoint.Y(),dpoint.Z());
1128  gset.points.push_back(tpoint);
1129  graphics.push_back(gset);
1130 
1131  } // end loop of cdc hits of track candidate
1132  vector<const DFDCPseudo*> fdcpseudos;
1133  trCand[n]->Get(fdcpseudos,"",1);
1134  DGraphicSet gsetp(color, kMarker, 0.5);
1135 
1136  for(unsigned int i=0; i<fdcpseudos.size(); i++){
1137  const DFDCWire *wire = fdcpseudos[i]->wire;
1138 
1139  // Pseudo point
1140  TVector3 pos(fdcpseudos[i]->xy.X(),
1141  fdcpseudos[i]->xy.Y(), wire->origin.Z());
1142  gsetp.points.push_back(pos);
1143  }
1144  graphics.push_back(gsetp);
1145 
1146  }
1147  }
1148 
1149  // Wire Based Track Hits and trajectory for Debuger Window
1150  for(unsigned int n=0; n<trWB.size(); n++){
1151  if (n>=MaxWireTracks)
1152  break;
1153  char str1[128];
1154  sprintf(str1,"WireBased%d",n+1);
1155 
1156  if(hdvmf->GetCheckButton(str1)){
1157 
1158  int color = trWB[n]->candidateid;
1159  if (color > 4)
1160  color++;
1161  if (color > 6)
1162  color++;
1163 
1164  AddKinematicDataTrack(trWB[n], color, 1.5);
1165 
1166  vector<const DCDCTrackHit*> cdctrackhits;
1167  trWB[n]->Get(cdctrackhits,"",1);
1168  for(unsigned int i=0; i<cdctrackhits.size(); i++){
1169  const DCDCWire *wire = cdctrackhits[i]->wire;
1170  DGraphicSet gset(color, kLine, 1.0);
1171  DVector3 dpoint=wire->origin-(wire->L/2.0)*wire->udir;
1172  TVector3 tpoint(dpoint.X(),dpoint.Y(),dpoint.Z());
1173  gset.points.push_back(tpoint);
1174  dpoint=wire->origin+(wire->L/2.0)*wire->udir;
1175  tpoint.SetXYZ(dpoint.X(),dpoint.Y(),dpoint.Z());
1176  gset.points.push_back(tpoint);
1177  graphics.push_back(gset);
1178 
1179  } // end loop of cdc hits of track candidate
1180  vector<const DFDCPseudo*> fdcpseudos;
1181  trWB[n]->Get(fdcpseudos,"",1);
1182  DGraphicSet gsetp(color, kMarker, 0.5);
1183 
1184  for(unsigned int i=0; i<fdcpseudos.size(); i++){
1185  const DFDCWire *wire = fdcpseudos[i]->wire;
1186 
1187  // Pseudo point
1188  TVector3 pos(fdcpseudos[i]->xy.X(),
1189  fdcpseudos[i]->xy.Y(), wire->origin.Z());
1190  gsetp.points.push_back(pos);
1191  }
1192  graphics.push_back(gsetp);
1193 
1194  }
1195  }
1196 
1197  // Time Based Track Hits and trajectory for Debuger Window
1198  for(unsigned int n=0; n<trTB.size(); n++){
1199  if (n>=MaxTimeTracks)
1200  break;
1201  char str1[128];
1202  sprintf(str1,"TimeBased%d",n+1);
1203 
1204  if(hdvmf->GetCheckButton(str1)){
1205 
1206  int color = trTB[n]->candidateid;
1207  if (color > 4)
1208  color++;
1209  if (color > 6)
1210  color++;
1211 
1212  AddKinematicDataTrack(trTB[n], color, 1.5);
1213 
1214  vector<const DCDCTrackHit*> cdctrackhits;
1215  trTB[n]->Get(cdctrackhits,"",1);
1216  for(unsigned int i=0; i<cdctrackhits.size(); i++){
1217  const DCDCWire *wire = cdctrackhits[i]->wire;
1218  DGraphicSet gset(color, kLine, 1.0);
1219  DVector3 dpoint=wire->origin-(wire->L/2.0)*wire->udir;
1220  TVector3 tpoint(dpoint.X(),dpoint.Y(),dpoint.Z());
1221  gset.points.push_back(tpoint);
1222  dpoint=wire->origin+(wire->L/2.0)*wire->udir;
1223  tpoint.SetXYZ(dpoint.X(),dpoint.Y(),dpoint.Z());
1224  gset.points.push_back(tpoint);
1225  graphics.push_back(gset);
1226 
1227  } // end loop of cdc hits of track candidate
1228  vector<const DFDCPseudo*> fdcpseudos;
1229  trTB[n]->Get(fdcpseudos,"",1);
1230  DGraphicSet gsetp(color, kMarker, 0.5);
1231 
1232  for(unsigned int i=0; i<fdcpseudos.size(); i++){
1233  const DFDCWire *wire = fdcpseudos[i]->wire;
1234 
1235  // Pseudo point
1236  TVector3 pos(fdcpseudos[i]->xy.X(),
1237  fdcpseudos[i]->xy.Y(), wire->origin.Z());
1238  gsetp.points.push_back(pos);
1239  }
1240  graphics.push_back(gsetp);
1241 
1242  }
1243  }
1244 
1245  // TOF reconstructed points
1246  if (hdvmf->GetCheckButton("tof")){
1247  vector<const DTOFPoint *>tofpoints;
1248  loop->Get(tofpoints);
1249  DGraphicSet gset(kRed, kMarker, 0.5);
1250  for(unsigned int i=0; i<tofpoints.size(); i++){
1251  const DTOFPoint *hit = tofpoints[i];
1252  TVector3 pos(hit->pos.x(),hit->pos.y(),hit->pos.z());
1253  gset.points.push_back(pos);
1254  }
1255  graphics.push_back(gset);
1256  }
1257 
1258  // TOF Truth points
1259  if(hdvmf->GetCheckButton("toftruth")){
1260  vector<const DMCTrackHit*> mctrackhits;
1261  loop->Get(mctrackhits);
1262  DGraphicSet gset(kBlack, kMarker, 0.5);
1263  for(unsigned int i=0; i<mctrackhits.size(); i++){
1264  const DMCTrackHit *hit = mctrackhits[i];
1265  if(hit->system != SYS_TOF)continue;
1266 
1267  TVector3 pos(hit->r*cos(hit->phi), hit->r*sin(hit->phi), hit->z);
1268  gset.points.push_back(pos);
1269  }
1270  graphics.push_back(gset);
1271  }
1272 
1273  // BCAL Truth points
1274  if(hdvmf->GetCheckButton("bcaltruth")){
1275  vector<const DMCTrackHit*> mctrackhits;
1276  loop->Get(mctrackhits);
1277  DGraphicSet gset(kBlack, kMarker, 1.0);
1278  for(unsigned int i=0; i<mctrackhits.size(); i++){
1279  const DMCTrackHit *hit = mctrackhits[i];
1280  if(hit->system != SYS_BCAL)continue;
1281 
1282  TVector3 pos(hit->r*cos(hit->phi), hit->r*sin(hit->phi), hit->z);
1283  gset.points.push_back(pos);
1284 
1285  TMarker *m = new TMarker(pos.X(), pos.Y(), 2);
1286  graphics_xyA.push_back(m);
1287  }
1288  //graphics.push_back(gset);
1289  }
1290 
1291  // FCAL Truth points
1292  if(hdvmf->GetCheckButton("fcaltruth")){
1293  vector<const DFCALGeometry*> fcalgeometries;
1294  vector<const DFCALHit*> mcfcalhits;
1295  loop->Get(fcalgeometries);
1296  loop->Get(mcfcalhits, "TRUTH");
1297  if(fcalgeometries.size()>0){
1298  const DFCALGeometry *fgeom = fcalgeometries[0];
1299 
1300  DGraphicSet gset(kBlack, kMarker, 0.25);
1301  for(unsigned int i=0; i<mcfcalhits.size(); i++){
1302  const DFCALHit *hit = mcfcalhits[i];
1303 
1304  DVector2 pos_face = fgeom->positionOnFace(hit->row, hit->column);
1305  TVector3 pos(pos_face.X(), pos_face.Y(), FCAL_Zmin);
1306  gset.points.push_back(pos);
1307 
1308  TMarker *m = new TMarker(pos.X(), pos.Y(), 2);
1309  //m->SetColor(kGreen);
1310  //m->SetLineWidth(1);
1311  graphics_xyB.push_back(m);
1312 
1313  TMarker *m1 = new TMarker(pos.Z(), pos.X(), 2);
1314  graphics_xz.push_back(m1);
1315  TMarker *m2 = new TMarker(pos.Z(), pos.Y(), 2);
1316  graphics_yz.push_back(m2);
1317  }
1318  //graphics.push_back(gset);
1319  }
1320  }
1321 
1322  // BCAL reconstructed photons
1323  if(hdvmf->GetCheckButton("recon_photons_bcal")){
1324  vector<const DNeutralParticle*> neutrals;
1325  loop->Get(neutrals);
1326 
1327  DGraphicSet gset(kYellow+2, kMarker, 1.25);
1328  gset.marker_style=21;
1329  for(unsigned int i=0; i<neutrals.size(); i++){
1330  auto locNeutralShower = neutrals[i]->dNeutralShower;
1331  DetectorSystem_t locDetectorSystem = locNeutralShower->dDetectorSystem;
1332  if(locDetectorSystem == SYS_BCAL){
1333  TVector3 pos( locNeutralShower->dSpacetimeVertex.X(),
1334  locNeutralShower->dSpacetimeVertex.Y(),
1335  locNeutralShower->dSpacetimeVertex.Z());
1336  gset.points.push_back(pos);
1337 
1338  double dist2 = 2.0 + 5.0*locNeutralShower->dEnergy;
1339  TEllipse *e = new TEllipse(pos.X(), pos.Y(), dist2, dist2);
1340  e->SetLineColor(kGreen);
1341  e->SetFillStyle(0);
1342  e->SetLineWidth(2);
1343  graphics_xyA.push_back(e);
1344  }
1345  }
1346  //graphics.push_back(gset);
1347  }
1348 
1349  // FCAL reconstructed photons
1350  if(hdvmf->GetCheckButton("recon_photons_fcal")){
1351  vector<const DNeutralParticle*> neutrals;
1352  loop->Get(neutrals);
1353  DGraphicSet gset(kOrange, kMarker, 1.25);
1354  gset.marker_style=2;
1355  for(unsigned int i=0; i<neutrals.size(); i++){
1356  auto locNeutralShower = neutrals[i]->dNeutralShower;
1357  DetectorSystem_t locDetectorSystem = locNeutralShower->dDetectorSystem;
1358  if(locDetectorSystem == SYS_FCAL){
1359 
1360  TVector3 pos( locNeutralShower->dSpacetimeVertex.X(),
1361  locNeutralShower->dSpacetimeVertex.Y(),
1362  locNeutralShower->dSpacetimeVertex.Z());
1363  gset.points.push_back(pos);
1364 
1365  double dist2 = 2.0 + 10.0*locNeutralShower->dEnergy;
1366  TEllipse *e = new TEllipse(pos.X(), pos.Y(), dist2, dist2);
1367  e->SetLineColor(kGreen);
1368  e->SetFillStyle(0);
1369  e->SetLineWidth(2);
1370  graphics_xyB.push_back(e);
1371 
1372  TEllipse *e1 = new TEllipse(pos.Z(), pos.X(), dist2, dist2);
1373  e1->SetLineColor(kGreen);
1374  e1->SetFillStyle(0);
1375  e1->SetLineWidth(2);
1376  graphics_xz.push_back(e1);
1377  TEllipse *e2 = new TEllipse(pos.Z(), pos.Y(), dist2, dist2);
1378  e2->SetLineColor(kGreen);
1379  e2->SetFillStyle(0);
1380  e2->SetLineWidth(2);
1381  graphics_yz.push_back(e2);
1382 
1383  }
1384  }
1385  //graphics.push_back(gset);
1386  }
1387 
1388  // Reconstructed photons matched with tracks
1389  if(hdvmf->GetCheckButton("recon_photons_track_match")){
1390  vector<const DChargedTrack*> ctracks;
1391  loop->Get(ctracks);
1392  for(unsigned int i=0; i<ctracks.size(); i++){
1393  const DChargedTrack *locCTrack = ctracks[i];
1394  vector<const DNeutralShower*> locNeutralShowers;
1395  locCTrack->GetT(locNeutralShowers);
1396 
1397  if (!locNeutralShowers.size()) continue;
1398  const DNeutralShower *locNeutralShower = locNeutralShowers[0];
1399 
1400 
1401  // Decide if this hit BCAL of FCAL based on z of position on calorimeter
1402  bool is_bcal = (locNeutralShower->dDetectorSystem == SYS_BCAL );
1403 
1404  // Draw on all frames except FCAL frame
1405  DGraphicSet gset(kRed, kMarker, 1.25);
1406  gset.marker_style = is_bcal ? 22:3;
1407  TVector3 tpos( locNeutralShower->dSpacetimeVertex.X(),
1408  locNeutralShower->dSpacetimeVertex.Y(),
1409  locNeutralShower->dSpacetimeVertex.Z());
1410  gset.points.push_back(tpos);
1411  graphics.push_back(gset);
1412 
1413  // For BCAL hits, don't draw them on FCAL pane
1414  if(is_bcal)continue;
1415 
1416  // Draw on FCAL pane
1417  double dist2 = 2.0 + 2.0*locNeutralShower->dEnergy;
1418  TEllipse *e = new TEllipse(tpos.X(), tpos.Y(), dist2, dist2);
1419  e->SetLineColor(gset.color);
1420  e->SetFillStyle(0);
1421  e->SetLineWidth(1);
1422  TMarker *m = new TMarker(tpos.X(), tpos.Y(), gset.marker_style);
1423  m->SetMarkerColor(gset.color);
1424  m->SetMarkerSize(1.75);
1425  graphics_xyB.push_back(e);
1426  graphics_xyB.push_back(m);
1427  }
1428  }
1429 
1430  // FCAL and BCAL thrown photon projections
1431  if(hdvmf->GetCheckButton("thrown_photons_fcal") || hdvmf->GetCheckButton("thrown_photons_bcal")){
1432  vector<const DMCThrown*> throwns;
1433  loop->Get(throwns);
1434  DGraphicSet gset(kSpring, kMarker, 1.25);
1435  for(unsigned int i=0; i<throwns.size(); i++){
1436  const DMCThrown *thrown = throwns[i];
1437  if(thrown->charge()!=0.0)continue;
1438 
1439  // This may seem a little funny, but it saves swimming the reference trajectory
1440  // multiple times. The GetIntersectionWithCalorimeter() method will find the
1441  // intersection point of the photon with whichever calorimeter it actually hits
1442  // or if it doesn't hit either of them. Then, we decide to draw the marker based
1443  // on whether the flag is set to draw the calorimeter it hit or not.
1444  DVector3 pos;
1445  DetectorSystem_t who;
1446  GetIntersectionWithCalorimeter(thrown, pos, who);
1447 
1448  if(who!=SYS_FCAL && who!=SYS_BCAL)continue;
1449  if(who==SYS_FCAL && !hdvmf->GetCheckButton("thrown_photons_fcal"))continue;
1450  if(who==SYS_BCAL && !hdvmf->GetCheckButton("thrown_photons_bcal"))continue;
1451  TVector3 tpos(pos.X(),pos.Y(),pos.Z());
1452  gset.points.push_back(tpos);
1453 
1454  // Only draw on FCAL pane if photon hits FCAL
1455  if(who==SYS_BCAL)continue;
1456 
1457  double dist2 = 2.0 + 2.0*thrown->energy();
1458  TEllipse *e = new TEllipse(pos.X(), pos.Y(), dist2, dist2);
1459  e->SetLineColor(kSpring);
1460  e->SetFillStyle(0);
1461  e->SetLineWidth(4);
1462  graphics_xyB.push_back(e);
1463  }
1464  graphics.push_back(gset);
1465  }
1466 
1467  // FCAL and BCAL thrown charged particle projections
1468  if(hdvmf->GetCheckButton("thrown_charged_fcal") || hdvmf->GetCheckButton("thrown_charged_bcal")){
1469  vector<const DMCThrown*> throwns;
1470  loop->Get(throwns);
1471 
1472  for(unsigned int i=0; i<throwns.size(); i++){
1473  const DMCThrown *thrown = throwns[i];
1474  if(thrown->charge()==0.0)continue;
1475 
1476  // This may seem a little funny, but it saves swimming the reference trajectory
1477  // multiple times. The GetIntersectionWithCalorimeter() method will find the
1478  // intersection point of the photon with whichever calorimeter it actually hits
1479  // or if it doesn't hit either of them. Then, we decide to draw the marker based
1480  // on whether the flag is set to draw the calorimeter it hit or not.
1481  DVector3 pos;
1482  DetectorSystem_t who;
1483  GetIntersectionWithCalorimeter(thrown, pos, who);
1484 
1485  if(who!=SYS_FCAL && who!=SYS_BCAL)continue;
1486  if(who==SYS_FCAL && !hdvmf->GetCheckButton("thrown_charged_fcal"))continue;
1487  if(who==SYS_BCAL && !hdvmf->GetCheckButton("thrown_charged_bcal"))continue;
1488 
1489  DGraphicSet gset(thrown->charge()>0.0 ? kBlue:kRed, kMarker, 1.25);
1490  TVector3 tpos(pos.X(),pos.Y(),pos.Z());
1491  gset.points.push_back(tpos);
1492  graphics.push_back(gset);
1493 
1494  //double dist2 = 6.0 + 2.0*thrown->momentum().Mag();
1495  double dist2 = DELTA_R_FCAL;
1496  TEllipse *e = new TEllipse(pos.X(), pos.Y(), dist2, dist2);
1497  e->SetLineColor(thrown->charge()>0.0 ? kBlue:kRed);
1498  e->SetFillStyle(0);
1499  e->SetLineWidth(4);
1500  graphics_xyB.push_back(e);
1501  }
1502  }
1503 
1504  // FCAL and BCAL reconstructed charged particle projections
1505  if(hdvmf->GetCheckButton("recon_charged_bcal") || hdvmf->GetCheckButton("recon_charged_fcal")){
1506 
1507  // Here we used the full time-based track list, even though it includes multiple
1508  // hypotheses for each candidate. This is because currently, the photon/track
1509  // matching code in PID/DPhoton_factory.cc uses the DTrackTimeBased objects and
1510  // the current purpose of drawing these is to see matching of reconstructed
1511  // charged tracks with calorimeter clusters.
1512  vector<const DTrackTimeBased*> tracks;
1513  loop->Get(tracks, hdvmf->GetFactoryTag("DTrackTimeBased"));
1514 
1515  for(unsigned int i=0; i<tracks.size(); i++){
1516  const DTrackTimeBased *track = tracks[i];
1517 
1518  // See notes for above sections
1519  DVector3 pos;
1520  DetectorSystem_t who;
1521  GetIntersectionWithCalorimeter(track, pos, who);
1522 
1523  if(who!=SYS_FCAL && who!=SYS_BCAL)continue;
1524  if(who==SYS_FCAL && !hdvmf->GetCheckButton("recon_charged_fcal"))continue;
1525  if(who==SYS_BCAL && !hdvmf->GetCheckButton("recon_charged_bcal"))continue;
1526 
1527  DGraphicSet gset(track->charge()>0.0 ? kBlue:kRed, kMarker, 1.25);
1528  TVector3 tpos(pos.X(),pos.Y(),pos.Z());
1529  gset.points.push_back(tpos);
1530  graphics.push_back(gset);
1531 
1532  if(who==SYS_BCAL)continue; // Don't draw tracks hitting BCAL on FCAL pane
1533 
1534  //double dist2 = 6.0 + 2.0*track->momentum().Mag();
1535  double dist2 = DELTA_R_FCAL;
1536  TEllipse *e = new TEllipse(pos.X(), pos.Y(), dist2, dist2);
1537  e->SetLineColor(track->charge()>0.0 ? kBlue:kRed);
1538  e->SetFillStyle(0);
1539  e->SetLineWidth(4);
1540  graphics_xyB.push_back(e);
1541  }
1542  }
1543 
1544  // CCAL reconstructed clusters
1545  if(hdvmf->GetCheckButton("recon_photons_ccal")){
1546  vector<const DCCALShower*> clusters;
1547  loop->Get(clusters);
1548  for(auto cluster : clusters){
1549 
1550  double E = cluster->E/1000.0; // divide by 1000 since energy does not seem to be calibrated to GeV at the moment. 2018-12-10 DL
1551  double dist2 = 1.0 + 0.5*E;
1552  DVector3 pos(cluster->x,cluster->y,cluster->z);
1553 
1554  TEllipse *e = new TEllipse(pos.X(), pos.Y(), dist2, dist2);
1555  e->SetLineColor(kGreen);
1556  e->SetFillStyle(0);
1557  e->SetLineWidth(2);
1558  graphics_xyB.push_back(e);
1559 
1560  TEllipse *e1 = new TEllipse(pos.Z(), pos.X(), dist2, dist2);
1561  e1->SetLineColor(kGreen);
1562  e1->SetFillStyle(0);
1563  e1->SetLineWidth(2);
1564 
1565  graphics_xz.push_back(e1);
1566  TEllipse *e2 = new TEllipse(pos.Z(), pos.Y(), dist2, dist2);
1567  e2->SetLineColor(kGreen);
1568  e2->SetFillStyle(0);
1569  e2->SetLineWidth(2);
1570  graphics_yz.push_back(e2);
1571 
1572  }
1573  //graphics.push_back(gset);
1574  }
1575 
1576  // DMCTrajectoryPoints
1577  if(hdvmf->GetCheckButton("trajectories")){
1578  vector<const DMCTrajectoryPoint*> mctrajectorypoints;
1579  loop->Get(mctrajectorypoints);
1580  //sort(mctrajectorypoints.begin(), mctrajectorypoints.end(), DMCTrajectoryPoint_track_cmp);
1581 
1582  poly_type drawtype = hdvmf->GetCheckButton("trajectories_lines") ? kLine:kMarker;
1583  double drawsize = hdvmf->GetCheckButton("trajectories_lines") ? 1.0:0.3;
1584  DGraphicSet gset(kBlack, drawtype, drawsize);
1585  //gset.marker_style = 7;
1586  TVector3 last_point;
1587  int last_track=-1;
1588  int last_part=-1;
1589  for(unsigned int i=0; i<mctrajectorypoints.size(); i++){
1590  const DMCTrajectoryPoint *pt = mctrajectorypoints[i];
1591 
1592  switch(pt->part){
1593  case Gamma:
1594  if(!hdvmf->GetCheckButton("trajectories_photon"))continue;
1595  break;
1596  case Electron:
1597  if(!hdvmf->GetCheckButton("trajectories_electron"))continue;
1598  break;
1599  case Positron:
1600  if(!hdvmf->GetCheckButton("trajectories_positron"))continue;
1601  break;
1602  case Proton:
1603  if(!hdvmf->GetCheckButton("trajectories_proton"))continue;
1604  break;
1605  case Neutron:
1606  if(!hdvmf->GetCheckButton("trajectories_neutron"))continue;
1607  break;
1608  case PiPlus:
1609  if(!hdvmf->GetCheckButton("trajectories_piplus"))continue;
1610  break;
1611  case PiMinus:
1612  if(!hdvmf->GetCheckButton("trajectories_piminus"))continue;
1613  break;
1614  default:
1615  if(!hdvmf->GetCheckButton("trajectories_other"))continue;
1616  break;
1617  }
1618 
1619  TVector3 v(pt->x, pt->y, pt->z);
1620 
1621  if(i>0){
1622  //if((v-last_point).Mag() > 10.0){
1623  if(pt->track!=last_track || pt->part!=last_part){
1624  if(hdvmf->GetCheckButton("trajectories_colors")){
1625  switch(last_part){
1626  case Gamma:
1627  gset.color = kOrange;
1628  break;
1629  case Electron:
1630  case PiMinus:
1631  gset.color = kRed+2;
1632  break;
1633  case Positron:
1634  case Proton:
1635  case PiPlus:
1636  gset.color = kBlue+1;
1637  break;
1638  case Neutron:
1639  gset.color = kGreen+2;
1640  break;
1641  default:
1642  gset.color = kBlack;
1643  break;
1644  }
1645  }else{
1646  gset.color = kBlack;
1647  }
1648  graphics.push_back(gset);
1649  gset.points.clear();
1650  }
1651  }
1652 
1653  gset.points.push_back(v);
1654  last_point = v;
1655  last_track = pt->track;
1656  last_part = pt->part;
1657  }
1658 
1659  if(hdvmf->GetCheckButton("trajectories_colors")){
1660  switch(last_part){
1661  case Gamma:
1662  gset.color = kOrange;
1663  break;
1664  case Electron:
1665  case PiMinus:
1666  gset.color = kRed+2;
1667  break;
1668  case Positron:
1669  case Proton:
1670  case PiPlus:
1671  gset.color = kBlue+1;
1672  break;
1673  case Neutron:
1674  gset.color = kGreen+2;
1675  break;
1676  default:
1677  gset.color = kBlack;
1678  break;
1679  }
1680  }else{
1681  gset.color = kBlack;
1682  }
1683  graphics.push_back(gset);
1684  }
1685 
1686  // DTrackCandidate
1687  if(hdvmf->GetCheckButton("candidates")){
1688  vector<const DTrackCandidate*> trackcandidates;
1689  loop->Get(trackcandidates, hdvmf->GetFactoryTag("DTrackCandidate"));
1690  for(unsigned int i=0; i<trackcandidates.size(); i++){
1691  int color=i+1;
1692  double size=2.0;
1693  //if(trackcandidates[i]->charge() >0.0) color += 100; // lighter shade
1694  //if(trackcandidates[i]->charge() <0.0) color += 150; // darker shade
1695  AddKinematicDataTrack(trackcandidates[i], color, size);
1696  }
1697  }
1698 
1699  // DTrackWireBased
1700  if(hdvmf->GetCheckButton("wiretracks")){
1701  vector<const DTrackWireBased*> wiretracks;
1702  loop->Get(wiretracks, hdvmf->GetFactoryTag("DTrackWireBased"));
1703  for(unsigned int i=0; i<wiretracks.size(); i++){
1704  AddKinematicDataTrack(wiretracks[i], (wiretracks[i]->charge()>0.0 ? kBlue:kRed)+2, 1.25);
1705  }
1706  }
1707 
1708  // DTrackTimeBased
1709  if(hdvmf->GetCheckButton("timetracks")){
1710  vector<const DTrackTimeBased*> timetracks;
1711  loop->Get(timetracks, hdvmf->GetFactoryTag("DTrackTimeBased"));
1712  for(unsigned int i=0; i<timetracks.size(); i++){
1713  AddKinematicDataTrack(timetracks[i], (timetracks[i]->charge()>0.0 ? kBlue:kRed)+0, 1.00);
1714  }
1715  }
1716 
1717  // DChargedTrack
1718  if(hdvmf->GetCheckButton("chargedtracks")){
1719  vector<const DChargedTrack*> chargedtracks;
1720  loop->Get(chargedtracks, hdvmf->GetFactoryTag("DChargedTrack"));
1721  for(unsigned int i=0; i<chargedtracks.size(); i++){
1722  int color=kViolet-3;
1723  double size=1.25;
1724  if (chargedtracks[i]->Get_Charge() > 0) color=kMagenta;
1725 
1726  if (chargedtracks[i]->Get_BestFOM()->mass() > 0.9) size=2.5;
1727  AddKinematicDataTrack(chargedtracks[i]->Get_BestFOM(),color,size);
1728  }
1729  }
1730 
1731  // DNeutralParticles
1732  if(hdvmf->GetCheckButton("neutrals")){
1733  vector<const DNeutralParticle*> photons;
1734  loop->Get(photons, hdvmf->GetFactoryTag("DNeutralParticle"));
1735 
1736  for(unsigned int i=0; i<photons.size(); i++){
1737  int color = kBlack;
1738  auto locNeutralShower = photons[i]->dNeutralShower;
1739  DetectorSystem_t locDetSys = locNeutralShower->dDetectorSystem;
1740  if(locDetSys==SYS_FCAL)color = kOrange;
1741  if(locDetSys==SYS_BCAL)color = kYellow+2;
1742  //if(locDetSys==DPhoton::kCharge)color = kRed;
1743  AddKinematicDataTrack(photons[i]->Get_BestFOM(), color, 1.00);
1744  }
1745  }
1746 }
1747 
1749 {
1750  // This routine is run only once when the canvas is created.
1751  BCALHitCanvas = hdvmf->GetBcalDispFrame();
1752  BCALHitMatrixU = new TH2F("BCALHitMatrixU","BCAL Hits Upstream Energy;Sector number;Layer;Energy (MeV)", 48*4+2, -1.5, 192.5, 10, 0., 10.);
1753  BCALHitMatrixD = new TH2F("BCALHitMatrixD","BCAL Hits Downstream Energy;Sector number;Layer;Energy (MeV)",48*4+2, -1.5, 192.5, 10, 0., 10.);
1754  BCALParticles = new TH2F("BCALParticles","BCAL Particle Hits Type;Phi angle [deg];;Particle Momentum",(48*4+2)*4, -1.87, 361.87, 1, 0., 1.);
1755  BCALHitMatrixU->SetStats(0);
1756  BCALHitMatrixD->SetStats(0);
1757  BCALParticles->SetStats(0);
1758  Float_t size = 0.06;
1759  BCALHitMatrixU->GetXaxis()->SetTitleSize(size);
1760  BCALHitMatrixU->GetXaxis()->SetTitleOffset(0.8);
1761  BCALHitMatrixU->GetXaxis()->SetLabelSize(size);
1762  BCALHitMatrixU->GetYaxis()->SetTitleSize(size);
1763  BCALHitMatrixU->GetYaxis()->SetTitleOffset(0.4);
1764  BCALHitMatrixU->GetYaxis()->SetLabelSize(size);
1765  BCALHitMatrixU->GetZaxis()->SetTitleSize(size);
1766  BCALHitMatrixU->GetZaxis()->SetTitleOffset(0.4);
1767  BCALHitMatrixU->GetZaxis()->SetLabelSize(size);
1768  BCALHitMatrixD->GetXaxis()->SetTitleSize(size);
1769  BCALHitMatrixD->GetXaxis()->SetTitleOffset(0.8);
1770  BCALHitMatrixD->GetXaxis()->SetLabelSize(size);
1771  BCALHitMatrixD->GetYaxis()->SetTitleSize(size);
1772  BCALHitMatrixD->GetYaxis()->SetTitleOffset(0.4);
1773  BCALHitMatrixD->GetYaxis()->SetLabelSize(size);
1774  BCALHitMatrixD->GetZaxis()->SetTitleSize(size);
1775  BCALHitMatrixD->GetZaxis()->SetTitleOffset(0.4);
1776  BCALHitMatrixD->GetZaxis()->SetLabelSize(size);
1777  BCALParticles->GetXaxis()->SetTitleSize(size);
1778  BCALParticles->GetXaxis()->SetTitleOffset(0.8);
1779  BCALParticles->GetXaxis()->SetLabelSize(size);
1780  BCALParticles->GetYaxis()->SetTitleSize(size);
1781  BCALParticles->GetYaxis()->SetTitleOffset(0.4);
1782  BCALParticles->GetYaxis()->SetLabelSize(size);
1783  BCALParticles->GetZaxis()->SetTitleSize(size);
1784  BCALParticles->GetZaxis()->SetTitleOffset(0.4);
1785  BCALParticles->GetZaxis()->SetLabelSize(size);
1786 
1787  if (BCALHitCanvas) {
1788  vector<const DBCALHit*> locBcalHits;
1789  loop->Get(locBcalHits);
1790  BCALHitMatrixU->Reset();
1791  BCALHitMatrixD->Reset();
1792  for (unsigned int k=0;k<locBcalHits.size();k++){
1793 
1794  const DBCALHit* hit = locBcalHits[k];
1795  float idxY = (float)hit->layer-1;
1796  float idxX = (float) (hit->sector-1 + (hit->module-1)*4);
1797  if (hit->end == DBCALGeometry::kUpstream){
1798  if (hit->layer==1){
1799  BCALHitMatrixU->Fill(idxX,idxY,hit->E);
1800  } else if (hit->layer==2){
1801  BCALHitMatrixU->Fill(idxX,idxY,hit->E);
1802  BCALHitMatrixU->Fill(idxX,idxY+1.,hit->E);
1803  } else if (hit->layer==3){
1804  BCALHitMatrixU->Fill(idxX,idxY+1,hit->E);
1805  BCALHitMatrixU->Fill(idxX,idxY+2.,hit->E);
1806  BCALHitMatrixU->Fill(idxX,idxY+3.,hit->E);
1807  } else if (hit->layer==4){
1808  BCALHitMatrixU->Fill(idxX,idxY+3,hit->E);
1809  BCALHitMatrixU->Fill(idxX,idxY+4.,hit->E);
1810  BCALHitMatrixU->Fill(idxX,idxY+5.,hit->E);
1811  BCALHitMatrixU->Fill(idxX,idxY+6.,hit->E);
1812  }
1813  } else {
1814  if (hit->layer==1){
1815  BCALHitMatrixD->Fill(idxX,idxY,hit->E);
1816  } else if (hit->layer==2){
1817  BCALHitMatrixD->Fill(idxX,idxY,hit->E);
1818  BCALHitMatrixD->Fill(idxX,idxY+1.,hit->E);
1819  } else if (hit->layer==3){
1820  BCALHitMatrixD->Fill(idxX,idxY+1,hit->E);
1821  BCALHitMatrixD->Fill(idxX,idxY+2.,hit->E);
1822  BCALHitMatrixD->Fill(idxX,idxY+3.,hit->E);
1823  } else if (hit->layer==4){
1824  BCALHitMatrixD->Fill(idxX,idxY+3,hit->E);
1825  BCALHitMatrixD->Fill(idxX,idxY+4.,hit->E);
1826  BCALHitMatrixD->Fill(idxX,idxY+5.,hit->E);
1827  BCALHitMatrixD->Fill(idxX,idxY+6.,hit->E);
1828  }
1829  }
1830  }
1831 
1832  LayerLegend = new TLegend(0.91,0.7,0.99,0.99);
1833  for (int layer=0; layer<4; layer++) {
1834  char name[255];
1835  sprintf(name,"PhiZLayer%i",layer+1);
1836  BCALPointZphiLayer[layer] = new TH2F(name,"BCAL Points vs Z position;Phi angle [deg];Z position (cm);Energy",48*4,0,360,48,-80,400);
1837  FormatHistogram(BCALPointZphiLayer[layer],layer+1);
1838  BCALPointZphiLayer[layer]->SetLineWidth(2);
1839  sprintf(name,"PhiTLayer%i",layer+1);
1840  BCALPointPhiTLayer[layer] = new TH2F(name,"BCAL Points vs time;Phi angle [deg];time (ns);Energy",48*4,0,360,50,-100,100);
1841  FormatHistogram(BCALPointPhiTLayer[layer],layer+1);
1842  BCALPointPhiTLayer[layer]->SetLineWidth(2);
1843 
1844  sprintf(name,"layer %i",layer+1);
1845  LayerLegend->AddEntry(BCALPointZphiLayer[layer],name);
1846  }
1847 
1848  vector<const DBCALIncidentParticle*> locBcalParticles;
1849  loop->Get(locBcalParticles);
1850  BCALParticles->Reset();
1851  BCALPLables.clear();
1852  for (unsigned int k=0;k<locBcalParticles.size();k++){
1853 
1854  const DBCALIncidentParticle* part = locBcalParticles[k];
1855 
1856  float p = TMath::Sqrt(part->px*part->px + part->py*part->py + part->pz*part->pz);
1857 
1858  float phi=999;
1859  if (part->x!=0){
1860  phi = TMath::ATan(TMath::Abs(part->y/part->x));
1861  //cout<<phi<<" "<<part->y<<" / "<< part->x<<endl;
1862  if (part->y>0){
1863  if (part->x<0.){
1864  phi = 3.1415926 - phi;
1865  }
1866  } else {
1867  if (part->x<0){
1868  phi += 3.1415926;
1869  } else {
1870  phi = 3.1415926*2. - phi;
1871  }
1872  }
1873 
1874  phi = phi*180./3.1415926;
1875  }
1876  BCALParticles->Fill(phi,0.5,p);
1877  char l[20];
1878  sprintf(l,"%d",part->ptype);
1879  TText *t = new TText(phi,1.01,l);
1880  t->SetTextSize(0.08);
1881  t->SetTextFont(72);
1882  t->SetTextAlign(21);
1883  BCALPLables.push_back(t);
1884  }
1885 
1886 
1887  BCALHitCanvas->Clear();
1888  BCALHitCanvas->Divide(1,3,0.001,0.001);
1889  BCALHitCanvas->cd(1);
1890  gPad->SetRightMargin(0.2);
1891  BCALHitMatrixU->Draw("colz");
1892  BCALHitCanvas->cd(2);
1893  gPad->SetRightMargin(0.2);
1894  BCALHitMatrixD->Draw("colz");
1895  BCALHitCanvas->cd(3);
1896  gPad->SetRightMargin(0.2);
1897  // BCALParticles->Draw("colz");
1898  // for (unsigned int n=0;n<BCALPLables.size();n++){
1899  // BCALPLables[n]->Draw("same");
1900  // }
1901  BCALHitCanvas->Update();
1902  }
1903 
1904 }
1905 //------------------------------------------------------------------
1906 // UpdateTrackLabels
1907 //------------------------------------------------------------------
1909 {
1910  // Get the label pointers
1911  string name, tag;
1912  map<string, vector<TGLabel*> > &thrownlabs = hdvmf->GetThrownLabels();
1913  map<string, vector<TGLabel*> > &reconlabs = hdvmf->GetReconstructedLabels();
1914  hdvmf->GetReconFactory(name, tag);
1915 
1916  // Get Thrown particles
1917  vector<const DMCThrown*> throwns;
1918  if(loop)loop->Get(throwns);
1919 
1920  // Get the track info as DKinematicData objects
1921  vector<const DKinematicData*> trks;
1922  vector<const DKinematicData*> TrksCand;
1923  vector<const DTrackWireBased*> TrksWireBased;
1924  vector<const DTrackTimeBased*> TrksTimeBased;
1925  vector<const DTrackCandidate*> cand;
1926  if(loop)loop->Get(cand);
1927  for(unsigned int i=0; i<cand.size(); i++)TrksCand.push_back(cand[i]);
1928 
1929  if(loop)loop->Get(TrksWireBased);
1930  if(loop)loop->Get(TrksTimeBased);
1931 
1932  if(name=="DChargedTrack"){
1933  vector<const DChargedTrack*> chargedtracks;
1934  if(loop)loop->Get(chargedtracks, tag.c_str());
1935  for(unsigned int i=0; i<chargedtracks.size(); i++){
1936  trks.push_back(chargedtracks[i]->Get_BestFOM());
1937  }
1938  }
1939  if(name=="DTrackTimeBased"){
1940  vector<const DTrackTimeBased*> timetracks;
1941  if(loop)loop->Get(timetracks, tag.c_str());
1942  for(unsigned int i=0; i<timetracks.size(); i++)trks.push_back(timetracks[i]);
1943  }
1944  if(name=="DTrackWireBased"){
1945  vector<const DTrackWireBased*> wiretracks;
1946  if(loop)loop->Get(wiretracks, tag.c_str());
1947  for(unsigned int i=0; i<wiretracks.size(); i++)trks.push_back(wiretracks[i]);
1948  }
1949  if(name=="DTrackCandidate"){
1950  vector<const DTrackCandidate*> candidates;
1951  if(loop)loop->Get(candidates, tag.c_str());
1952  for(unsigned int i=0; i<candidates.size(); i++)trks.push_back(candidates[i]);
1953  }
1954  if(name=="DNeutralParticle"){
1955  vector<const DNeutralParticle*> photons;
1956  if(loop)loop->Get(photons, tag.c_str());
1957  for(unsigned int i=0; i<photons.size(); i++) {
1958  trks.push_back(photons[i]->Get_BestFOM());
1959  }
1960  }
1961 
1962  // Clear all labels (i.e. draw ---- in them)
1963  map<string, vector<TGLabel*> >::iterator iter;
1964  for(iter=reconlabs.begin(); iter!=reconlabs.end(); iter++){
1965  vector<TGLabel*> &labs = iter->second;
1966  for(unsigned int i=1; i<labs.size(); i++){
1967  labs[i]->SetText("--------");
1968  }
1969  }
1970  for(iter=thrownlabs.begin(); iter!=thrownlabs.end(); iter++){
1971  vector<TGLabel*> &labs = iter->second;
1972  for(unsigned int i=1; i<labs.size(); i++){
1973  labs[i]->SetText("--------");
1974  }
1975  }
1976 
1977  // Loop over thrown particles and fill in labels
1978  int ii=0;
1979  for(unsigned int i=0; i<throwns.size(); i++){
1980  const DMCThrown *trk = throwns[i];
1981  if(trk->type==0)continue;
1982  int row = thrownlabs["trk"].size()-(ii++)-1;
1983  if(row<1)break;
1984 
1985  stringstream trkno, type, p, theta, phi, z;
1986  trkno<<setprecision(4)<<i+1;
1987  thrownlabs["trk"][row]->SetText(trkno.str().c_str());
1988 
1989  thrownlabs["type"][row]->SetText(ParticleType((Particle_t)trk->type));
1990 
1991  p<<setprecision(4)<<trk->momentum().Mag();
1992  thrownlabs["p"][row]->SetText(p.str().c_str());
1993 
1994  theta<<setprecision(4)<<trk->momentum().Theta()*TMath::RadToDeg();
1995  thrownlabs["theta"][row]->SetText(theta.str().c_str());
1996 
1997  double myphi = trk->momentum().Phi();
1998  if(myphi<0.0)myphi+=2.0*M_PI;
1999  phi<<setprecision(4)<<myphi;
2000  thrownlabs["phi"][row]->SetText(phi.str().c_str());
2001 
2002  z<<setprecision(4)<<trk->position().Z();
2003  thrownlabs["z"][row]->SetText(z.str().c_str());
2004  }
2005 
2006  // Loop over tracks and fill in labels
2007  for(unsigned int i=0; i<trks.size(); i++){
2008  const DKinematicData *trk = trks[i];
2009  int row = reconlabs["trk"].size()-i-1;
2010  if(row<1)break;
2011 
2012  stringstream trkno, type, p, theta, phi, z, chisq_per_dof, Ndof,cand;
2013  stringstream fom;
2014  trkno<<setprecision(4)<<i+1;
2015  reconlabs["trk"][row]->SetText(trkno.str().c_str());
2016 
2017  double mass = trk->mass();
2018  if(fabs(mass-0.13957)<1.0E-4)type<<"pi";
2019  else if(fabs(mass-0.93827)<1.0E-4)type<<"proton";
2020  else if(fabs(mass-0.493677)<1.0E-4)type<<"K";
2021  else if(fabs(mass-0.000511)<1.0E-4)type<<"e";
2022  else if (fabs(mass)<1.e-4 && fabs(trk->charge())<1.e-4) type << "gamma";
2023  else type<<"q=";
2024  if (fabs(trk->charge())>1.e-4){
2025  type<<(trk->charge()>0 ? "+":"-");
2026  }
2027  reconlabs["type"][row]->SetText(type.str().c_str());
2028 
2029  p<<setprecision(4)<<trk->momentum().Mag();
2030  reconlabs["p"][row]->SetText(p.str().c_str());
2031 
2032  theta<<setprecision(4)<<trk->momentum().Theta()*TMath::RadToDeg();
2033  reconlabs["theta"][row]->SetText(theta.str().c_str());
2034 
2035  phi<<setprecision(4)<<trk->momentum().Phi()*TMath::RadToDeg();
2036  reconlabs["phi"][row]->SetText(phi.str().c_str());
2037 
2038  z<<setprecision(4)<<trk->position().Z();
2039  reconlabs["z"][row]->SetText(z.str().c_str());
2040 
2041  // Get chisq and Ndof for DTrackTimeBased or DTrackWireBased objects
2042  const DTrackTimeBased *timetrack=dynamic_cast<const DTrackTimeBased*>(trk);
2043  const DTrackWireBased *track=dynamic_cast<const DTrackWireBased*>(trk);
2044 
2045  const DTrackCandidate *candidate=dynamic_cast<const DTrackCandidate*>(trk);
2046  const DChargedTrackHypothesis *chargedtrack=dynamic_cast<const DChargedTrackHypothesis *>(trk);
2047 
2048  if(timetrack){
2049  chisq_per_dof<<setprecision(4)<<timetrack->chisq/timetrack->Ndof;
2050  Ndof<<timetrack->Ndof;
2051  fom << timetrack->FOM;
2052  }else if(track){
2053  chisq_per_dof<<setprecision(4)<<track->chisq/track->Ndof;
2054  Ndof<<track->Ndof;
2055  fom << "N/A";
2056  }else if (candidate){
2057  chisq_per_dof<<setprecision(4)<<candidate->chisq/candidate->Ndof;
2058  Ndof<<candidate->Ndof;
2059  fom << "N/A";
2060  }
2061  else if (chargedtrack){
2062  auto locTrackTimeBased = chargedtrack->Get_TrackTimeBased();
2063  chisq_per_dof<<setprecision(4)<<locTrackTimeBased->chisq/locTrackTimeBased->Ndof;
2064  Ndof<<locTrackTimeBased->Ndof;
2065  fom << locTrackTimeBased->FOM;
2066  }
2067  else{
2068  chisq_per_dof << "--------";
2069  Ndof << "--------";
2070  fom << "--------";
2071  }
2072 
2073  reconlabs["chisq/Ndof"][row]->SetText(chisq_per_dof.str().c_str());
2074  reconlabs["Ndof"][row]->SetText(Ndof.str().c_str());
2075  reconlabs["FOM"][row]->SetText(fom.str().c_str());
2076 
2077  if (timetrack){
2078  cand << timetrack->candidateid;
2079  }
2080  else if (track){
2081  cand << track->candidateid;
2082  }
2083  else {
2084  cand << "--------";
2085  }
2086  reconlabs["cand"][row]->SetText(cand.str().c_str());
2087  }
2088 
2089  // Have the pop-up window with the full particle list update it's labels
2090  if( fulllistmf ) fulllistmf->UpdateTrackLabels(throwns, trks);
2091  if( debugermf ) {
2092  debugermf->SetTrackCandidates(TrksCand);
2093  debugermf->SetTrackWireBased(TrksWireBased);
2094  debugermf->SetTrackTimeBased(TrksTimeBased);
2095  debugermf->UpdateTrackLabels();
2096  }
2097 }
2098 
2099 //------------------------------------------------------------------
2100 // AddKinematicDataTrack
2101 //------------------------------------------------------------------
2102 void MyProcessor::AddKinematicDataTrack(const DKinematicData* kd, int color, double size)
2103 {
2104  // Create a reference trajectory with the given kinematic data and swim
2105  // it through the detector.
2106  DReferenceTrajectory rt(Bfield);
2107  rt.Rsqmax_interior = RMAX_INTERIOR*RMAX_INTERIOR;
2108  rt.Rsqmax_exterior = RMAX_EXTERIOR*RMAX_EXTERIOR;
2109 
2110  if(MATERIAL_MAP_MODEL=="DRootGeom"){
2111  rt.SetDRootGeom(RootGeom);
2112  rt.SetDGeometry(NULL);
2113  }else if(MATERIAL_MAP_MODEL=="DGeometry"){
2114  rt.SetDRootGeom(NULL);
2115  rt.SetDGeometry(geom);
2116  }else if(MATERIAL_MAP_MODEL!="NONE"){
2117  _DBG_<<"WARNING: Invalid value for TRKFIT:MATERIAL_MAP_MODEL (=\""<<MATERIAL_MAP_MODEL<<"\")"<<endl;
2118  }
2119 
2120  rt.SetMass(kd->mass());
2121  rt.Swim(kd->position(), kd->momentum(), kd->charge());
2122 
2123  // Create a new graphics set and fill it with all of the trajectory points
2124  DGraphicSet gset(color, kLine, size);
2126  for(int i=0; i<rt.Nswim_steps; i++, step++){
2127  TVector3 tpoint(step->origin.X(),step->origin.Y(),step->origin.Z());
2128  gset.points.push_back(tpoint);
2129  }
2130 
2131  // Push the graphics set onto the stack
2132  graphics.push_back(gset);
2133 }
2134 
2135 //------------------------------------------------------------------
2136 // GetIntersectionWithCalorimeter
2137 //------------------------------------------------------------------
2139 {
2140  // Create a reference trajectory with the given kinematic data and swim
2141  // it through the detector.
2142  DReferenceTrajectory rt(Bfield);
2143  rt.Rsqmax_interior = RMAX_INTERIOR*RMAX_INTERIOR;
2144  rt.Rsqmax_exterior = RMAX_EXTERIOR*RMAX_EXTERIOR;
2145 
2146  if(MATERIAL_MAP_MODEL=="DRootGeom"){
2147  rt.SetDRootGeom(RootGeom);
2148  rt.SetDGeometry(NULL);
2149  }else if(MATERIAL_MAP_MODEL=="DGeometry"){
2150  rt.SetDRootGeom(NULL);
2151  rt.SetDGeometry(geom);
2152  }else if(MATERIAL_MAP_MODEL!="NONE"){
2153  _DBG_<<"WARNING: Invalid value for TRKFIT:MATERIAL_MAP_MODEL (=\""<<MATERIAL_MAP_MODEL<<"\")"<<endl;
2154  }
2155 
2156  rt.SetMass(kd->mass());
2157  rt.Swim(kd->position(), kd->momentum(), kd->charge());
2158 
2159  // Find intersection with FCAL
2160  DVector3 pos_fcal;
2161  double s_fcal = 1.0E6;
2162  DVector3 origin(0.0, 0.0, FCAL_Zmin);
2163  DVector3 norm(0.0, 0.0, -1.0);
2164  //rt.GetIntersectionWithPlane(origin, norm, pos_fcal, &s_fcal); // This gives different answer than below!
2165  DVector3 p_at_intersection;
2166  rt.GetIntersectionWithPlane(origin, norm, pos_fcal, p_at_intersection, &s_fcal, NULL, NULL, SYS_FCAL);
2167  if(pos_fcal.Perp()<FCAL_Rmin || pos_fcal.Perp()>FCAL_Rmax || !isfinite(pos_fcal.Z()))s_fcal = 1.0E6;
2168 
2169  // Find intersection with BCAL
2170  DVector3 pos_bcal;
2171  double s_bcal = 1.0E6;
2172  rt.GetIntersectionWithRadius(BCAL_Rmin, pos_bcal, &s_bcal);
2173  if(pos_bcal.Z()<BCAL_Zmin || pos_bcal.Z()>(BCAL_Zmin+BCAL_Zlen) || !isfinite(pos_bcal.Z()))s_bcal = 1.0E6;
2174 
2175  if(s_fcal>10000.0 && s_bcal>10000.0){
2176  // neither calorimeter hit
2177  who = SYS_NULL;
2178  pos.SetXYZ(0.0,0.0,0.0);
2179  }else if(s_fcal<s_bcal){
2180  // FCAL hit
2181  who = SYS_FCAL;
2182  pos = pos_fcal;
2183  }else{
2184  // BCAL hit
2185  who = SYS_BCAL;
2186  pos = pos_bcal;
2187  }
2188 }
2189 
2190 //------------------------------------------------------------------
2191 // GetFactoryNames
2192 //------------------------------------------------------------------
2193 void MyProcessor::GetFactoryNames(vector<string> &facnames)
2194 {
2195  vector<JEventLoop*> loops = app->GetJEventLoops();
2196  if(loops.size()>0){
2197  vector<string> facnames;
2198  loops[0]->GetFactoryNames(facnames);
2199  }
2200 }
2201 
2202 //------------------------------------------------------------------
2203 // GetFactories
2204 //------------------------------------------------------------------
2205 void MyProcessor::GetFactories(vector<JFactory_base*> &factories)
2206 {
2207  vector<JEventLoop*> loops = app->GetJEventLoops();
2208  if(loops.size()>0){
2209  factories = loops[0]->GetFactories();
2210  }
2211 }
2212 
2213 //------------------------------------------------------------------
2214 // GetNrows
2215 //------------------------------------------------------------------
2216 unsigned int MyProcessor::GetNrows(const string &factory, string tag)
2217 {
2218  vector<JEventLoop*> loops = app->GetJEventLoops();
2219  if(loops.size()>0){
2220  // Here is something a little tricky. The GetFactory() method of JEventLoop
2221  // gets the factory of the specified data name and tag, but without trying
2222  // to substitute a user-specified tag (a'la -PDEFTAG:XXX=YYY) as is done
2223  // on normal Get() method calls. Therefore, we have to check for the default
2224  // tags ourselves and substitute it "by hand".
2225  if(tag==""){
2226  map<string,string> default_tags = loops[0]->GetDefaultTags();
2227  map<string, string>::const_iterator iter = default_tags.find(factory);
2228  if(iter!=default_tags.end())tag = iter->second.c_str();
2229  }
2230  JFactory_base *fac = loops[0]->GetFactory(factory, tag.c_str());
2231 
2232  // Since calling GetNrows will cause the program to quit if there is
2233  // not a valid event, then first check that there is one before calling it
2234  if(loops[0]->GetJEvent().GetJEventSource() == NULL)return 0;
2235 
2236  return fac==NULL ? 0:(unsigned int)fac->GetNrows();
2237  }
2238 
2239  return 0;
2240 }
2241 
2242 //------------------------------------------------------------------
2243 // GetDReferenceTrajectory
2244 //------------------------------------------------------------------
2245 void MyProcessor::GetDReferenceTrajectory(string dataname, string tag, unsigned int index, DReferenceTrajectory* &rt, vector<const DCDCTrackHit*> &cdchits)
2246 {
2247 _DBG__;
2248  // initialize rt to NULL in case we don't find the one requested
2249  rt = NULL;
2250  cdchits.clear();
2251 
2252  // Get pointer to the JEventLoop so we can get at the data
2253  vector<JEventLoop*> loops = app->GetJEventLoops();
2254  if(loops.size()==0)return;
2255  JEventLoop* &loop = loops[0];
2256 
2257  // Variables to hold track parameters
2258  DVector3 pos, mom(0,0,0);
2259  double q=0.0;
2260  double mass = 0.13957018;
2261 
2262  // Find the specified track
2263  if(dataname=="DChargedTrack"){
2264  vector<const DChargedTrack*> chargedtracks;
2265  vector<const DTrackTimeBased*> timebasedtracks;
2266  loop->Get(chargedtracks, tag.c_str());
2267  if(index>=chargedtracks.size())return;
2268  q = chargedtracks[index]->Get_Charge();
2269  pos = chargedtracks[index]->Get_BestFOM()->position();
2270  mom = chargedtracks[index]->Get_BestFOM()->momentum();
2271  chargedtracks[index]->Get_BestFOM()->GetT(timebasedtracks);
2272  if(!timebasedtracks.empty()){
2273  timebasedtracks[0]->Get(cdchits);
2274  mass = chargedtracks[index]->Get_BestFOM()->mass();
2275  }
2276  }
2277 
2278  if(dataname=="DTrackTimeBased"){
2279  vector<const DTrackTimeBased*> timetracks;
2280  loop->Get(timetracks, tag.c_str());
2281  if(index>=timetracks.size())return;
2282  q = timetracks[index]->charge();
2283  pos = timetracks[index]->position();
2284  mom = timetracks[index]->momentum();
2285  timetracks[index]->Get(cdchits);
2286  mass = timetracks[index]->mass();
2287  }
2288 
2289  if(dataname=="DTrackWireBased"){
2290  vector<const DTrackWireBased*> wiretracks;
2291  loop->Get(wiretracks, tag.c_str());
2292  if(index>=wiretracks.size())return;
2293  q = wiretracks[index]->charge();
2294  pos = wiretracks[index]->position();
2295  mom = wiretracks[index]->momentum();
2296  wiretracks[index]->Get(cdchits);
2297  mass = wiretracks[index]->mass();
2298  }
2299 
2300  if(dataname=="DTrackCandidate"){
2301  vector<const DTrackCandidate*> tracks;
2302  loop->Get(tracks, tag.c_str());
2303  if(index>=tracks.size())return;
2304  q = tracks[index]->charge();
2305  pos = tracks[index]->position();
2306  mom = tracks[index]->momentum();
2307  tracks[index]->Get(cdchits);
2308  mass = tracks[index]->mass();
2309  }
2310 
2311  if(dataname=="DMCThrown"){
2312  vector<const DMCThrown*> tracks;
2313  loop->Get(tracks, tag.c_str());
2314  if(index>=tracks.size())return;
2315  const DMCThrown *t = tracks[index];
2316  q = t->charge();
2317  pos = t->position();
2318  mom = t->momentum();
2319  tracks[index]->Get(cdchits);
2320  mass = tracks[index]->mass();
2321 _DBG_<<"mass="<<mass<<endl;
2322  }
2323 
2324  // Make sure we found a charged particle we can track
2325  if(q==0.0 || mom.Mag()<0.01)return;
2326 
2327  // Create a new DReference trajectory object. The caller takes
2328  // ownership of this and so they are responsible for deleting it.
2329  rt = new DReferenceTrajectory(Bfield);
2330  rt->Rsqmax_interior = RMAX_INTERIOR*RMAX_INTERIOR;
2331  rt->Rsqmax_exterior = RMAX_EXTERIOR*RMAX_EXTERIOR;
2332  rt->SetMass(mass);
2333  if(MATERIAL_MAP_MODEL=="DRootGeom"){
2334  rt->SetDRootGeom(RootGeom);
2335  rt->SetDGeometry(NULL);
2336  }else if(MATERIAL_MAP_MODEL=="DGeometry"){
2337  rt->SetDRootGeom(NULL);
2338  rt->SetDGeometry(geom);
2339  }else if(MATERIAL_MAP_MODEL!="NONE"){
2340  _DBG_<<"WARNING: Invalid value for TRKFIT:MATERIAL_MAP_MODEL (=\""<<MATERIAL_MAP_MODEL<<"\")"<<endl;
2341  }
2342  rt->Swim(pos, mom, q);
2343 }
2344 
2345 //------------------------------------------------------------------
2346 // GetAllWireHits
2347 //------------------------------------------------------------------
2348 void MyProcessor::GetAllWireHits(vector<pair<const DCoordinateSystem*,double> > &allhits)
2349 {
2350  /// Argument is vector of pairs that contain a pointer to the
2351  /// DCoordinateSystem representing a wire and a double that
2352  /// represents the drift distance. To get info on the specific
2353  /// wire, one needs to attempt a dynamic_cast to both a DCDCWire
2354  /// and a DFDCWire and access the parameters of whichever one succeeds.
2355 
2356  // Get pointer to the JEventLoop so we can get at the data
2357  vector<JEventLoop*> loops = app->GetJEventLoops();
2358  if(loops.size()==0)return;
2359  JEventLoop* &loop = loops[0];
2360 
2361  // Get CDC wire hits
2362  vector<const DCDCTrackHit*> cdchits;
2363  loop->Get(cdchits);
2364  for(unsigned int i=0; i<cdchits.size(); i++){
2365  pair<const DCoordinateSystem*,double> hit;
2366  hit.first = cdchits[i]->wire;
2367  hit.second = cdchits[i]->dist;
2368  allhits.push_back(hit);
2369  }
2370 
2371  // Get FDC wire hits
2372  vector<const DFDCPseudo*> fdchits;
2373  loop->Get(fdchits);
2374  for(unsigned int i=0; i<fdchits.size(); i++){
2375  pair<const DCoordinateSystem*,double> hit;
2376  hit.first = fdchits[i]->wire;
2377  hit.second = 0.0055*fdchits[i]->time;
2378  allhits.push_back(hit);
2379  }
2380 }
2381 
2382 
2383 
2384 //------------------------------------------------------------------
2385 // FormatHistogram
2386 //------------------------------------------------------------------
2387 void MyProcessor::FormatHistogram(TH2* histo, int color=1)
2388 {
2389  histo->SetStats(0);
2390  histo->SetLineColor(color);
2391  histo->SetMarkerColor(color);
2392  Float_t size = 0.06;
2393  histo->SetTitleSize(size);
2394  histo->GetXaxis()->SetTitleSize(size);
2395  histo->GetXaxis()->SetTitleOffset(0.8);
2396  histo->GetXaxis()->SetLabelSize(size);
2397  histo->GetYaxis()->SetTitleSize(size);
2398  histo->GetYaxis()->SetTitleOffset(0.5);
2399  histo->GetYaxis()->SetLabelSize(size);
2400  histo->GetZaxis()->SetTitleSize(size);
2401  histo->GetZaxis()->SetTitleOffset(0.4);
2402  histo->GetZaxis()->SetLabelSize(size);
2403 }
void SetNTrWireBased(Int_t d)
hdv_fulllistframe * GetFullListFrame(void)
Definition: track.h:16
void SetNTrCand(Int_t d)
void SetReconstructedFactories(vector< string > &facnames)
float E
Definition: DBCALHit.h:30
DApplication * dapp
void UpdateBcalDisp(void)
int row
Definition: DCCALHit.h:22
float chisq
Chi-squared for the track (not chisq/dof!)
jerror_t init(void)
int end
Definition: DTOFHit.h:23
int module() const
Definition: DBCALPoint.h:64
TCanvas * GetBcalDispFrame(void)
float chisq
Chi-squared for the track (not chisq/dof!)
void SetTrig(char *trigstring)
int type
Definition: DFDCHit.h:44
Int_t layer
TVector2 DVector2
Definition: DVector2.h:9
int column
Definition: DCCALHit.h:23
bool GetCheckButton(string who)
TVector3 DVector3
Definition: DVector3.h:14
Definition: GlueX.h:16
int plane
Definition: DTOFHit.h:21
bool has_fADC
Definition: DTOFHit.h:29
uint32_t trig_mask
Definition: DL1Trigger.h:18
double energy(void) const
float rho() const
Definition: DBCALCluster.h:59
int module
Definition: DBCALHit.h:25
Double_t x[NCHANNELS]
Definition: st_tw_resols.C:39
sprintf(text,"Post KinFit Cut")
TPolyLine * GetCCALPolyLine(int row, int col)
#define y
DMagneticFieldMap * GetBfield(unsigned int run_number=1)
void SetEvent(ULong64_t id)
float phi() const
Definition: DBCALPoint.h:56
int sector() const
Definition: DBCALPoint.h:66
void SetMass(double mass)
JEventLoop * eventloop
Definition: hdview2.cc:16
const DVector3 & position(void) const
int gLayer
Definition: DFDCHit.h:28
void Swim(const DVector3 &pos, const DVector3 &mom, double q=-1000.0, const TMatrixFSym *cov=NULL, double smax=2000.0, const DCoordinateSystem *wire=NULL)
const DTrackTimeBased * Get_TrackTimeBased(void) const
void SetTimeBasedTrackFactories(vector< string > &facnames)
void UpdateTrackLabels(void)
static char * ParticleType(Particle_t p)
Definition: particleType.h:142
float t
Definition: DTOFHit.h:28
void FillGraphics(void)
void SetNTrTimeBased(Int_t d)
int layer
Definition: DBCALHit.h:26
DetectorSystem_t
Definition: GlueX.h:15
float t() const
Definition: DBCALCluster.h:49
void AddKinematicDataTrack(const DKinematicData *kd, int color, double size)
static float FCAL_Rmin
trig[33-1]
static char index(char c)
Definition: base64.cpp:115
jerror_t brun(JEventLoop *eventLoop, int32_t runnumber)
Called once at program start.
DVector2 positionOnFace(int row, int column) const
oid_t candidateid
id of DTrackCandidate corresponding to this track
static float BCAL_Rmin
int layer() const
Definition: DBCALPoint.h:65
void SetDGeometry(const DGeometry *geom)
Definition: GlueX.h:17
Definition: GlueX.h:19
JApplication * japp
TF1 * f
Definition: FitGains.C:21
DetectorSystem_t system
particle type
Definition: DMCTrackHit.h:25
DetectorSystem_t dDetectorSystem
DLorentzVector dSpacetimeVertex
static float BCAL_Zmin
DBCALGeometry::End end
Definition: DBCALHit.h:28
MyProcessor * gMYPROC
DRootGeom * GetRootGeom(unsigned int run_number)
TEllipse * e
map< string, vector< TGLabel * > > & GetThrownLabels(void)
float E() const
Definition: DBCALPoint.h:38
DVector3 pos
Definition: DTOFPoint.h:33
float z
coordinates of hit in cm and rad
Definition: DMCTrackHit.h:20
int GetColor(int number)
Definition: GlueX.h:20
Definition: GlueX.h:18
TPolyLine * GetBCALPolyLine(int mod, int layer, int sector)
int row
Definition: DFCALHit.h:23
bool DMCTrajectoryPoint_track_cmp(const DMCTrajectoryPoint *a, const DMCTrajectoryPoint *b)
TPolyLine * GetFCALPolyLine(int channel)
DGeometry * GetDGeometry(unsigned int run_number)
Definition: GlueX.h:22
#define _DBG_
Definition: HDEVIO.h:12
void FormatHistogram(TH2 *, int)
void SetChargedTrackFactories(vector< string > &facnames)
void GetDReferenceTrajectory(string dataname, string tag, unsigned int index, DReferenceTrajectory *&rt, vector< const DCDCTrackHit * > &cdchits)
static vector< vector< DFDCWire * > > fdcwires
float z() const
Definition: DBCALPoint.h:60
void SetDRootGeom(const DRootGeom *RootGeom)
#define MaxWireTracks
int Ndof
Number of degrees of freedom in the fit.
double charge(void) const
void GetReconFactory(string &name, string &tag)
bool has_TDC
Definition: DTOFHit.h:30
int element
Definition: DFDCHit.h:25
void GetFactories(vector< JFactory_base * > &factories)
void GetAllWireHits(vector< pair< const DCoordinateSystem *, double > > &allhits)
#define _DBG__
Definition: HDEVIO.h:13
void SetRun(Int_t id)
int sector
Definition: DBCALHit.h:27
float y
Definition: DFCALHit.h:26
void GetIntersectionWithCalorimeter(const DKinematicData *kd, DVector3 &pos, DetectorSystem_t &who)
&lt;A href=&quot;index.html#legend&quot;&gt; &lt;IMG src=&quot;CORE.png&quot; width=&quot;100&quot;&gt; &lt;/A&gt;
static float BCAL_Zlen
int Ndof
Number of degrees of freedom in the fit.
jerror_t GetIntersectionWithRadius(double R, DVector3 &mypos, double *s=NULL, double *t=NULL, DVector3 *dir=NULL) const
map< string, vector< TGLabel * > > & GetReconstructedLabels(void)
double sqrt(double)
float E
Definition: DCCALHit.h:26
void DoMyRedraw(void)
hdv_debugerframe * GetDebugerFrame(void)
double sin(double)
float chisq
Chi-squared for the track (not chisq/dof!)
float phi() const
Definition: DBCALCluster.h:65
float theta() const
Definition: DBCALCluster.h:62
const DVector3 & momentum(void) const
static float FCAL_Zmin
float t
Definition: DFDCHit.h:32
int bar
Definition: DTOFHit.h:22
bool GetFDCWires(vector< vector< DFDCWire * > > &fdcwires) const
Definition: DGeometry.cc:1078
oid_t candidateid
which DTrackCandidate this came from
float t() const
Definition: DBCALPoint.h:41
File: DTOFHit.h Created: Tue Jan 18 16:15:26 EST 2011 Creator: B. Zihlmann Purpose: Container class t...
Definition: DTOFHit.h:16
float x
Definition: DFCALHit.h:25
float E
Definition: DFCALHit.h:27
void GetFactoryNames(vector< string > &facnames)
int Ndof
Number of degrees of freedom in the fit.
hdv_mainframe * hdvmf
#define MaxTimeTracks
int type
GEANT particle ID.
Definition: DMCThrown.h:20
printf("string=%s", string)
jerror_t evnt(JEventLoop *eventLoop, uint64_t eventnumber)
Called every event.
int column
Definition: DFCALHit.h:24
static float FCAL_Rmax
void SetSource(string source)
void SetWireBasedTrackFactories(vector< string > &facnames)
TPolyLine * GetTOFPolyLine(int translate_side, int tof_ch)
void SetCandidateFactories(vector< string > &facnames)
const char * GetFactoryTag(string who)
jerror_t GetIntersectionWithPlane(const DVector3 &origin, const DVector3 &norm, DVector3 &pos, double *s=NULL, double *t=NULL, double *var_t=NULL, DetectorSystem_t detector=SYS_NULL) const
float stereo
Definition: DCDCWire.h:82
double mass(void) const
class DFDCHit: definition for a basic FDC hit data type.
Definition: DFDCHit.h:20
Particle_t
Definition: particleType.h:12
float E() const
Definition: DBCALCluster.h:41
unsigned int GetNrows(const string &factory, string tag)