Merge pull request #65 from JeffersonLab/hdview2_primex

Hdview2 primex

src/libraries/CCAL/DCCALCluster.h

@@ -18,7 +18,7 @@ using namespace std;
 using namespace jana;
 
 #define CCAL_USER_HITS_MAX  1000
-#define MAX_SHOWER_RADIUS   5
+//#define MAX_SHOWER_RADIUS   5   <-- This appears to have been copied from DFCALCluster.h and then not used  2018-12-10  DL
 
 class DCCALCluster : public JObject {
    public:

src/libraries/HDGEOMETRY/DGeometry.cc

@@ -1666,6 +1666,24 @@ bool DGeometry::GetBCALPhiShift(float &bcal_phi_shift) const
    }
 }
 
+//---------------------------------
+// GetCCALZ
+//---------------------------------
+bool DGeometry::GetCCALZ(double &z_ccal) const
+{
+   vector<double> ComptonEMcalpos;
+   bool good = Get("//section/composition/posXYZ[@volume='ComptonEMcal']/@X_Y_Z", ComptonEMcalpos);
+
+   if(!good){
+      _DBG_<<"Unable to retrieve ComptonEMcal position."<<endl;
+      z_ccal=876.106; // from some version of HDDS that may be out of date.  2018-12-10 DL
+      return false;
+   }else{
+	   z_ccal = ComptonEMcalpos[2];
+      return true;
+   }
+}
+
 //---------------------------------
 // GetFCALZ
 //---------------------------------

src/libraries/HDGEOMETRY/DGeometry.h

@@ -149,6 +149,8 @@ class DGeometry{
       bool GetBCALDepth(float &bcal_depth) const; ///< depth (or height) of BCAL module in cm
       bool GetBCALPhiShift(float &bcal_phi_shift) const; ///< phi angle in degrees that first BCAL module is shifted from being centered at ph=0.0
 
+      bool GetCCALZ(double &z_ccal) const; /// z-location of front face of CCAL in cm
+
       bool GetFCALZ(double &z_fcal) const; ///< z-location of front face of FCAL in cm
       bool GetDIRCZ(double &z_dirc) const; ///< z-location of DIRC in cm
       bool GetTOFZ(vector<double> &z_tof) const; ///< z-location of front face of each of TOF in cm

src/plugins/monitoring/occupancy_online/CCAL_occupancy.C

@@ -8,8 +8,8 @@
 // The following are special comments used by RootSpy to know
 // which histograms to fetch for the macro.
 //
-// hnamepath: /occupancy/fcal_occ
-// hnamepath: /occupancy/fcal_num_events
+// hnamepath: /occupancy/ccal_occ
+// hnamepath: /occupancy/ccal_num_events
 //
 // e-mail: davidl@jlab.org
 // e-mail: staylor@jlab.org

src/programs/Analysis/hdview2/MyProcessor.cc

@@ -38,6 +38,8 @@ using namespace std;
 #include "JANA/JGeometry.h"
 #include "TRACKING/DMCTrajectoryPoint.h"
 #include "FCAL/DFCALHit.h"
+#include "CCAL/DCCALHit.h"
+#include "CCAL/DCCALCluster.h"
 #include "TOF/DTOFGeometry.h"
 #include "TOF/DTOFHit.h"
 #include "TOF/DTOFTDCDigiHit.h" 
@@ -668,6 +670,61 @@ void MyProcessor::FillGraphics(void)
 	    poly->SetFillColor(TColor::GetColor(r,g,b));
 	  }
 	}
