Merge pull request #860 from JeffersonLab/AddFDCPositionsToRESTv2

Add fdc positions to res tv2

src/libraries/ANALYSIS/DEventWriterROOT.cc

@@ -5,6 +5,7 @@ static bool BCAL_VERBOSE_OUTPUT = false;
 static bool FCAL_VERBOSE_OUTPUT = false;
 static bool CCAL_VERBOSE_OUTPUT = false;
 static bool DIRC_OUTPUT = true;
+static bool FDC_VERBOSE_OUTPUT = true;
 
 static bool STORE_PULL_INFO = false;
 static bool STORE_ERROR_MATRIX_INFO = false;
@@ -55,6 +56,14 @@ void DEventWriterROOT::Initialize(JEventLoop* locEventLoop)
 	dTargetCenterZ = 65.0;
 	locGeometry->GetTargetZ(dTargetCenterZ);
 
+	// Get upstream positions of each FDC package
+	vector<double>fdc_z_wires;
+	locGeometry->GetFDCZ(fdc_z_wires);
+	dFdcPackages[0]=fdc_z_wires[0]-1.; // just upstream
+	dFdcPackages[1]=fdc_z_wires[6]-1.;
+	dFdcPackages[2]=fdc_z_wires[12]-1.;
+	dFdcPackages[3]=fdc_z_wires[18]-1.;
+
 	//CREATE TTREES
 	for(auto& locVertexInfo : locVertexInfos)
 	{
@@ -271,6 +280,8 @@ TMap* DEventWriterROOT::Create_UserInfoMaps(DTreeBranchRegister& locBranchRegist
 		locMiscInfoMap->Add(new TObjString("REST:JANACALIBCONTEXT"), new TObjString(REST_JANA_CALIB_CONTEXT.c_str()));
 
 	// Note: adding these parameters (e.g. in hd_root) will create warnings with "<-- NO DEFAULT! (TYPO?)". Safe to ignore these.
+	if(gPARMS->Exists("ANALYSIS:FDC_VERBOSE_ROOT_OUTPUT"))
+		{gPARMS->GetParameter("ANALYSIS:FDC_VERBOSE_ROOT_OUTPUT", FDC_VERBOSE_OUTPUT); cout << "ANALYSIS:FDC_VERBOSE_ROOT_OUTPUT set to " << FDC_VERBOSE_OUTPUT << ", IGNORE the \"<-- NO DEFAULT! (TYPO?)\" message " << endl;}
 	if(gPARMS->Exists("ANALYSIS:BCAL_VERBOSE_ROOT_OUTPUT"))
 		{gPARMS->GetParameter("ANALYSIS:BCAL_VERBOSE_ROOT_OUTPUT", BCAL_VERBOSE_OUTPUT); cout << "ANALYSIS:BCAL_VERBOSE_ROOT_OUTPUT set to " << BCAL_VERBOSE_OUTPUT << ", IGNORE the \"<-- NO DEFAULT! (TYPO?)\" message " << endl;}
 	if(gPARMS->Exists("ANALYSIS:FCAL_VERBOSE_ROOT_OUTPUT"))
@@ -697,6 +708,12 @@ void DEventWriterROOT::Create_Branches_ChargedHypotheses(DTreeBranchRegister& lo
 	locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "dEdx_CDC"), locArraySizeString, dInitNumTrackArraySize);
 	locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "dEdx_CDC_integral"), locArraySizeString, dInitNumTrackArraySize);
 	locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "dEdx_FDC"), locArraySizeString, dInitNumTrackArraySize);
+	if (FDC_VERBOSE_OUTPUT){
+	  for (unsigned int j=0;j<4;j++){
+	    locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, dFDCxLeaves[j].c_str()), locArraySizeString, dInitNumTrackArraySize); 
+	    locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, dFDCyLeaves[j].c_str()), locArraySizeString, dInitNumTrackArraySize);  
+	  }
+	}
 
 	//TIMING INFO
 	locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "HitTime"), locArraySizeString, dInitNumTrackArraySize);
