Cm hit id tracking

offline/packages/tpc/TpcClusterizer.cc

@@ -6,6 +6,7 @@
 #include <trackbase/TrkrClusterv3.h>
 #include <trackbase/TrkrClusterv4.h>
 #include <trackbase/TrkrClusterv5.h>
+#include <trackbase/TrkrClusterv6.h>
 #include <trackbase/TrkrClusterHitAssocv3.h>
 #include <trackbase/TrkrDefs.h>  // for hitkey, getLayer
 #include <trackbase/TrkrHit.h>
@@ -64,6 +65,7 @@ namespace
     unsigned short it = 0;
     unsigned short adc = 0;
     unsigned short edge = 0;
+    unsigned long long parentID = -1;
   };
 
   using vec_dVerbose = std::vector<std::vector<std::pair<int,int>>>;
@@ -94,7 +96,7 @@ namespace
     unsigned short maxHalfSizePhi = 0;
     double m_tdriftmax = 0;
     double sampa_tbias = 0;
-    int cluster_version = 4;
+    int cluster_version = 6;
     std::vector<assoc> association_vector;
     std::vector<TrkrCluster*> cluster_vector;
     int verbosity = 0;
@@ -258,7 +260,7 @@ namespace
     return;
   }
 
-  void get_cluster(int phibin, int tbin, const thread_data& my_data, const std::vector<std::vector<unsigned short>> &adcval, std::vector<ihit> &ihit_list, int &touch, int &edge)
+  void get_cluster(int phibin, int tbin, const thread_data& my_data, const std::vector<std::vector<unsigned short>> &adcval, unsigned long long parentID, std::vector<ihit> &ihit_list, int &touch, int &edge)
 	{
 	  // search along phi at the peak in t
 
@@ -277,6 +279,7 @@ namespace
 		hit.iphi = iphi;
 		hit.it = it;
 		hit.adc = adcval[iphi][it];
+		hit.parentID = parentID;
 		if(touch>0){
 		  if((iphi == (phibin - phidown))||
 		     (iphi == (phibin + phiup))){
@@ -317,6 +320,8 @@ namespace
       int clus_size = ihit_list.size();
       int max_adc  = 0;
       if(clus_size == 1) return;
+
+      std::vector<unsigned long long> parentIDs;
       //      std::cout << "process list" << std::endl;    
       std::vector<TrkrDefs::hitkey> hitkeyvec;
 
@@ -336,6 +341,8 @@ namespace
 	if(iphi < phibinlo) phibinlo = iphi;
 	if(it > tbinhi) tbinhi = it;
 	if(it < tbinlo) tbinlo = it;
+
+	parentIDs.push_back(iter->parentID);
 
 	// update phi sums
 	//	double phi_center = my_data.layergeom->get_phicenter(iphi);
@@ -369,8 +376,27 @@ namespace
       //      std::cout << "done process list" << std::endl;
       if (adc_sum < 10){
 	hitkeyvec.clear();
+	parentIDs.clear();
 	return;  // skip obvious noise "clusters"
       }  
+
+      unsigned long long commonParentID = -1;
+      int maxNpIDs = 0;
+      for(int pID=0; pID<(int)parentIDs.size(); pID++){
+
+	int nPIDs = 0;
+	for(int pID2=0; pID2<(int)parentIDs.size(); pID2++){
+	  if(pID == pID2) continue;
+	  if(parentIDs[pID] == parentIDs[pID2]) nPIDs++;
+	}
+	if(nPIDs > maxNpIDs){
+	  maxNpIDs = nPIDs;
+	  commonParentID = parentIDs[pID];
+	}
+      }
+      double pctParentID = 1.0*maxNpIDs/parentIDs.size();
+
+
       // This is the global position
       double clusiphi = iphi_sum / adc_sum;
       double clusphi = my_data.layergeom->get_phi(clusiphi);
@@ -493,6 +519,26 @@ namespace
 	my_data.cluster_vector.push_back(clus);
   b_made_cluster = true;
 	}
+      }else if(my_data.cluster_version==6){
+	//std::cout << "ver5" << std::endl;
+	//	std::cout << "clus num" << my_data.cluster_vector.size() << " X " << local(0) << " Y " << clust << std::endl;
+	if(sqrt(phi_err_square) > 0.01){
+	auto clus = new TrkrClusterv6;
+	//auto clus = std::make_unique<TrkrClusterv3>();
+	clus->setAdc(adc_sum);  
+	clus->setMaxAdc(max_adc); 
+	clus->setEdge(nedge);
+	clus->setPhiSize(phisize);
+	clus->setZSize(tsize);
+	clus->setSubSurfKey(subsurfkey);      
+	clus->setLocalX(local(0));
+	clus->setLocalY(clust);
+	clus->setPhiError(sqrt(phi_err_square));
+	clus->setZError(sqrt(t_err_square * pow(my_data.tGeometry->get_drift_velocity(),2)));
+	clus->setParentID(commonParentID);
+	clus->setPctParentID(pctParentID);
+	my_data.cluster_vector.push_back(clus);
+	}
       }
 
       if (my_data.fillClusHitsVerbose && b_made_cluster) {
@@ -557,7 +603,7 @@ namespace
       for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
 	   hitr != hitrangei.second;
 	   ++hitr){
-	
+
 	if( TpcDefs::getPad(hitr->first) - phioffset < 0 ){
 	  //std::cout << "WARNING phibin out of range: " << TpcDefs::getPad(hitr->first) - phioffset << " | " << phibins << std::endl;
 	  continue;
@@ -592,6 +638,7 @@ namespace
 	    thisHit.it = tbin;
 	    thisHit.adc = adc;
 	    thisHit.edge = 0;
+	    thisHit.parentID = hitr->second->getParentID();
 	    all_hit_map.insert(std::make_pair(adc, thisHit));
 	  }
 	  if(adc>my_data->threshold){
@@ -675,13 +722,14 @@ namespace
       ihit hiHit = iter->second;
       int iphi = hiHit.iphi;
       int it = hiHit.it;
+      unsigned long long iparentID = hiHit.parentID;
       //put all hits in the all_hit_map (sorted by adc)
       //start with highest adc hit
       // -> cluster around it and get vector of hits
       std::vector<ihit> ihit_list;
       int ntouch = 0;
       int nedge  =0;
-      get_cluster(iphi, it, *my_data, adcval, ihit_list, ntouch, nedge );
+      get_cluster(iphi, it, *my_data, adcval, iparentID, ihit_list, ntouch, nedge );
 
       // -> calculate cluster parameters
       // -> add hits to truth association

offline/packages/tpc/TpcClusterizer.h

@@ -67,7 +67,7 @@ class TpcClusterizer : public SubsysReco
   double SectorFiducialCut = 0.5;
   unsigned short MaxClusterHalfSizePhi = 3;
   unsigned short MaxClusterHalfSizeT = 5;
-  int cluster_version = 4;
+  int cluster_version = 6;
   double m_tdriftmax = 0;
   double AdcClockPeriod = 53.0;   // ns 
 

offline/packages/tpccalib/PHTpcCentralMembraneClusterizer.cc

@@ -21,7 +21,8 @@
 #include <trackbase/TrkrDefs.h>
 #include <trackbase/TpcDefs.h>
 #include <trackbase/TrkrCluster.h>
-#include <trackbase/CMFlashClusterv3.h>
+#include <trackbase/TrkrClusterv6.h>
+#include <trackbase/CMFlashClusterv4.h>
 #include <trackbase/CMFlashClusterContainerv1.h>
 #include <trackbase/ActsGeometry.h>
 #include <trackbase/TrkrCluster.h>
@@ -59,8 +60,8 @@ int PHTpcCentralMembraneClusterizer::InitRun(PHCompositeNode *topNode)
     hxy = new TH2F("hxy","cluster x:y",800,-100,100,800,-80,80);
     hz = new TH1F("hz","cluster z", 220, -2,2);
 
-    hz_pos = new TH1F("hz_pos","cluster z>0", 400, 0,110);
-    hz_neg = new TH1F("hz_neg","cluster z<0", 400, -110,0);
+    hz_pos = new TH1F("hz_pos","cluster z>0", 750, 0,200);
+    hz_neg = new TH1F("hz_neg","cluster z<0", 750, -200,0);
 
     hClustE[0]= new TH1F("hRawClusterEnergy","Cluster Energy Before Merging;E[?]",200,0,2000);
     hClustE[1] = new TH1F("hMatchedClusterEnergy","Pair Cluster Energy After Merging;E[?]",200,0,2000);
@@ -110,11 +111,15 @@ int PHTpcCentralMembraneClusterizer::process_event(PHCompositeNode *topNode)
   std::vector<int>i_pair; //vector for pair matching
   std::vector<float>energy;//vector for energy values of clusters
   std::vector<bool>isAcrossGap;
+  std::vector<unsigned long long>parentIDs;
+  std::vector<double>pctParentIDs;
   int nTpcClust = 0;
 
   double mean_z_content_plus = 0.0;
   double mean_z_content_minus = 0.0;
 
+  int clusCount = 0;
+
   //first loop over clusters to make mod phi histograms of each layer and each pair of layers
   for(const auto& hitsetkey:_cluster_map->getHitSetKeys(TrkrDefs::TrkrId::tpcId))
   {
@@ -124,6 +129,12 @@ int PHTpcCentralMembraneClusterizer::process_event(PHCompositeNode *topNode)
     {
 
       const auto& [cluskey, cluster] = *clusiter;
+
+      //      TrkrClusterv6 *clusterv6 = new TrkrClusterv6();
+      //      clusterv6->CopyFrom(clusiter->second);
+
+      //std::cout << "cluster " << clusCount << "   parentID: " << cluster->getParentID() << "   pctParentID: " << cluster->getPctParentID() << std::endl;
+
       auto glob = tgeometry->getGlobalPosition(cluskey, cluster);
       TVector3 tmp_pos(glob(0),glob(1),glob(2));
 
@@ -162,9 +173,9 @@ int PHTpcCentralMembraneClusterizer::process_event(PHCompositeNode *topNode)
 	if(layer < 54) hphi_reco_pair_neg[layer-7]->Fill(phiMod);
 	if(layer > 7) hphi_reco_pair_neg[layer-8]->Fill(phiMod);
       }
-
-    }
-  }
+      clusCount++;
+    }//end loop over clusters
+  }//end loop over hitsets
 
   for(int i=1; i<hz_pos->GetNbinsX(); i++){
     mean_z_content_plus += hz_pos->GetBinContent(i);
@@ -177,6 +188,8 @@ int PHTpcCentralMembraneClusterizer::process_event(PHCompositeNode *topNode)
   mean_z_content_minus = mean_z_content_minus/hz_neg->GetNbinsX();
 
 
+  std::cout << "mean_z_content_plus " << mean_z_content_plus << "   mean_z_content_minus " << mean_z_content_minus << std::endl;
+
   //use peak in z distributions to identify if there is a CM Flash. Peak should be >20% of entries (needs to be tuned)
   if(hz_pos->GetMaximum() < 5*mean_z_content_plus || hz_neg->GetMaximum() < 5*mean_z_content_minus) return Fun4AllReturnCodes::EVENT_OK;
 
@@ -240,21 +253,37 @@ int PHTpcCentralMembraneClusterizer::process_event(PHCompositeNode *topNode)
       //Find which phiMod bin cluster is in (with given layer) and if counts below threshold don't use
       if(z >= 0){
 	phiBin = hphi_reco_pos[lay-7]->GetXaxis()->FindBin(phiMod);
-	if(hphi_reco_pos[lay-7]->GetBinContent(phiBin) < m_moduloThreshold) aboveThreshold = false;
+	if(hphi_reco_pos[lay-7]->GetBinContent(phiBin) < (m_moduloThreshold + (mean_z_content_plus/18.0))) aboveThreshold = false;
       }else{
 	phiBin = hphi_reco_neg[lay-7]->GetXaxis()->FindBin(phiMod);
-	if(hphi_reco_neg[lay-7]->GetBinContent(phiBin) < m_moduloThreshold) aboveThreshold = false;
+	if(hphi_reco_neg[lay-7]->GetBinContent(phiBin) < (m_moduloThreshold + (mean_z_content_minus/18.0))) aboveThreshold = false;
       }
 
-      if( !aboveThreshold ) continue;
-
-      if(cluster->getAdc() < _min_adc_value) continue;
-      if( std::abs(z) < _min_z_value) continue;
+      if( !aboveThreshold ){
+	std::cout << "below threshold" << std::endl;
+	continue;
+      }      
+
+      if(cluster->getAdc() < _min_adc_value){
+	std::cout << "adc too low" << std::endl;
+	continue;
+      }
+      if( std::abs(z) < _min_z_value){
+	std::cout << "z too low" << std::endl;
+	continue;
+      }
 
       //require that z be within 1cm of peak in z distributions (separate for each side)
-      if( (z>=0 && std::abs(z - hz_pos->GetBinCenter(hz_pos->GetMaximumBin())) > 1.0) || (z<0 && std::abs(z - hz_neg->GetBinCenter(hz_neg->GetMaximumBin())) > 1.0) ) continue;
+      if( (z>=0 && std::abs(z - hz_pos->GetBinCenter(hz_pos->GetMaximumBin())) > 1.0) || (z<0 && std::abs(z - hz_neg->GetBinCenter(hz_neg->GetMaximumBin())) > 1.0) ){
+	std::cout << "z not close enough to peak" << std::endl;
+	continue;
+      }
+      //std::cout << "accepted cluster " << m_accepted_clusters << "   parentID: " << cluster->getParentID() << "   pctParentID: " << cluster->getPctParentID() << std::endl;
+
 
       ++m_accepted_clusters;
+
+
 
       i_pair.push_back(-1);
       energy.push_back(cluster->getAdc());
@@ -263,6 +292,8 @@ int PHTpcCentralMembraneClusterizer::process_event(PHCompositeNode *topNode)
       layer.push_back((int)(TrkrDefs::getLayer(cluskey)));
       side.push_back(TpcDefs::getSide(cluskey));
       isAcrossGap.push_back(false);
+      parentIDs.push_back(cluster->getParentID());
+      pctParentIDs.push_back(cluster->getPctParentID());
       if(Verbosity() > 0) std::cout << ":\t" << x << "\t" << y << "\t" << z <<std::endl;
     }
   }