+
+	// CCAL hits
+	if(hdvmf->GetCheckButton("ccal")){
+	  vector<const DCCALHit*> ccalhits;
+	  loop->Get(ccalhits);
+
+	  for(unsigned int i=0; i<ccalhits.size(); i++){
+	    const DCCALHit *hit = ccalhits[i];
+	    TPolyLine *poly = hdvmf->GetCCALPolyLine(hit->row, hit->column);
+	    if(!poly)continue;
+
+	    // The aim is to have a log scale in energy (see BCAL)
+	    double E = 1000*hit->E;      // Change Energy to MeV
+	    E /= 1000.0; // CCAL is currently not calibrated and appears to be at least 1E3 too large  2018-12-10 DL
+	    if(E<0.0) continue;
+	    double logE = log10(E);
+
+	    float r,g,b;
+	    if (logE<0){
+	    	r = 1.;
+	    	g = 1.;
+	    	b = 1.;
+	    } else {
+	    	if (logE<1){
+	    		r = 1.;
+	    		g = 1.;
+	    		b = 1.-logE;
+	    	} else {
+	    		if (logE<2){
+	    			r = 1.;
+	    			g = 1.-(logE-1);
+	    			b = 0.;
+	    		} else {
+	    			if (logE<3){
+	    				r = 1.;
+	    				g = 0.;
+	    				b = 1.-(logE-2);
+	    			} else {
+	    				if (logE<4){
+	    					r = 1.-(logE-3);
+	    					g = 0.;
+	    					b = 1.;
+	    				} else {
+	    					r = 0;
+	    					g = 0;
+	    					b = 0;
+	    				}
+	    			}
+	    		}
+	    	}
+	    }
+	    poly->SetFillColor(TColor::GetColor(r,g,b));
+	  }
+	}
+
 	// TOF hits
 	if(hdvmf->GetCheckButton("tof")){
 
@@ -1484,6 +1541,38 @@ void MyProcessor::FillGraphics(void)
 		}
 	}
 