@@ -1822,7 +1839,31 @@ void DEventWriterROOT::Fill_ChargedHypo(DTreeFillData* locTreeFillData, unsigned
 	locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "dEdx_CDC"), locChargedTrackHypothesis->Get_dEdx_CDC_amp(), locArrayIndex);
 	locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "dEdx_CDC_integral"), locChargedTrackHypothesis->Get_dEdx_CDC_int(), locArrayIndex);
 	locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "dEdx_FDC"), locTrackTimeBased->ddEdx_FDC, locArrayIndex);
-
+
+	if (FDC_VERBOSE_OUTPUT){
+	  if (locTrackTimeBased->extrapolations.find(SYS_FDC) != locTrackTimeBased->extrapolations.end()) {
+	    vector<DTrackFitter::Extrapolation_t>locExtraps=locTrackTimeBased->extrapolations.at(SYS_FDC);
+	    int locPackageIndex=-1,locOldIndex=-1;
+	    for (unsigned int j=0;j<locExtraps.size();j++){
+	      DVector3 locPos=locExtraps[j].position;
+	      double locZ=locPos.z();
+
+	      if (locZ>dFdcPackages[3]) locPackageIndex=3;
+	      else if (locZ>dFdcPackages[2]) locPackageIndex=2;
+	      else if (locZ>dFdcPackages[1]) locPackageIndex=1;
+	      else locPackageIndex=0;
+
+	      if (locPackageIndex!=locOldIndex){
+		locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, dFDCxLeaves[locPackageIndex].c_str()), locPos.x(), locArrayIndex);
+		locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, dFDCyLeaves[locPackageIndex].c_str()), locPos.y(), locArrayIndex);
+	      }
+
+	      if (locPackageIndex==3) break;
+	      locOldIndex=locPackageIndex;
+	    }
+	  }
+	}
+
 	//HIT ENERGY
 	double locTOFdEdx = (locChargedTrackHypothesis->Get_TOFHitMatchParams() != NULL) ? locChargedTrackHypothesis->Get_TOFHitMatchParams()->dEdx : 0.0;
 	locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "dEdx_TOF"), locTOFdEdx, locArrayIndex);

src/libraries/ANALYSIS/DEventWriterROOT.h

@@ -90,6 +90,9 @@ class DEventWriterROOT : public JObject
 		unsigned int dInitNumComboArraySize;
 
 		double dTargetCenterZ;
+  double dFdcPackages[4];
+  vector<string>dFDCxLeaves={"FDC1_X","FDC2_X","FDC3_X","FDC4_X"};
+  vector<string>dFDCyLeaves={"FDC1_Y","FDC2_Y","FDC3_Y","FDC4_Y"};
 
 		//DEFAULT ACTIONS LISTED SEPARATELY FROM CUSTOM (in case in derived class user does something bizarre)
 		map<const DReaction*, DCutAction_ThrownTopology*> dCutActionMap_ThrownTopology;

src/libraries/HDDM/DEventSourceREST.cc

