12 #include <TLorentzVector.h>
13 #include <TLorentzRotation.h>
38 #include "JANA/JGeometry.h"
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");
107 gPARMS->SetDefaultParameter(
"BCALVERBOSE", BCALVERBOSE,
"Verbosity level for BCAL objects and display");
130 vector<JEventLoop*> loops = app->GetJEventLoops();
133 vector<string> facnames;
134 loops[0]->GetFactoryNames(facnames);
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]");
177 MATERIAL_MAP_MODEL=
"DGeometry";
178 gPARMS->SetDefaultParameter(
"TRKFIT:MATERIAL_MAP_MODEL", MATERIAL_MAP_MODEL);
180 eventloop->GetCalib(
"PID/photon_track_matching", photon_track_matching);
181 DELTA_R_FCAL = photon_track_matching[
"DELTA_R_FCAL"];
191 if(!eventLoop)
return NOERROR;
193 last_jevent.FreeEvent();
194 last_jevent = loop->GetJEvent();
200 loop->GetSingle(trig);
209 string source =
"<no source>";
210 if(last_jevent.GetJEventSource())source = last_jevent.GetJEventSource()->GetSourceName();
212 cout<<
"----------- New Event "<<eventnumber<<
" (run " << last_jevent.GetRunNumber()<<
") -------------"<<endl;
218 japp->SetSequentialEventComplete();
239 graphics_xyA.clear();
240 graphics_xyB.clear();
243 graphics_tof_hits.clear();
247 vector<const DSCHit *>schits;
249 vector<const DTrackCandidate*> trCand;
251 vector<const DTrackTimeBased*> trTB;
253 vector<const DTrackWireBased*> trWB;
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++){
267 const DBCALHit* hit = locBcalHits[k];
268 float idxY = (float)hit->
layer-1;
269 float idxX = (
float) (hit->
sector-1 + (hit->
module-1)*4);
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);
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);
306 vector<const DBCALPoint*> locBcalPoint;
307 loop->Get(locBcalPoint);
309 BCALPointZphiLayer[
layer]->Reset();
310 BCALPointPhiTLayer[
layer]->Reset();
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;
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());
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());
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());
335 BCALpoints_z.clear();
336 BCALpoints_phi.clear();
338 for (
unsigned int k=0;k<locBcalPoint.size();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",
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;
352 if (point->
t()<-100) time=-99;
353 BCALPointPhiTLayer[point->
layer()-1]->Fill(pointphi,time,weight);
357 vector<const DBCALCluster*> locBcalCluster;
358 loop->Get(locBcalCluster);
360 unsigned int oldsize = BCALClusterZphiHistos.size();
361 for (
unsigned int k=0;k<oldsize;k++){
362 delete BCALClusterZphiHistos[k];
364 BCALClusterZphiHistos.clear();
366 delete ClusterLegend;
367 ClusterLegend =
new TLegend(0.91,0.01,0.99,0.99);
369 for(
unsigned int k = 0 ; k < locBcalCluster.size() ; k++){
372 vector<const DBCALPoint*> clusterpoints_range;
373 cluster_range->Get(clusterpoints_range);
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());
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());
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());
392 BCALpoints_z.clear();
393 BCALpoints_phi.clear();
395 for (
unsigned int k=0;k<locBcalCluster.size();k++){
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));
402 FormatHistogram(BCALClusterZphiHistos[k],color);
403 BCALClusterZphiHistos[k]->SetLineWidth(2);
404 sprintf(name,
"Cluster %i",k+1);
405 ClusterLegend->AddEntry(BCALClusterZphiHistos[k], name);
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);
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));
427 vector<const DBCALIncidentParticle*> locBcalParticles;
428 loop->Get(locBcalParticles);
429 BCALParticles->Reset();
431 for (
unsigned int k=0;k<locBcalParticles.size();k++){
434 float p = TMath::Sqrt(part->
px*part->
px + part->
py*part->
py + part->
pz*part->
pz);
437 phi = TMath::ATan(TMath::Abs(part->
y/part->
x));
441 phi = 3.1415926 - phi;
447 phi = 3.1415926*2. - phi;
451 phi = phi*180./3.1415926;
454 BCALParticles->Fill(phi,0.5,p);
457 TText *t =
new TText(phi,1.01,l);
458 t->SetTextSize(0.08);
461 BCALPLables.push_back(t);
466 if (BCALClusterZphiHistos.size()>0) {
467 float max = BCALClusterZphiHistos[0]->GetMaximum();
468 if (max>clustermax) clustermax=
max;
471 BCALHitCanvas->Clear();
472 BCALHitCanvas->Divide(1,3,0.001,0.001);
473 float leftmargin = 0.07;
474 BCALHitCanvas->cd(1);
477 gPad->SetLeftMargin(leftmargin);
489 if (BCALPointPhiTLayer[
layer]->GetEntries()>0) {
490 BCALPointPhiTLayer[
layer]->SetMaximum(clustermax);
492 BCALPointPhiTLayer[
layer]->Draw(
"box");
495 BCALPointPhiTLayer[
layer]->Draw(
"box,same");
501 BCALHitCanvas->cd(2);
504 gPad->SetLeftMargin(leftmargin);
513 if (BCALPointZphiLayer[
layer]->GetEntries()>0) {
514 BCALPointZphiLayer[
layer]->SetMaximum(clustermax);
516 BCALPointZphiLayer[
layer]->Draw(
"box");
519 BCALPointZphiLayer[
layer]->Draw(
"box,same");
525 BCALHitCanvas->cd(3);
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");
537 ClusterLegend->Draw();
542 BCALHitCanvas->Update();
548 vector<const DBCALHit*> bcalhits;
551 for(
unsigned int i=0; i<bcalhits.size(); i++){
558 double E = 1000.0*hit->
E;
559 double logE = log10(E);
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);
605 poly->SetLineWidth(1);
606 poly->SetFillStyle(1001);
612 vector<const DFCALHit*> fcalhits;
615 for(
unsigned int i=0; i<fcalhits.size(); i++){
621 double a = hit->
E/0.005;
622 double f =
sqrt(a>1.0 ? 1.0:a<0.0 ? 0.0:a);
628 float b = f*(1.0-grey) + grey;
632 double E = 1000*hit->
E;
634 double logE = log10(E);
676 vector<const DCCALHit*> ccalhits;
679 for(
unsigned int i=0; i<ccalhits.size(); i++){
685 double E = 1000*hit->
E;
688 double logE = log10(E);
731 double hit_north[45];
732 double hit_south[45];
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));
741 vector<const DTOFHit*> tofhits;
744 for(
unsigned int i=0; i<tofhits.size(); i++){
745 const DTOFHit *tof_hit = tofhits[i];
747 int plane = tof_hit->
plane;
748 int bar = tof_hit->
bar;
749 int end = tof_hit->
end;
750 float t = tof_hit->
t;
772 pmtPline->SetFillColor(2);
778 pmtPline->SetFillColor(2);
781 cerr <<
"Out of range TOF end" << endl;
789 pmtPline->SetFillColor(2);
795 pmtPline->SetFillColor(2);
798 cerr <<
"Out of range TOF end" << endl;
802 cerr <<
"Out of range TOF plane" << endl;
816 if (hit_down[bar]<=0 || (t < hit_down[bar]) ){
821 if (hit_up[bar]<=0 || (t < hit_up[bar]) ){
826 cerr <<
"Out of range TOF end" << endl;
831 if (hit_north[bar]<=0 || (t < hit_north[bar]) ){
836 if (hit_south[bar]<=0 || (t < hit_south[bar]) ){
841 cerr <<
"Out of range TOF end" << endl;
845 cerr <<
"Out of range TOF plane" << endl;
856 Float_t distY_Horz = -126;
857 Float_t distX_Vert = -126;
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;
866 distY_Horz = distY_Horz + 3.0;
867 distX_Vert = distX_Vert + 3.0;
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);
874 graphics_tof_hits.push_back(tdc_cir);
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);
882 graphics_tof_hits.push_back(tdc_cir);
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;
890 distY_Horz = distY_Horz + 3.0;
891 distX_Vert = distX_Vert + 3.0;
898 for (
unsigned int i=0;i<schits.size();i++){
900 double r_start=7.7493;
901 double phi0=0.2094395*(schits[i]->sector-1);
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);
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);
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);
941 graphics.push_back(gset);
946 vector<const DCDCTrackHit*> cdctrackhits;
947 loop->Get(cdctrackhits);
949 for(
unsigned int i=0; i<cdctrackhits.size(); i++){
950 const DCDCWire *wire = cdctrackhits[i]->wire;
952 int color = (cdctrackhits[i]->tdrift>-20 && cdctrackhits[i]->tdrift<800) ? kCyan:kYellow;
955 TVector3 tpoint(dpoint.X(),dpoint.Y(),dpoint.Z());
956 gset.
points.push_back(tpoint);
958 tpoint.SetXYZ(dpoint.X(),dpoint.Y(),dpoint.Z());
959 gset.
points.push_back(tpoint);
960 graphics.push_back(gset);
967 double dist1 = cdctrackhits[i]->dist;
968 TEllipse *
e =
new TEllipse(x, y, dist1, dist1);
971 graphics_xyA.push_back(e);
973 double dist2 = dist1 - 4.0*55.0E-4;
974 e =
new TEllipse(x, y, dist2, dist2);
977 graphics_xyA.push_back(e);
984 vector<const DFDCHit*> fdchits;
987 for(
unsigned int i=0; i<fdchits.size(); i++){
988 const DFDCHit *fdchit = fdchits[i];
989 if(fdchit->
type!=0)
continue;
992 _DBG_<<
"Couldn't find wire for gLayer="<<fdchit->
gLayer<<
" and element="<<fdchit->
element<<endl;
997 int color = (fdchit->
t>-50 && fdchit->
t<2000) ? kCyan:kYellow;
1000 TVector3 tpoint(dpoint.X(),dpoint.Y(),dpoint.Z());
1001 gset.
points.push_back(tpoint);
1003 tpoint.SetXYZ(dpoint.X(),dpoint.Y(),dpoint.Z());
1004 gset.
points.push_back(tpoint);
1005 graphics.push_back(gset);
1011 vector<const DFDCIntersection*> fdcints;
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);
1020 graphics.push_back(gsetp);
1025 vector<const DFDCPseudo*> fdcpseudos;
1026 loop->Get(fdcpseudos);
1029 for(
unsigned int i=0; i<fdcpseudos.size(); i++){
1030 const DFDCWire *wire = fdcpseudos[i]->wire;
1033 TVector3 pos(fdcpseudos[i]->xy.X(),
1034 fdcpseudos[i]->xy.Y(), wire->
origin.Z());
1035 gsetp.
points.push_back(pos);
1037 graphics.push_back(gsetp);
1042 vector<const DMCThrown*> mcthrown;
1043 loop->Get(mcthrown);
1044 for(
unsigned int i=0; i<mcthrown.size(); i++){
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){
1067 AddKinematicDataTrack(mcthrown[i], color, size);
1073 vector<const DMCTrackHit*> mctrackhits;
1074 loop->Get(mctrackhits);
1076 for(
unsigned int i=0; i<mctrackhits.size(); i++){
1080 TVector3 pos(hit->
r*cos(hit->
phi), hit->
r*
sin(hit->
phi), hit->
z);
1081 gset.
points.push_back(pos);
1083 graphics.push_back(gset);
1088 vector<const DMCTrackHit*> mctrackhits;
1089 loop->Get(mctrackhits);
1091 for(
unsigned int i=0; i<mctrackhits.size(); i++){
1095 TVector3 pos(hit->
r*cos(hit->
phi), hit->
r*
sin(hit->
phi), hit->
z);
1096 gset.
points.push_back(pos);
1098 graphics.push_back(gset);
1102 for(
unsigned int n=0; n<trCand.size(); n++){
1106 sprintf(str1,
"Candidate%d",n+1);
1116 AddKinematicDataTrack(trCand[n], color, 1.5);
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;
1124 TVector3 tpoint(dpoint.X(),dpoint.Y(),dpoint.Z());
1125 gset.
points.push_back(tpoint);
1127 tpoint.SetXYZ(dpoint.X(),dpoint.Y(),dpoint.Z());
1128 gset.
points.push_back(tpoint);
1129 graphics.push_back(gset);
1132 vector<const DFDCPseudo*> fdcpseudos;
1133 trCand[n]->Get(fdcpseudos,
"",1);
1136 for(
unsigned int i=0; i<fdcpseudos.size(); i++){
1137 const DFDCWire *wire = fdcpseudos[i]->wire;
1140 TVector3 pos(fdcpseudos[i]->xy.X(),
1141 fdcpseudos[i]->xy.Y(), wire->
origin.Z());
1142 gsetp.
points.push_back(pos);
1144 graphics.push_back(gsetp);
1150 for(
unsigned int n=0; n<trWB.size(); n++){
1154 sprintf(str1,
"WireBased%d",n+1);
1158 int color = trWB[n]->candidateid;
1164 AddKinematicDataTrack(trWB[n], color, 1.5);
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;
1172 TVector3 tpoint(dpoint.X(),dpoint.Y(),dpoint.Z());
1173 gset.
points.push_back(tpoint);
1175 tpoint.SetXYZ(dpoint.X(),dpoint.Y(),dpoint.Z());
1176 gset.
points.push_back(tpoint);
1177 graphics.push_back(gset);
1180 vector<const DFDCPseudo*> fdcpseudos;
1181 trWB[n]->Get(fdcpseudos,
"",1);
1184 for(
unsigned int i=0; i<fdcpseudos.size(); i++){
1185 const DFDCWire *wire = fdcpseudos[i]->wire;
1188 TVector3 pos(fdcpseudos[i]->xy.X(),
1189 fdcpseudos[i]->xy.Y(), wire->
origin.Z());
1190 gsetp.
points.push_back(pos);
1192 graphics.push_back(gsetp);
1198 for(
unsigned int n=0; n<trTB.size(); n++){
1202 sprintf(str1,
"TimeBased%d",n+1);
1206 int color = trTB[n]->candidateid;
1212 AddKinematicDataTrack(trTB[n], color, 1.5);
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;
1220 TVector3 tpoint(dpoint.X(),dpoint.Y(),dpoint.Z());
1221 gset.
points.push_back(tpoint);
1223 tpoint.SetXYZ(dpoint.X(),dpoint.Y(),dpoint.Z());
1224 gset.
points.push_back(tpoint);
1225 graphics.push_back(gset);
1228 vector<const DFDCPseudo*> fdcpseudos;
1229 trTB[n]->Get(fdcpseudos,
"",1);
1232 for(
unsigned int i=0; i<fdcpseudos.size(); i++){
1233 const DFDCWire *wire = fdcpseudos[i]->wire;
1236 TVector3 pos(fdcpseudos[i]->xy.X(),
1237 fdcpseudos[i]->xy.Y(), wire->
origin.Z());
1238 gsetp.
points.push_back(pos);
1240 graphics.push_back(gsetp);
1247 vector<const DTOFPoint *>tofpoints;
1248 loop->Get(tofpoints);
1250 for(
unsigned int i=0; i<tofpoints.size(); i++){
1252 TVector3 pos(hit->
pos.x(),hit->
pos.y(),hit->
pos.z());
1253 gset.
points.push_back(pos);
1255 graphics.push_back(gset);
1260 vector<const DMCTrackHit*> mctrackhits;
1261 loop->Get(mctrackhits);
1263 for(
unsigned int i=0; i<mctrackhits.size(); i++){
1267 TVector3 pos(hit->
r*cos(hit->
phi), hit->
r*
sin(hit->
phi), hit->
z);
1268 gset.
points.push_back(pos);
1270 graphics.push_back(gset);
1275 vector<const DMCTrackHit*> mctrackhits;
1276 loop->Get(mctrackhits);
1278 for(
unsigned int i=0; i<mctrackhits.size(); i++){
1282 TVector3 pos(hit->
r*cos(hit->
phi), hit->
r*
sin(hit->
phi), hit->
z);
1283 gset.
points.push_back(pos);
1285 TMarker *m =
new TMarker(pos.X(), pos.Y(), 2);
1286 graphics_xyA.push_back(m);
1293 vector<const DFCALGeometry*> fcalgeometries;
1294 vector<const DFCALHit*> mcfcalhits;
1295 loop->Get(fcalgeometries);
1296 loop->Get(mcfcalhits,
"TRUTH");
1297 if(fcalgeometries.size()>0){
1301 for(
unsigned int i=0; i<mcfcalhits.size(); i++){
1302 const DFCALHit *hit = mcfcalhits[i];
1305 TVector3 pos(pos_face.X(), pos_face.Y(),
FCAL_Zmin);
1306 gset.
points.push_back(pos);
1308 TMarker *m =
new TMarker(pos.X(), pos.Y(), 2);
1311 graphics_xyB.push_back(m);
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);
1324 vector<const DNeutralParticle*> neutrals;
1325 loop->Get(neutrals);
1329 for(
unsigned int i=0; i<neutrals.size(); i++){
1330 auto locNeutralShower = neutrals[i]->dNeutralShower;
1333 TVector3 pos( locNeutralShower->dSpacetimeVertex.X(),
1334 locNeutralShower->dSpacetimeVertex.Y(),
1335 locNeutralShower->dSpacetimeVertex.Z());
1336 gset.
points.push_back(pos);
1338 double dist2 = 2.0 + 5.0*locNeutralShower->dEnergy;
1339 TEllipse *
e =
new TEllipse(pos.X(), pos.Y(), dist2, dist2);
1340 e->SetLineColor(kGreen);
1343 graphics_xyA.push_back(e);
1351 vector<const DNeutralParticle*> neutrals;
1352 loop->Get(neutrals);
1355 for(
unsigned int i=0; i<neutrals.size(); i++){
1356 auto locNeutralShower = neutrals[i]->dNeutralShower;
1360 TVector3 pos( locNeutralShower->dSpacetimeVertex.X(),
1361 locNeutralShower->dSpacetimeVertex.Y(),
1362 locNeutralShower->dSpacetimeVertex.Z());
1363 gset.
points.push_back(pos);
1365 double dist2 = 2.0 + 10.0*locNeutralShower->dEnergy;
1366 TEllipse *
e =
new TEllipse(pos.X(), pos.Y(), dist2, dist2);
1367 e->SetLineColor(kGreen);
1370 graphics_xyB.push_back(e);
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);
1390 vector<const DChargedTrack*> ctracks;
1392 for(
unsigned int i=0; i<ctracks.size(); i++){
1394 vector<const DNeutralShower*> locNeutralShowers;
1395 locCTrack->GetT(locNeutralShowers);
1397 if (!locNeutralShowers.size())
continue;
1410 gset.
points.push_back(tpos);
1411 graphics.push_back(gset);
1414 if(is_bcal)
continue;
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);
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);
1432 vector<const DMCThrown*> throwns;
1435 for(
unsigned int i=0; i<throwns.size(); i++){
1437 if(thrown->
charge()!=0.0)
continue;
1446 GetIntersectionWithCalorimeter(thrown, pos, who);
1451 TVector3 tpos(pos.X(),pos.Y(),pos.Z());
1452 gset.
points.push_back(tpos);
1457 double dist2 = 2.0 + 2.0*thrown->
energy();
1458 TEllipse *
e =
new TEllipse(pos.X(), pos.Y(), dist2, dist2);
1459 e->SetLineColor(kSpring);
1462 graphics_xyB.push_back(e);
1464 graphics.push_back(gset);
1469 vector<const DMCThrown*> throwns;
1472 for(
unsigned int i=0; i<throwns.size(); i++){
1474 if(thrown->
charge()==0.0)
continue;
1483 GetIntersectionWithCalorimeter(thrown, pos, who);
1490 TVector3 tpos(pos.X(),pos.Y(),pos.Z());
1491 gset.
points.push_back(tpos);
1492 graphics.push_back(gset);
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);
1500 graphics_xyB.push_back(e);
1512 vector<const DTrackTimeBased*> tracks;
1515 for(
unsigned int i=0; i<tracks.size(); i++){
1521 GetIntersectionWithCalorimeter(track, pos, who);
1528 TVector3 tpos(pos.X(),pos.Y(),pos.Z());
1529 gset.
points.push_back(tpos);
1530 graphics.push_back(gset);
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);
1540 graphics_xyB.push_back(e);
1546 vector<const DCCALShower*> clusters;
1547 loop->Get(clusters);
1548 for(
auto cluster : clusters){
1550 double E = cluster->E/1000.0;
1551 double dist2 = 1.0 + 0.5*E;
1552 DVector3 pos(cluster->x,cluster->y,cluster->z);
1554 TEllipse *
e =
new TEllipse(pos.X(), pos.Y(), dist2, dist2);
1555 e->SetLineColor(kGreen);
1558 graphics_xyB.push_back(e);
1560 TEllipse *e1 =
new TEllipse(pos.Z(), pos.X(), dist2, dist2);
1561 e1->SetLineColor(kGreen);
1562 e1->SetFillStyle(0);
1563 e1->SetLineWidth(2);
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);
1578 vector<const DMCTrajectoryPoint*> mctrajectorypoints;
1579 loop->Get(mctrajectorypoints);
1586 TVector3 last_point;
1589 for(
unsigned int i=0; i<mctrajectorypoints.size(); i++){
1619 TVector3 v(pt->
x, pt->
y, pt->
z);
1623 if(pt->
track!=last_track || pt->
part!=last_part){
1627 gset.
color = kOrange;
1631 gset.
color = kRed+2;
1636 gset.
color = kBlue+1;
1639 gset.
color = kGreen+2;
1642 gset.
color = kBlack;
1646 gset.
color = kBlack;
1648 graphics.push_back(gset);
1653 gset.
points.push_back(v);
1655 last_track = pt->
track;
1656 last_part = pt->
part;
1662 gset.
color = kOrange;
1666 gset.
color = kRed+2;
1671 gset.
color = kBlue+1;
1674 gset.
color = kGreen+2;
1677 gset.
color = kBlack;
1681 gset.
color = kBlack;
1683 graphics.push_back(gset);
1688 vector<const DTrackCandidate*> trackcandidates;
1690 for(
unsigned int i=0; i<trackcandidates.size(); i++){
1695 AddKinematicDataTrack(trackcandidates[i], color, size);
1701 vector<const DTrackWireBased*> wiretracks;
1703 for(
unsigned int i=0; i<wiretracks.size(); i++){
1704 AddKinematicDataTrack(wiretracks[i], (wiretracks[i]->charge()>0.0 ? kBlue:kRed)+2, 1.25);
1710 vector<const DTrackTimeBased*> timetracks;
1712 for(
unsigned int i=0; i<timetracks.size(); i++){
1713 AddKinematicDataTrack(timetracks[i], (timetracks[i]->charge()>0.0 ? kBlue:kRed)+0, 1.00);
1719 vector<const DChargedTrack*> chargedtracks;
1721 for(
unsigned int i=0; i<chargedtracks.size(); i++){
1722 int color=kViolet-3;
1724 if (chargedtracks[i]->Get_Charge() > 0) color=kMagenta;
1726 if (chargedtracks[i]->Get_BestFOM()->mass() > 0.9) size=2.5;
1727 AddKinematicDataTrack(chargedtracks[i]->Get_BestFOM(),color,size);
1733 vector<const DNeutralParticle*> photons;
1736 for(
unsigned int i=0; i<photons.size(); i++){
1738 auto locNeutralShower = photons[i]->dNeutralShower;
1740 if(locDetSys==
SYS_FCAL)color = kOrange;
1741 if(locDetSys==
SYS_BCAL)color = kYellow+2;
1743 AddKinematicDataTrack(photons[i]->Get_BestFOM(), color, 1.00);
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);
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++){
1794 const DBCALHit* hit = locBcalHits[k];
1795 float idxY = (float)hit->
layer-1;
1796 float idxX = (
float) (hit->
sector-1 + (hit->
module-1)*4);
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);
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);
1832 LayerLegend =
new TLegend(0.91,0.7,0.99,0.99);
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);
1844 sprintf(name,
"layer %i",layer+1);
1845 LayerLegend->AddEntry(BCALPointZphiLayer[layer],name);
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++){
1856 float p = TMath::Sqrt(part->
px*part->
px + part->
py*part->
py + part->
pz*part->
pz);
1860 phi = TMath::ATan(TMath::Abs(part->
y/part->
x));
1864 phi = 3.1415926 - phi;
1870 phi = 3.1415926*2. - phi;
1874 phi = phi*180./3.1415926;
1876 BCALParticles->Fill(phi,0.5,p);
1879 TText *t =
new TText(phi,1.01,l);
1880 t->SetTextSize(0.08);
1882 t->SetTextAlign(21);
1883 BCALPLables.push_back(t);
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);
1901 BCALHitCanvas->Update();
1917 vector<const DMCThrown*> throwns;
1918 if(loop)loop->Get(throwns);
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]);
1929 if(loop)loop->Get(TrksWireBased);
1930 if(loop)loop->Get(TrksTimeBased);
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());
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]);
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]);
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]);
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());
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(
"--------");
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(
"--------");
1979 for(
unsigned int i=0; i<throwns.size(); i++){
1981 if(trk->
type==0)
continue;
1982 int row = thrownlabs[
"trk"].size()-(ii++)-1;
1985 stringstream trkno, type, p, theta, phi, z;
1986 trkno<<setprecision(4)<<i+1;
1987 thrownlabs[
"trk"][row]->SetText(trkno.str().c_str());
1991 p<<setprecision(4)<<trk->
momentum().Mag();
1992 thrownlabs[
"p"][row]->SetText(p.str().c_str());
1994 theta<<setprecision(4)<<trk->
momentum().Theta()*TMath::RadToDeg();
1995 thrownlabs[
"theta"][row]->SetText(theta.str().c_str());
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());
2002 z<<setprecision(4)<<trk->
position().Z();
2003 thrownlabs[
"z"][row]->SetText(z.str().c_str());
2007 for(
unsigned int i=0; i<trks.size(); i++){
2009 int row = reconlabs[
"trk"].size()-i-1;
2012 stringstream trkno, type, p, theta, phi, z, chisq_per_dof, Ndof,cand;
2014 trkno<<setprecision(4)<<i+1;
2015 reconlabs[
"trk"][row]->SetText(trkno.str().c_str());
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";
2024 if (fabs(trk->
charge())>1.
e-4){
2025 type<<(trk->
charge()>0 ?
"+":
"-");
2027 reconlabs[
"type"][row]->SetText(type.str().c_str());
2029 p<<setprecision(4)<<trk->
momentum().Mag();
2030 reconlabs[
"p"][row]->SetText(p.str().c_str());
2032 theta<<setprecision(4)<<trk->
momentum().Theta()*TMath::RadToDeg();
2033 reconlabs[
"theta"][row]->SetText(theta.str().c_str());
2035 phi<<setprecision(4)<<trk->
momentum().Phi()*TMath::RadToDeg();
2036 reconlabs[
"phi"][row]->SetText(phi.str().c_str());
2038 z<<setprecision(4)<<trk->
position().Z();
2039 reconlabs[
"z"][row]->SetText(z.str().c_str());
2049 chisq_per_dof<<setprecision(4)<<timetrack->
chisq/timetrack->
Ndof;
2050 Ndof<<timetrack->
Ndof;
2051 fom << timetrack->
FOM;
2053 chisq_per_dof<<setprecision(4)<<track->
chisq/track->
Ndof;
2056 }
else if (candidate){
2057 chisq_per_dof<<setprecision(4)<<candidate->
chisq/candidate->
Ndof;
2058 Ndof<<candidate->
Ndof;
2061 else if (chargedtrack){
2063 chisq_per_dof<<setprecision(4)<<locTrackTimeBased->
chisq/locTrackTimeBased->Ndof;
2064 Ndof<<locTrackTimeBased->Ndof;
2065 fom << locTrackTimeBased->FOM;
2068 chisq_per_dof <<
"--------";
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());
2086 reconlabs[
"cand"][row]->SetText(cand.str().c_str());
2090 if( fulllistmf ) fulllistmf->UpdateTrackLabels(throwns, trks);
2092 debugermf->SetTrackCandidates(TrksCand);
2093 debugermf->SetTrackWireBased(TrksWireBased);
2094 debugermf->SetTrackTimeBased(TrksTimeBased);
2095 debugermf->UpdateTrackLabels();
2110 if(MATERIAL_MAP_MODEL==
"DRootGeom"){
2113 }
else if(MATERIAL_MAP_MODEL==
"DGeometry"){
2116 }
else if(MATERIAL_MAP_MODEL!=
"NONE"){
2117 _DBG_<<
"WARNING: Invalid value for TRKFIT:MATERIAL_MAP_MODEL (=\""<<MATERIAL_MAP_MODEL<<
"\")"<<endl;
2128 gset.
points.push_back(tpoint);
2132 graphics.push_back(gset);
2146 if(MATERIAL_MAP_MODEL==
"DRootGeom"){
2149 }
else if(MATERIAL_MAP_MODEL==
"DGeometry"){
2152 }
else if(MATERIAL_MAP_MODEL!=
"NONE"){
2153 _DBG_<<
"WARNING: Invalid value for TRKFIT:MATERIAL_MAP_MODEL (=\""<<MATERIAL_MAP_MODEL<<
"\")"<<endl;
2161 double s_fcal = 1.0E6;
2167 if(pos_fcal.Perp()<
FCAL_Rmin || pos_fcal.Perp()>
FCAL_Rmax || !isfinite(pos_fcal.Z()))s_fcal = 1.0E6;
2171 double s_bcal = 1.0E6;
2175 if(s_fcal>10000.0 && s_bcal>10000.0){
2178 pos.SetXYZ(0.0,0.0,0.0);
2179 }
else if(s_fcal<s_bcal){
2195 vector<JEventLoop*> loops = app->GetJEventLoops();
2197 vector<string> facnames;
2198 loops[0]->GetFactoryNames(facnames);
2207 vector<JEventLoop*> loops = app->GetJEventLoops();
2209 factories = loops[0]->GetFactories();
2218 vector<JEventLoop*> loops = app->GetJEventLoops();
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();
2230 JFactory_base *fac = loops[0]->GetFactory(factory, tag.c_str());
2234 if(loops[0]->GetJEvent().GetJEventSource() == NULL)
return 0;
2236 return fac==NULL ? 0:(
unsigned int)fac->GetNrows();
2253 vector<JEventLoop*> loops = app->GetJEventLoops();
2254 if(loops.size()==0)
return;
2255 JEventLoop* &loop = loops[0];
2260 double mass = 0.13957018;
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();
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();
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();
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();
2311 if(dataname==
"DMCThrown"){
2312 vector<const DMCThrown*> tracks;
2313 loop->Get(tracks, tag.c_str());
2314 if(index>=tracks.size())
return;
2319 tracks[
index]->Get(cdchits);
2320 mass = tracks[
index]->mass();
2321 _DBG_<<
"mass="<<mass<<endl;
2325 if(q==0.0 || mom.Mag()<0.01)
return;
2333 if(MATERIAL_MAP_MODEL==
"DRootGeom"){
2336 }
else if(MATERIAL_MAP_MODEL==
"DGeometry"){
2339 }
else if(MATERIAL_MAP_MODEL!=
"NONE"){
2340 _DBG_<<
"WARNING: Invalid value for TRKFIT:MATERIAL_MAP_MODEL (=\""<<MATERIAL_MAP_MODEL<<
"\")"<<endl;
2342 rt->
Swim(pos, mom, q);
2357 vector<JEventLoop*> loops = app->GetJEventLoops();
2358 if(loops.size()==0)
return;
2359 JEventLoop* &loop = loops[0];
2362 vector<const DCDCTrackHit*> 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);
2372 vector<const DFDCPseudo*> 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);
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);
void SetNTrWireBased(Int_t d)
hdv_fulllistframe * GetFullListFrame(void)
void SetReconstructedFactories(vector< string > &facnames)
void UpdateBcalDisp(void)
float chisq
Chi-squared for the track (not chisq/dof!)
TCanvas * GetBcalDispFrame(void)
float chisq
Chi-squared for the track (not chisq/dof!)
void SetTrig(char *trigstring)
bool GetCheckButton(string who)
double energy(void) const
sprintf(text,"Post KinFit Cut")
vector< TVector3 > points
TPolyLine * GetCCALPolyLine(int row, int col)
DMagneticFieldMap * GetBfield(unsigned int run_number=1)
void SetEvent(ULong64_t id)
void SetMass(double mass)
const DVector3 & position(void) const
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)
void SetNTrTimeBased(Int_t d)
void AddKinematicDataTrack(const DKinematicData *kd, int color, double size)
static char index(char c)
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
void SetDGeometry(const DGeometry *geom)
DetectorSystem_t system
particle type
DetectorSystem_t dDetectorSystem
DLorentzVector dSpacetimeVertex
DRootGeom * GetRootGeom(unsigned int run_number)
map< string, vector< TGLabel * > > & GetThrownLabels(void)
float z
coordinates of hit in cm and rad
TPolyLine * GetBCALPolyLine(int mod, int layer, int sector)
bool DMCTrajectoryPoint_track_cmp(const DMCTrajectoryPoint *a, const DMCTrajectoryPoint *b)
TPolyLine * GetFCALPolyLine(int channel)
DGeometry * GetDGeometry(unsigned int run_number)
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
void SetDRootGeom(const DRootGeom *RootGeom)
int Ndof
Number of degrees of freedom in the fit.
double charge(void) const
void GetReconFactory(string &name, string &tag)
void GetFactories(vector< JFactory_base * > &factories)
void GetAllWireHits(vector< pair< const DCoordinateSystem *, double > > &allhits)
void GetIntersectionWithCalorimeter(const DKinematicData *kd, DVector3 &pos, DetectorSystem_t &who)
<A href="index.html#legend"> <IMG src="CORE.png" width="100"> </A>
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)
hdv_debugerframe * GetDebugerFrame(void)
float chisq
Chi-squared for the track (not chisq/dof!)
const DVector3 & momentum(void) const
bool GetFDCWires(vector< vector< DFDCWire * > > &fdcwires) const
oid_t candidateid
which DTrackCandidate this came from
File: DTOFHit.h Created: Tue Jan 18 16:15:26 EST 2011 Creator: B. Zihlmann Purpose: Container class t...
void GetFactoryNames(vector< string > &facnames)
int Ndof
Number of degrees of freedom in the fit.
int type
GEANT particle ID.
printf("string=%s", string)
jerror_t evnt(JEventLoop *eventLoop, uint64_t eventnumber)
Called every event.
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
class DFDCHit: definition for a basic FDC hit data type.
unsigned int GetNrows(const string &factory, string tag)