@@ -306,9 +337,13 @@ int PHTpcCentralMembraneClusterizer::process_event(PHCompositeNode *topNode)
 
     int layerMatch = -1;
 
+    /*
     int nRowsMatch = 2;
     if(layer[i] >= 39) nRowsMatch = 4;
     else if(layer[i] >= 23) nRowsMatch = 3;
+    */
+
+    int nRowsMatch = 1;
 
     if( pos[i].Z() >= 0 ){
       if(layer[i] == 7){
@@ -417,6 +452,8 @@ int PHTpcCentralMembraneClusterizer::process_event(PHCompositeNode *topNode)
   std::vector<unsigned int> nclusters;
   std::vector<bool> isREdge;
   std::vector<std::pair<int,int> > pairNums;
+  std::vector<unsigned long long> parentID;
+  std::vector<double> pctParentID;
 
   //  int nR2 = 0;
   //int nR3 = 0;
@@ -496,7 +533,15 @@ int PHTpcCentralMembraneClusterizer::process_event(PHCompositeNode *topNode)
         nclusters.push_back(2);
 	if(isAcrossGap[i] && isAcrossGap[i_pair[i]]) isREdge.push_back(true);
 	else isREdge.push_back(false);
-
+
+	if(energy[i] > energy[i_pair[i]]){
+	  parentID.push_back(parentIDs[i]);
+	  pctParentID.push_back(energy[i]*pctParentIDs[i]/(energy[i]+energy[i_pair[i]]));
+	}else{
+	  parentID.push_back(parentIDs[i_pair[i]]);
+	  pctParentID.push_back(energy[i_pair[i]]*pctParentIDs[i_pair[i]]/(energy[i]+energy[i_pair[i]]));
+	}
+
 	tmp_pair.first = i;
 	tmp_pair.second = i_pair[i];
 	pairNums.push_back(tmp_pair);
@@ -512,14 +557,18 @@ int PHTpcCentralMembraneClusterizer::process_event(PHCompositeNode *topNode)
       }
     } else {
 
+      if(isAcrossGap[i]) continue;
+
       if(_histos)  hClustE[2]->Fill(energy[i]);
       // These single cluster cases have good phi, but do not have a good radius centroid estimate - may want to skip them, record nclusters and identify if across gap
       //      if(layer[i] == 7) isAcrossGap[i] = true;
-
+      
       aveenergy.push_back(energy[i]);
       avepos.push_back(pos[i]);
       nclusters.push_back(1);
       isREdge.push_back(isAcrossGap[i]);
+      parentID.push_back(parentIDs[i]);
+      pctParentID.push_back(pctParentIDs[i]);
       tmp_pair.first = i;
       tmp_pair.second = -1;
       pairNums.push_back(tmp_pair);
@@ -533,15 +582,16 @@ int PHTpcCentralMembraneClusterizer::process_event(PHCompositeNode *topNode)
 
   for(unsigned int iv = 0; iv <avepos.size(); ++iv)
   {
-    auto cmfc = new CMFlashClusterv3();
-
+    auto cmfc = new CMFlashClusterv4();
 
     cmfc->setX(avepos[iv].X());
     cmfc->setY(avepos[iv].Y());
     cmfc->setZ(avepos[iv].Z());
     cmfc->setAdc(aveenergy[iv]);
     cmfc->setNclusters(nclusters[iv]);
     cmfc->setIsRGap(isREdge[iv]);
+    cmfc->setParentID(parentID[iv]);
+    cmfc->setPctParentID(pctParentID[iv]);
 
     int pair1 = pairNums[iv].first;
     int pair2 = pairNums[iv].second;
@@ -626,6 +676,9 @@ int PHTpcCentralMembraneClusterizer::End(PHCompositeNode * /*topNode*/ )
     hDist2Adj->Write();
 
     for(int i=0; i<47; i++){
+      hphi_reco_pos[i]->Write();
+      hphi_reco_neg[i]->Write();
+
       hphi_reco_pair_pos[i]->Write();
       hphi_reco_pair_neg[i]->Write();
     }