@@ -1408,6 +1408,22 @@ jerror_t DEventSourceREST::Extract_DTrackTimeBased(hddm_r::HDDM *record,
       tra->setErrorMatrix(loc7x7ErrorMatrix);
       (*loc7x7ErrorMatrix)(6, 6) = fit.getT0err()*fit.getT0err();
 
+      // Positions at each FDC package
+      const hddm_r::FdcTrackPosList locFdcTrackPosList = iter->getFdcTrackPoses();
+      hddm_r::FdcTrackPosList::iterator locFdcTrackPosIterator = locFdcTrackPosList.begin();
+      if (locFdcTrackPosIterator!=locFdcTrackPosList.end()){
+	// Create the extrapolation vector
+	vector<DTrackFitter::Extrapolation_t>myvector;
+	tra->extrapolations.emplace(SYS_FDC,myvector);
+	for(; locFdcTrackPosIterator != locFdcTrackPosList.end(); ++locFdcTrackPosIterator){
+	  DVector3 pos(locFdcTrackPosIterator->getX(),
+		       locFdcTrackPosIterator->getY(),
+		       locFdcTrackPosIterator->getZ());
+	  DVector3 mom;
+	  tra->extrapolations[SYS_FDC].push_back(DTrackFitter::Extrapolation_t(pos,mom,0.,0.));
+	}
+      }
+
       // Track parameters at exit of tracking volume
       const hddm_r::ExitParamsList& locExitParamsList = iter->getExitParamses();
       hddm_r::ExitParamsList::iterator locExitParamsIterator = locExitParamsList.begin();	

src/libraries/HDDM/DEventWriterREST.cc

@@ -28,6 +28,9 @@ DEventWriterREST::DEventWriterREST(JEventLoop* locEventLoop, string locOutputFil
 	}
 	japp->Unlock("RESTWriter");
 
+	REST_WRITE_FDC_TRACK_POS=true;
+	gPARMS->SetDefaultParameter("REST:WRITE_FDC_TRACK_POS", REST_WRITE_FDC_TRACK_POS,"Add track positions at each FDC package");
+
 	REST_WRITE_TRACK_EXIT_PARAMS=true;
 	gPARMS->SetDefaultParameter("REST:WRITE_TRACK_EXIT_PARAMS", REST_WRITE_TRACK_EXIT_PARAMS,"Add track parameters at exit to tracking volume");
 
@@ -549,6 +552,24 @@ bool DEventWriterREST::Write_RESTEvent(JEventLoop* locEventLoop, string locOutpu
                elo3().setDEdxAmp(tracks[i]->ddEdx_FDC_amp_trunc[it]);
             }
 		}
+		if (REST_WRITE_FDC_TRACK_POS){
+		  if (tracks[i]->extrapolations.find(SYS_FDC) != tracks[i]->extrapolations.end()) {
+		    vector<DTrackFitter::Extrapolation_t>extraps=tracks[i]->extrapolations.at(SYS_FDC);
+		    double oldz=0.;
+		    for (unsigned int k=0;k<extraps.size();k++){
+		      DVector3 pos=extraps[k].position;
+		      // Write out one position per package
+		      if (pos.z()>oldz+20.){
+			hddm_r::FdcTrackPosList locFdcTrackPoses = tra().addFdcTrackPoses(1);
+			locFdcTrackPoses().setX(pos.x());
+			locFdcTrackPoses().setY(pos.y());
+			locFdcTrackPoses().setZ(pos.z());
+			oldz=pos.z();
+		      }
+		    }
+		  }
+		}
+
 		if (REST_WRITE_TRACK_EXIT_PARAMS){
 		  if (tracks[i]->extrapolations.find(SYS_NULL) != tracks[i]->extrapolations.end()) {
 		    vector<DTrackFitter::Extrapolation_t>extraps=tracks[i]->extrapolations.at(SYS_NULL);

src/libraries/HDDM/DEventWriterREST.h

@@ -62,6 +62,7 @@ class DEventWriterREST : public JObject
 		bool REST_WRITE_FMWPC_HITS;
 		bool REST_WRITE_CCAL_SHOWERS;
 		bool REST_WRITE_TRACK_EXIT_PARAMS;
+		bool REST_WRITE_FDC_TRACK_POS;
 		bool ADD_FCAL_DATA_FOR_CPP;
 		bool REST_WRITE_FCAL_HITS;
 

src/libraries/HDDM/rest.xml

@@ -82,6 +82,7 @@
         <CDCdEdxTrunc minOccurs="0" maxOccurs="unbounded" ntrunc="int" dx="float" dEdx="float" dxAmp="float" dEdxAmp="float" lunit="cm" dEdx_unit="GeV/cm"/>
         <FDCdEdxTrunc minOccurs="0" maxOccurs="unbounded" ntrunc="int" dx="float" dEdx="float" dxAmp="float" dEdxAmp="float" lunit="cm" dEdx_unit="GeV/cm"/>
       </dEdxDC>
+      <fdcTrackPos maxOccurs="4" minOccurs="0" x="float" y="float" z="float" lunit="cm"/>
       <exitParams maxOccurs="1" minOccurs="0" 
 		  x1="float" y1="float" z1="float"
                   px1="float" py1="float" pz1="float"