Merge pull request #245 from JeffersonLab/trd_fullfield

Options for parsing GEM SRS hits and updates to DIRC plugins

src/libraries/DAQ/DEVIOWorkerThread.cc

@@ -71,6 +71,7 @@ DEVIOWorkerThread::DEVIOWorkerThread(
 	PARSE_TRIGGER       = true;
 	PARSE_SSP           = true;
 	PARSE_GEMSRS        = true;
+        NSAMPLES_GEMSRS     = 9;
 
 	LINK_TRIGGERTIME    = true;
 }
@@ -2329,7 +2330,7 @@ void DEVIOWorkerThread::MakeDGEMSRSWindowRawData(DParsedEvent *pe, uint32_t roci
 	rawDataTS.clear();
 
 	Int_t fAPVHeaderLevel = 1500;
-	Int_t fNbOfTimeSamples = 21; // hard coded maximum number of time samples
+	Int_t fNbOfTimeSamples = NSAMPLES_GEMSRS; // hard coded maximum number of time samples
 
 	uint8_t NCH = 128;
 

src/libraries/DAQ/DEVIOWorkerThread.h

@@ -91,7 +91,8 @@ class DEVIOWorkerThread{
 		bool  PARSE_TRIGGER;
 		bool  PARSE_SSP;
 		bool  PARSE_GEMSRS;
-
+                int   NSAMPLES_GEMSRS;
+
 		bool  LINK_TRIGGERTIME;
 		bool  LINK_CONFIG;
 

src/libraries/DAQ/JEventSource_EVIOpp.cc

@@ -93,6 +93,7 @@ JEventSource_EVIOpp::JEventSource_EVIOpp(const char* source_name):JEventSource(s
 	PARSE_TRIGGER = true;
 	PARSE_SSP = true;
 	PARSE_GEMSRS = false;
+        NSAMPLES_GEMSRS = 9;
 	APPLY_TRANSLATION_TABLE = true;
 	IGNORE_EMPTY_BOR = false;
 	F250_EMULATION_MODE = kEmulationAuto;
@@ -133,6 +134,7 @@ JEventSource_EVIOpp::JEventSource_EVIOpp(const char* source_name):JEventSource(s
 	gPARMS->SetDefaultParameter("EVIO:PARSE_TRIGGER", PARSE_TRIGGER, "Set this to 0 to disable parsing of the built trigger bank from CODA (for benchmarking/debugging)");
 	gPARMS->SetDefaultParameter("EVIO:PARSE_SSP", PARSE_SSP, "Set this to 0 to disable parsing of the SSP (DIRC data) bank from CODA (for benchmarking/debugging)");
 	gPARMS->SetDefaultParameter("EVIO:PARSE_GEMSRS", PARSE_GEMSRS, "Set this to 0 to disable parsing of the SRS (GEM data) bank from CODA (for benchmarking/debugging)");
+        gPARMS->SetDefaultParameter("EVIO:NSAMPLES_GEMSRS", NSAMPLES_GEMSRS, "Set this to number of readout samples for SRS (GEM data) bank from CODA (for benchmarking/debugging)");
 	gPARMS->SetDefaultParameter("EVIO:APPLY_TRANSLATION_TABLE", APPLY_TRANSLATION_TABLE, "Apply the translation table to create DigiHits (you almost always want this on)");
 	gPARMS->SetDefaultParameter("EVIO:IGNORE_EMPTY_BOR", IGNORE_EMPTY_BOR, "Set to non-zero to continue processing data even if an empty BOR event is encountered.");
 	gPARMS->SetDefaultParameter("EVIO:TREAT_TRUNCATED_AS_ERROR", TREAT_TRUNCATED_AS_ERROR, "Set to non-zero to have a truncated EVIO file the JANA return code to non-zero indicating the program errored.");
@@ -218,6 +220,7 @@ JEventSource_EVIOpp::JEventSource_EVIOpp(const char* source_name):JEventSource(s
 		w->PARSE_TRIGGER       = PARSE_TRIGGER;
 		w->PARSE_SSP           = PARSE_SSP;
 		w->PARSE_GEMSRS        = PARSE_GEMSRS;
+                w->NSAMPLES_GEMSRS     = NSAMPLES_GEMSRS;
 		w->LINK_TRIGGERTIME    = LINK_TRIGGERTIME;
 		w->LINK_CONFIG         = LINK_CONFIG;
 		w->run_number_seed     = run_number_seed;

src/libraries/DAQ/JEventSource_EVIOpp.h

@@ -190,6 +190,7 @@ class JEventSource_EVIOpp: public jana::JEventSource{
 		bool     PARSE_TRIGGER;
 		bool     PARSE_SSP;
 		bool     PARSE_GEMSRS;
+                int      NSAMPLES_GEMSRS;
 		bool     APPLY_TRANSLATION_TABLE;
 		int      ET_STATION_NEVENTS;
 		bool     ET_STATION_CREATE_BLOCKING;

src/libraries/TRD/DGEMHit_factory.cc

@@ -81,37 +81,52 @@ jerror_t DGEMHit_factory::evnt(JEventLoop *loop, uint64_t eventnumber)
 
     vector<const DGEMDigiWindowRawData*> windowrawdata;
     eventLoop->Get(windowrawdata);
+    if(windowrawdata.size() == 0) 
+      return NOERROR;
 
-    // loop over GEM channels
+    // determine pedestals for each APV and time slice
+    const DGEMSRSWindowRawData* srswindow;
+    windowrawdata[windowrawdata.size()-1]->GetSingle(srswindow);
+    uint nSamples = srswindow->samples.size();
+    uint nAPV = srswindow->apv_id + 1;
+
+    // vector of each time slice for given APV
+    vector<double> pedestalAPV[nAPV]; 
+    for(uint i=0; i<nAPV; i++) pedestalAPV[i].resize(nSamples);
+
+    for (const auto& window : windowrawdata) {
+	const DGEMSRSWindowRawData* srswindow;
+	window->GetSingle(srswindow);
+
+        vector<uint16_t> samples = srswindow->samples;
+        for(uint isample=0; isample<samples.size(); isample++) {
+          pedestalAPV[srswindow->apv_id][isample] += samples[isample]/128.;
+        }
+    }
+
+    // loop over GEM channels and subtract pedestal
     for (const auto& window : windowrawdata) {
 	const DGEMSRSWindowRawData* srswindow;
 	window->GetSingle(srswindow);
 
 	int plane = window->plane;
 	int strip = window->strip;
-
-	// loop over samples for each channel
-	vector<uint16_t> samples = srswindow->samples;
-	double pedestal = 0;
-	const int nsample_pedestal = 4;
-	for(int i=0; i<nsample_pedestal; i++) {
-		pedestal += samples[i];
-	}
-	pedestal /= nsample_pedestal;
-
+        vector<uint16_t> samples = srswindow->samples;
+
 	// loop over samples to get pulse peak and time
 	double pulse_peak = 0;
 	double pulse_time = 0;
-	for(uint isample=nsample_pedestal; isample<samples.size(); isample++) {
-		int adc_zs = -1 * (samples[isample]-pedestal); // invert zero suppressed ADC
+	for(uint isample=1; isample<samples.size(); isample++) {
+                int adc_zs = -1. * (samples[isample]-pedestalAPV[srswindow->apv_id][isample]); // invert zero suppressed ADC
+                //cout<<samples[isample]<<" "<<pedestalAPV[srswindow->apv_id][isample]<<" "<<adc_zs<<endl;
 		if(adc_zs > pulse_peak) {
 			pulse_peak = adc_zs;
 			pulse_time = isample * ns_per_sample;
 		}
 	}
 
 	if(pulse_peak < pulse_peak_threshold) continue;
-	    
+
 	// Build hit object
 	DGEMHit *hit = new DGEMHit;
 	hit->pulse_height = pulse_peak;
@@ -121,7 +136,6 @@ jerror_t DGEMHit_factory::evnt(JEventLoop *loop, uint64_t eventnumber)
 	// Apply calibration constants
 	double T = pulse_time;
 	hit->t = T;
-	//hit->t = hit->t + t_base[plane] - time_offsets[plane][strip];
 
 	hit->AddAssociatedObject(window);
 	_data.push_back(hit);

src/libraries/TRD/DTRDHit_factory.cc

@@ -92,7 +92,7 @@ jerror_t DTRDHit_factory::evnt(JEventLoop *loop, uint64_t eventnumber)
 
 	    // separate theresholds for Wire and GEM TRD
 	    if((digihit->plane == 4 || digihit->plane == 5 || digihit->plane == 0 || digihit->plane == 1) && pulse_height < 200) continue;
-	    if((digihit->plane == 2 || digihit->plane == 6) && pulse_height < 300) continue;
+	    if((digihit->plane == 2 || digihit->plane == 6) && pulse_height < 350) continue;
 
 	    // Time cut now
 	    double T = (double)digihit->pulse_time * 0.8;

src/plugins/Analysis/TRD_hists/JEventProcessor_TRD_hists.cc

@@ -36,11 +36,15 @@ const int NAPVchannels = 128;
 const int NGEMsamples = 21;
 
 static TH1I *trd_num_events;
+static TH1I *num_tracks;
 
-static TH2I *hWireTRDPoint_TrackX, *hWireTRDPoint_TrackY, *hGEMTRDHit_TrackX;
-static TH2I *hWireTRDPoint_DeltaXY, *hGEMTRDHit_DeltaX_T;
+static TH2I *hWireTRDPoint_TrackX, *hWireTRDPoint_TrackY;
+static TH2I *hGEMTRDHit_TrackX, *hGEMTRDHit_TrackY, *hGEMTRDHit_DeltaXY;
+static TH2I *hWireTRDPoint_DeltaXY, *hWireTRDPoint_DeltaXY_Pos, *hWireTRDPoint_DeltaXY_Neg;
+static TH2I *hGEMTRDHit_DeltaX_T;
 static TH1I *hWireTRDPoint_Time;
-static TH2I *hGEMSRSPoint_TrackX[5], *hGEMSRSPoint_TrackY[5], *hGEMSRSPoint_DeltaXY[5];
+static TH2I *hGEMSRSPoint_TrackX[5], *hGEMSRSPoint_TrackY[5];
+static TH2I *hGEMSRSPoint_DeltaXY[5], *hGEMSRSPoint_DeltaXY_Pos[5], *hGEMSRSPoint_DeltaXY_Neg[5];
 
 static TH2I *hWire_GEMTRDXstrip, *hWire_GEMTRDX_DeltaT;
 
@@ -80,25 +84,33 @@ jerror_t JEventProcessor_TRD_hists::init(void) {
     trdDir->cd();
     // book hists
     trd_num_events = new TH1I("trd_hists_num_events","TRD number of events",1,0.5,1.5);
+    num_tracks = new TH1I("num_tracks", "; # of straight (0) and time based (1) tracks",2,-0.5,1.5);
 
-    // digihit-level hists
     trdDir->cd();
+
+    // Straight track histograms...
     gDirectory->mkdir("StraightTracks")->cd();
 
     hWireTRDPoint_TrackX = new TH2I("WireTRDPoint_TrackX","; Wire TRD X (cm); Extrapolated track X (cm)",200,-55,55,200,-55,55);
     hWireTRDPoint_TrackY = new TH2I("WireTRDPoint_TrackY","; Wire TRD Y (cm); Extrapolated track Y (cm)",200,-85,-65,200,-85,-65);
     hWireTRDPoint_DeltaXY = new TH2I("WireTRDPoint_DeltaXY","; #Delta X (cm); #Delta Y (cm)",100,-5,5,100,-5,5);
+    hWireTRDPoint_DeltaXY_Pos = new TH2I("WireTRDPoint_DeltaXY_Pos","; #Delta X (cm); #Delta Y (cm)",100,-5,5,100,-5,5);
+    hWireTRDPoint_DeltaXY_Neg = new TH2I("WireTRDPoint_DeltaXY_Neg","; #Delta X (cm); #Delta Y (cm)",100,-5,5,100,-5,5);
     hWireTRDPoint_Time = new TH1I("WireTRDPoint_Time","; hit time",1000,0,1000);
 
     // GEM TRD
     hGEMTRDHit_TrackX = new TH2I("GEMTRDHit_TrackX","; GEM TRD X; Extrapolated track X",200,-55,55,200,-55,55);
+    hGEMTRDHit_TrackY = new TH2I("GEMTRDHit_TrackY","; GEM TRD Y; Extrapolated track Y",200,-85,-65,200,-85,-65);
+    hGEMTRDHit_DeltaXY = new TH2I("GEMTRDHit_DeltaXY","; #Delta X (cm); #Delta Y (cm)",100,-5,5,100,-5,5);
     hGEMTRDHit_DeltaX_T = new TH2I("GEMTRDHit_DeltaX_T","; #Delta X; hit time",1000,-20,20,1000,0,1000);
 
     // GEM SRS
     for(int i=0; i<5; i++) {
 	    hGEMSRSPoint_TrackX[i] = new TH2I(Form("GEMSRSPoint_TrackX_%d",i),Form("Package %d; GEM SRS X (cm); Extrapolated track X (cm)",i),200,-55,55,200,-55,55);
 	    hGEMSRSPoint_TrackY[i] = new TH2I(Form("GEMSRSPoint_TrackY_%d",i),Form("Package %d; GEM SRS Y (cm); Extrapolated track Y (cm)",i),200,-85,-65,200,-85,-65);
 	    hGEMSRSPoint_DeltaXY[i] = new TH2I(Form("GEMSRSPoint_DeltaXY_%d",i),Form("Package %d; #Delta X (cm); #Delta Y (cm)",i),500,-5,5,500,-5,5);
+	     hGEMSRSPoint_DeltaXY_Pos[i] = new TH2I(Form("GEMSRSPoint_DeltaXY_Pos_%d",i),Form("Package %d; #Delta X (cm); #Delta Y (cm)",i),500,-5,5,500,-5,5);
+	     hGEMSRSPoint_DeltaXY_Neg[i] = new TH2I(Form("GEMSRSPoint_DeltaXY_Neg_%d",i),Form("Package %d; #Delta X (cm); #Delta Y (cm)",i),500,-5,5,500,-5,5);
     }
 
     // GEM-Wire TRD correlatioin
@@ -126,6 +138,9 @@ jerror_t JEventProcessor_TRD_hists::brun(JEventLoop *eventLoop, int32_t runnumbe
     vector<double> z_trd;
     geom->GetTRDZ(z_trd);
 
+    const DMagneticFieldMap *bfield = dapp->GetBfield(runnumber);
+    dIsNoFieldFlag = ((dynamic_cast<const DMagneticFieldMapNoField*>(bfield)) != NULL);
+
     return NOERROR;
 }
 
@@ -181,7 +196,14 @@ jerror_t JEventProcessor_TRD_hists::evnt(JEventLoop *eventLoop, uint64_t eventnu
     eventLoop->Get(gem_points);
 
     vector<const DTrackWireBased*> straight_tracks;
-    eventLoop->Get(straight_tracks, "StraightLine");
+    vector<const DTrackTimeBased*> tracks;
+    if (dIsNoFieldFlag)
+	    eventLoop->Get(straight_tracks, "StraightLine");
+    else
+	    eventLoop->Get(tracks);
+
+    num_tracks->Fill(0.,(double)straight_tracks.size());
+    num_tracks->Fill(1.,(double)tracks.size());
 
     // FILL HISTOGRAMS
     // Since we are filling histograms local to this plugin, it will not interfere with other ROOT operations: can use plugin-wide ROOT fill lock
@@ -196,6 +218,10 @@ jerror_t JEventProcessor_TRD_hists::evnt(JEventLoop *eventLoop, uint64_t eventnu
     // check if have good extrapolated track with TRD wire point
     bool goodTrack = false;
 
+    /////////////////////////////
+    // Straight track analysis //
+    /////////////////////////////
+
     for (const auto& straight_track : straight_tracks) {
 	    vector<DTrackFitter::Extrapolation_t>extrapolations=straight_track->extrapolations.at(SYS_TRD);
 	    //cout<<"found straight track with "<<extrapolations.size()<<" extrapolations"<<endl;
@@ -221,6 +247,82 @@ jerror_t JEventProcessor_TRD_hists::evnt(JEventLoop *eventLoop, uint64_t eventnu
 			    }
 		    }
 
+		    // correlate GEM TRD with extrapolated tracks
+		    for (const auto& hit : hits) {
+			    if(hit->plane == 6 && fabs(extrapolation.position.Z() - 570.7) < 5.) {
+
+				    for (const auto& gem_hit : gem_hits) {
+					    if(gem_hit->plane == 7 && fabs(extrapolation.position.Z() - 570.7) < 5.) {		    
+						    // only look at tracks with good wire hit
+						    if(!goodTrack) continue;
+
+						    // choose particular region of GEMTRD XY plane
+						    //if(abs(hit->strip-17) > 3 || abs(gem_hit->strip-240) > 3)
+						    //	    continue;
+
+						    double locStripX = -5.25 + hit->strip*0.04;
+						    double locStripY = -80.0 + hit->strip*0.04;
+						    hGEMTRDHit_TrackX->Fill(locStripX, extrapolation.position.X());
+						    hGEMTRDHit_TrackY->Fill(locStripY, extrapolation.position.Y());
+
+						    double locDeltaX = locStripX - extrapolation.position.X();
+						    double locDeltaY = locStripY - extrapolation.position.Y();
+						    hGEMTRDHit_DeltaXY->Fill(locDeltaX, locDeltaY);
+
+						    hGEMTRDHit_DeltaX_T->Fill(locDeltaX, hit->t);
+					    }
+				    }
+			    }
+		    }
+
+		    // correlate GEM SRS with extrapolated tracks
+		    for (const auto& point : gem_points) {
+			    if(point->detector == 4 && fabs(extrapolation.position.Z() - 576.7) < 5.) { 
+
+				    hGEMSRSPoint_TrackX[point->detector]->Fill(point->x, extrapolation.position.X());
+				    hGEMSRSPoint_TrackY[point->detector]->Fill(point->y, extrapolation.position.Y());
+
+				    double locDeltaX = point->x - extrapolation.position.X();
+				    double locDeltaY = point->y - extrapolation.position.Y();
+				    hGEMSRSPoint_DeltaXY[point->detector]->Fill(locDeltaX, locDeltaY);
+			    }
+		    }
+	    }
+    }
+
+    /////////////////////////////
+    // Field on track analysis //
+    /////////////////////////////
+
+    for (const auto& track : tracks) {
+
+	    int charge = track->charge();
+	    vector<DTrackFitter::Extrapolation_t>extrapolations=track->extrapolations.at(SYS_TRD);
+	    //cout<<"found straight track with "<<extrapolations.size()<<" extrapolations"<<endl;
+
+	    for (const auto& extrapolation : extrapolations) {
+
+		    // correlate wire TRD with extrapolated tracks
+		    for (const auto& point : points) {
+
+			    if(point->detector == 0 && fabs(extrapolation.position.Z() - 548.8) < 5.) { 
+
+				    hWireTRDPoint_TrackX->Fill(point->x, extrapolation.position.X());
+				    hWireTRDPoint_TrackY->Fill(point->y, extrapolation.position.Y());
+
+				    double locDeltaX = point->x - extrapolation.position.X();
+				    double locDeltaY = point->y - extrapolation.position.Y();
+				    hWireTRDPoint_DeltaXY->Fill(locDeltaX, locDeltaY);
+				    if(charge > 0) hWireTRDPoint_DeltaXY_Pos->Fill(locDeltaX, locDeltaY);
+				    else hWireTRDPoint_DeltaXY_Neg->Fill(locDeltaX, locDeltaY);
+
+				    if(fabs(locDeltaX) < 5. && fabs(locDeltaY) < 5.) {
+					    hWireTRDPoint_Time->Fill(point->time);
+					    goodTrack = true;
+				    }
+			    }
+		    }
+
 		    // correlate GEM TRD with extrapolated tracks
 		    for (const auto& hit : hits) {
 			    if(hit->plane == 6 && fabs(extrapolation.position.Z() - 570.7) < 5.) {
@@ -246,11 +348,17 @@ jerror_t JEventProcessor_TRD_hists::evnt(JEventLoop *eventLoop, uint64_t eventnu
 				    double locDeltaX = point->x - extrapolation.position.X();
 				    double locDeltaY = point->y - extrapolation.position.Y();
 				    hGEMSRSPoint_DeltaXY[point->detector]->Fill(locDeltaX, locDeltaY);
+				    if(charge > 0) hGEMSRSPoint_DeltaXY_Pos[point->detector]->Fill(locDeltaX, locDeltaY);
+				    else hGEMSRSPoint_DeltaXY_Neg[point->detector]->Fill(locDeltaX, locDeltaY);
 			    }
 		    }
 	    }
     }
 
+    ///////////////////////////
+    // Plots for good tracks //
+    ///////////////////////////
+
     if(goodTrack) {
 	    for (const auto& hit : hits) {
 		    if(hit->plane != 0 && hit->plane != 4) continue; // only Wire TRD

src/plugins/Analysis/TRD_hists/JEventProcessor_TRD_hists.h

@@ -5,6 +5,8 @@
 #include <JANA/JEventProcessor.h>
 #include <TRACKING/DTrackFitter.h>
 #include <TRACKING/DTrackWireBased.h>
+#include <TRACKING/DTrackTimeBased.h>
+#include <HDGEOMETRY/DMagneticFieldMapNoField.h>
 
 class JEventProcessor_TRD_hists:public jana::JEventProcessor{
 public:
@@ -19,6 +21,7 @@ class JEventProcessor_TRD_hists:public jana::JEventProcessor{
     jerror_t erun(void); ///< Called everytime run number changes, provided brun has been called.
     jerror_t fini(void); ///< Called after last event of last event source has been processed.
     int wirePlaneOffset;
+    bool dIsNoFieldFlag;
 };
 
 #endif // _JEventProcessor_TRD_hists_