+	// CCAL reconstructed clusters
+	if(hdvmf->GetCheckButton("recon_photons_ccal")){
+		vector<const DCCALCluster*> clusters;
+		loop->Get(clusters);
+		for(auto cluster : clusters){
+
+			double E = cluster->getEnergy()/1000.0; // divide by 1000 since energy does not seem to be calibrated to GeV at the moment.  2018-12-10 DL
+			double dist2 = 1.0 + 0.5*E;
+			DVector3 pos = cluster->getCentroid();
+
+			TEllipse *e = new TEllipse(pos.X(), pos.Y(), dist2, dist2);
+			e->SetLineColor(kGreen);
+			e->SetFillStyle(0);
+			e->SetLineWidth(2);
+			graphics_xyB.push_back(e);
+
+			TEllipse *e1 = new TEllipse(pos.Z(), pos.X(), dist2, dist2);
+			e1->SetLineColor(kGreen);
+			e1->SetFillStyle(0);
+			e1->SetLineWidth(2);
+
+			graphics_xz.push_back(e1);
+			TEllipse *e2 = new TEllipse(pos.Z(), pos.Y(), dist2, dist2);
+			e2->SetLineColor(kGreen);
+			e2->SetFillStyle(0);
+			e2->SetLineWidth(2);
+			graphics_yz.push_back(e2);
+
+		}
+		//graphics.push_back(gset);
+	}
+
 	// DMCTrajectoryPoints
 	if(hdvmf->GetCheckButton("trajectories")){
 		vector<const DMCTrajectoryPoint*> mctrajectorypoints;

src/programs/Analysis/hdview2/hdv_mainframe.cc

@@ -65,6 +65,10 @@ static float FCAL_Zlen = 45.0;
 static float FCAL_Zmin = 622.8;
 static float FCAL_Rmin = 6.0;
 static float FCAL_Rmax = 212.0/2.0;
+static float CCAL_Zlen = 18.0;
+static float CCAL_Zmin = 876.106;
+static float CCAL_Rmin = 2.0;
+static float CCAL_Rmax = 24.0;
 static float CDC_Rmin = 9.0;
 static float CDC_Rmax = 59.0;
 static float CDC_Zlen = 150.0;
@@ -371,6 +375,7 @@ hdv_mainframe::hdv_mainframe(const TGWindow *p, UInt_t w, UInt_t h):TGMainFrame(
   checkbuttons["toftruth"]	= new TGCheckButton(hitdrawopts,	"TOFTruth");
   checkbuttons["fcal"]		= new TGCheckButton(hitdrawopts,	"FCAL");
   checkbuttons["bcal"]		= new TGCheckButton(hitdrawopts,	"BCAL");
+  checkbuttons["ccal"]		= new TGCheckButton(hitdrawopts,	"CCAL");
 
   hitdrawopts->AddFrame(checkbuttons["cdc"], lhints);
   hitdrawopts->AddFrame(checkbuttons["cdcdrift"], lhints);
@@ -382,6 +387,7 @@ hdv_mainframe::hdv_mainframe(const TGWindow *p, UInt_t w, UInt_t h):TGMainFrame(
   hitdrawopts->AddFrame(checkbuttons["toftruth"], lhints);
   hitdrawopts->AddFrame(checkbuttons["fcal"], lhints);
   hitdrawopts->AddFrame(checkbuttons["bcal"], lhints);
+  hitdrawopts->AddFrame(checkbuttons["ccal"], lhints);
 
   TGTextButton *moreOptions	= new TGTextButton(hitdrawopts,	"More options");
   hitdrawopts->AddFrame(moreOptions, lhints);
@@ -568,6 +574,7 @@ hdv_mainframe::hdv_mainframe(const TGWindow *p, UInt_t w, UInt_t h):TGMainFrame(
   checkbuttons["bcaltruth"]->Connect("Clicked()","hdv_mainframe", this, "DoMyRedraw()");
   checkbuttons["fcal"]->Connect("Clicked()","hdv_mainframe", this, "DoMyRedraw()");
   checkbuttons["fcaltruth"]->Connect("Clicked()","hdv_mainframe", this, "DoMyRedraw()");
+  checkbuttons["ccal"]->Connect("Clicked()","hdv_mainframe", this, "DoMyRedraw()");
   checkbuttons["trajectories"]->Connect("Clicked()","hdv_mainframe", this, "DoMyRedraw()");
 
   for (Int_t n=1;n<debugermf->GetNTrCand();n++){
@@ -743,13 +750,13 @@ void hdv_mainframe::SetRange(void)
 		// define range in each direction in cm
 		double x_width = 350.0/zoom_factor;
 		double y_width = x_width;
-		double z_width = 2.0*x_width;
+		double z_width = 2.70*x_width;
 		double xlo = x0 - x_width/2.0;
 		double xhi = x0 + x_width/2.0;
 		double ylo = y0 - y_width/2.0;
 		double yhi = y0 + y_width/2.0;
-		double zlo = z0 - z_width/2.0;
-		double zhi = z0 + z_width/2.0;
+		double zlo = z0 - z_width*0.40;
+		double zhi = z0 + z_width*0.60;
 
 		sideviewA->GetCanvas()->cd();
 		sideviewA->GetCanvas()->Range(zlo, xlo, zhi, xhi);
@@ -766,11 +773,11 @@ void hdv_mainframe::SetRange(void)
 		endviewA->GetCanvas()->cd();
 		endviewA->GetCanvas()->Range(xlo, ylo, xhi, yhi);
 		endviewB->GetCanvas()->cd();
-		//endviewB->GetCanvas()->Range(xlo*1.3, ylo*1.3, xhi*1.3, yhi*1.3);
-		endviewB->GetCanvas()->Range(-158, -158, 158, 158);
+		endviewB->GetCanvas()->Range(xlo*1.3, ylo*1.3, xhi*1.3, yhi*1.3);
+		//endviewB->GetCanvas()->Range(-158, -158, 158, 158);
 		endviewAmf->SetRange(xlo, ylo, xhi, yhi);
-		//endviewBmf->SetRange(xlo*1.3, ylo*1.3, xhi*1.3, yhi*1.3);
-		endviewBmf->SetRange(-158, -158, 158, 158);
+		endviewBmf->SetRange(xlo*1.3, ylo*1.3, xhi*1.3, yhi*1.3);
+		//endviewBmf->SetRange(-158, -158, 158, 158);
 
 	}else{
 		// define range in each direction in cm, radians
@@ -1368,6 +1375,14 @@ void hdv_mainframe::DrawDetectorsXY(void)
 		graphics_sideA.push_back(fcal1);
 		graphics_sideA.push_back(fcal2);
 
+		// ----- CCAL ------
+		TBox *ccal1 = new TBox(CCAL_Zmin,  CCAL_Rmin, CCAL_Zmin+CCAL_Zlen,  CCAL_Rmax);
+		TBox *ccal2 = new TBox(CCAL_Zmin, -CCAL_Rmin, CCAL_Zmin+CCAL_Zlen, -CCAL_Rmax);
+		ccal1->SetFillColor(42);
+		ccal2->SetFillColor(42);
+		graphics_sideA.push_back(ccal1);
+		graphics_sideA.push_back(ccal2);
+
 		// ------ scale ------
 		DrawScale(sideviewA->GetCanvas(), graphics_sideA);
 	}
@@ -1609,7 +1624,7 @@ void hdv_mainframe::DrawDetectorsXY(void)
 
 		// Set up 4 2-D vectors that point from the center of a block to its
 		// corners. This makes it easier to represent each corner as a vector
-		// in lab corrdinate whch we can extract r, phi from.
+		// in lab coordinate which we can extract r, phi from.
 		double blocksize = fcalgeom->blockSize();
 		DVector2 shift[4];
 		shift[0].Set(-blocksize/2, -blocksize/2);  // these are ordered such that they
@@ -1623,14 +1638,16 @@ void hdv_mainframe::DrawDetectorsXY(void)
 				int row = fcalgeom->row(chan);
 				int col = fcalgeom->column(chan);
 				if(!fcalgeom->isBlockActive(row, col))continue;
-				double x[4], y[4];
+				double x[5], y[5];
 				for(int i=0; i<4; i++){
 					DVector2 pos = shift[i] + fcalgeom->positionOnFace(chan);
 					x[i] = pos.X();
 					y[i] = pos.Y();
 				}
+				x[4] = x[0];
+				y[4] = y[0];
 
-				TPolyLine *poly = new TPolyLine(4, x, y);
+				TPolyLine *poly = new TPolyLine(5, x, y);
 				poly->SetFillColor(0);
 				poly->SetLineColor(kBlack);
 				graphics_endB.push_back(poly);
@@ -1639,6 +1656,60 @@ void hdv_mainframe::DrawDetectorsXY(void)
 			}
 		}
 
+		// ----- CCAL ------
+		// Get list of blocks. Loop over all getting x,y coordinates of corners for all active ones.
+
+		// Set up 4 2-D vectors that point from the center of a block to its
+		// corners. This makes it easier to represent each corner as a vector
+		// in lab coordinate which we can extract r, phi from.
+		blocksize = 2.0;
+		shift[0].Set(-blocksize/2, -blocksize/2);  // these are ordered such that they
+		shift[1].Set(-blocksize/2, +blocksize/2);  // go in a clockwise manner. This
+		shift[2].Set(+blocksize/2, +blocksize/2);  // ensures the r/phi coordinates also
+		shift[3].Set(+blocksize/2, -blocksize/2);  // define a single enclosed space
+		ccalblocks.clear();
+
+		if(GetCheckButton("ccal")){
+			for(int irow=0; irow<12; irow++){
+				for(int icol=0; icol<12; icol++){
+					if( (irow==5 || irow==6) && (icol==5 || icol==6) ) continue;
+
+					double center_x = (-6.0 + (icol+0.5))*2.0;
+					double center_y = (-6.0 + (irow+0.5))*2.0;
+					DVector2 mypos(center_x,center_y);
+
+					double x[5], y[5];
+					for(int i=0; i<4; i++){
+						DVector2 pos = shift[i] + mypos;
+						x[i] = pos.X();
+						y[i] = pos.Y();
+					}
+					x[4] = x[0];
+					y[4] = y[0];
+
+					TPolyLine *poly = new TPolyLine(5, x, y);
+					poly->SetFillColor(0);
+					poly->SetLineColor(kBlack);
+					graphics_endB.push_back(poly);
+
+					int channel = icol + 12*(irow);
+					ccalblocks[channel] = poly; // record so we can set the color later
+				}
+			}
+
+			// If also drawing FCAL, draw an extra outline of the FCAL beam hole
+			if(GetCheckButton("fcal")){
+				double x[5] = {-6.0, -6.0, 6.0, 6.0, -6.0};
+				double y[5] = {-6.0,  6.0, 6.0,-6.0, -6.0};
+				TPolyLine *poly = new TPolyLine(5, x, y);
+				poly->SetFillColor(0);
+				poly->SetLineColor(kBlack);
+				poly->SetLineWidth(3);
+				poly->SetLineStyle(2);
+				graphics_endB.push_back(poly);
+			}
+		}
+
 		// ------- TOF ---------// 
 		tofblocks.clear();
 		if(GetCheckButton("tof")){
@@ -2419,6 +2490,17 @@ TPolyLine* hdv_mainframe::GetFCALPolyLine(float x, float y)
 	return GetFCALPolyLine(fcalgeom->channel(row, column));
 }
 
+//-------------------
+// GetCCALPolyLine
+//-------------------
+TPolyLine* hdv_mainframe::GetCCALPolyLine(int row, int col)
+{
+	int channel = col + 12*(row);
+	map<int, TPolyLine*>::iterator iter = ccalblocks.find(channel);
+	if(iter==ccalblocks.end())return NULL;
+	return iter->second;
+}
+
 //------------------- 
 // GetTOFPolyLine
 //-------------------

src/programs/Analysis/hdview2/hdv_mainframe.h

@@ -139,6 +139,7 @@ class hdv_mainframe:public TGMainFrame {
 		void GetReconFactory(string &name, string &tag);
 		TPolyLine* GetFCALPolyLine(int channel);
 		TPolyLine* GetFCALPolyLine(float x, float y);
+		TPolyLine* GetCCALPolyLine(int row, int col);
 		TPolyLine* GetBCALPolyLine(int mod, int layer, int sector);
 		TPolyLine* GetTOFPolyLine(int translate_side, int tof_ch);
 
@@ -201,6 +202,7 @@ class hdv_mainframe:public TGMainFrame {
 		map<string, TGCheckButton*> checkbuttons;
 		map<int, TPolyLine*> fcalblocks;
 		map<int, TPolyLine*> bcalblocks;
+		map<int, TPolyLine*> ccalblocks;
 		map<int, map<int, TPolyLine*> > tofblocks;
 
 		TTimer *timer;

src/programs/Analysis/hdview2/hdv_optionsframe.cc

@@ -77,6 +77,7 @@ hdv_optionsframe::hdv_optionsframe(hdv_mainframe *hdvmf, const TGWindow *p, UInt
 	checkbuttons["recon_charged_fcal"]			= new TGCheckButton(hitdrawopts,	"Draw reconstructed charged track projections on FCAL");
 	checkbuttons["thrown_photons_fcal"]			= new TGCheckButton(hitdrawopts,	"Draw thrown photon projections on FCAL");
 	checkbuttons["recon_photons_fcal"]			= new TGCheckButton(hitdrawopts,	"Draw reconstructed photon projections on FCAL");
+	checkbuttons["recon_photons_ccal"]			= new TGCheckButton(hitdrawopts,	"Draw reconstructed photon projections on CCAL");
 	TGHorizontal3DLine *separator2				= new TGHorizontal3DLine(hitdrawopts, 250);
 	checkbuttons["recon_photons_track_match"]	= new TGCheckButton(hitdrawopts,	"Draw reconstructed photons matched to charged tracks");
 	TGHorizontal3DLine *separator3				= new TGHorizontal3DLine(hitdrawopts, 250);
@@ -105,6 +106,8 @@ hdv_optionsframe::hdv_optionsframe(hdv_mainframe *hdvmf, const TGWindow *p, UInt
 	hitdrawopts->AddFrame(checkbuttons["thrown_photons_fcal"], lhints);
 	hitdrawopts->AddFrame(checkbuttons["recon_photons_fcal"], lhints);
 	hitdrawopts->AddFrame(separator2, chints);
+	hitdrawopts->AddFrame(checkbuttons["recon_photons_ccal"], lhints);
+	hitdrawopts->AddFrame(separator2, chints);
 	hitdrawopts->AddFrame(checkbuttons["recon_photons_track_match"], lhints);
 	hitdrawopts->AddFrame(separator3, chints);
 	hitdrawopts->AddFrame(checkbuttons["trajectories_photon"], lhints);