ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHTpcCentralMembraneClusterizer.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHTpcCentralMembraneClusterizer.cc
2 
4 
6 #include <phool/getClass.h>
7 #include <phool/phool.h>
8 #include <TNtuple.h>
9 #include <TFile.h>
10 #include <TH1F.h>
11 #include <TH2F.h>
12 #include <string>
13 #include <TVector3.h>
14 #include <TCanvas.h>
15 #include <TGraph.h>
16 #include <cmath>
17 
19 #include <trackbase/TrkrDefs.h>
20 #include <tpc/TpcDefs.h>
32 #include <trackbase/TrkrHitv2.h>
33 #include <trackbase/TrkrHitSetv1.h>
34 
38 
39 
40 using namespace std;
41 
42 //____________________________________________________________________________..
44  SubsysReco(name)
45 {
46 
47 
48  // Make some histograms on an output file for diagnostics
49  char temp[500];
50  sprintf(temp,
51  "./eval_output/Energy_Histograms_%i.root", _process);
52  fout = new TFile(temp,"RECREATE");
53 
54  henergy = new TH1F("henergy", "cluster energy", 200, 0, 2000);
55  hxy = new TH2F("hxy","cluster x:y",800,-100,100,800,-80,80);
56  hz = new TH1F("hz","cluster z", 220, -2,2);
57 
58  hClustE[0]= new TH1F("hRawClusterEnergy","Cluster Energy Before Merging;E[?]",200,0,2000);
59  hClustE[1] = new TH1F("hMatchedClusterEnergy","Pair Cluster Energy After Merging;E[?]",200,0,2000);
60  hClustE[2] = new TH1F("hSoloClusterEnergy","Lone Cluster Energy After Merging;E[?]",200,0,2000);
61 
62  hDist=new TH1F("hDist","3D distance to nearby clusters on same padrow;dist[cm]",100,-1,10);
63  hDistRow=new TH2F("hDistRow","phi distance to nearby clusters vs (lower)row;dist[rad];padrow",100,-0.001,0.01,60,-0.5,59.5);
64  hDist2=new TH1F("hDist2","phi distance to nearby clusters on same padrow;dist[rad]",100,-0.001,0.01);
65  hDistRowAdj=new TH2F("hDistRowAdj","phi distance to nearby clusters vs (lower)row;dist[rad];(lower) padrow",100,-0.001,0.01,60,-0.5,59.5);
66  hDist2Adj=new TH1F("hDist2Adj","phi distance to nearby clusters on adjacent padrow;dist[rad]",100,-0.001,0.01);
67 
68 }
69 
70 //____________________________________________________________________________..
72 {
73 
74 }
75 
76 //____________________________________________________________________________..
78 {
79  int ret = GetNodes(topNode);
80  return ret;
81 }
82 
83 //____________________________________________________________________________..
85 {
86 
87  //local coord conversion below
88  ActsTrackingGeometry *tgeometry = findNode::getClass<ActsTrackingGeometry>(topNode,"ActsTrackingGeometry");
89  std::cout << "got tgeom" << std::endl;
90  ActsSurfaceMaps *surfmaps = findNode::getClass<ActsSurfaceMaps>(topNode,"ActsSurfaceMaps");
91  std::cout << "got surfmaps" << std::endl;
92  if(!tgeometry or !surfmaps)
93  {
94  std::cout << PHWHERE << "No Acts geometry on node tree. Can't continue."
95  << std::endl;
96  }
97 
98  if(Verbosity() > 0) std::cout << std::endl << "original size of cluster map: " << _cluster_map->size() << std::endl;
100 
101  std::vector<TVector3>pos; //position vector in cartesian
102  std::vector<int>layer; //cluster layer number
103  std::vector<int>i_pair; //vector for pair matching
104  std::vector<float>energy;//vector for energy values of clusters
105  int nTpcClust = 0;
106 
107  for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
108  hitsetitr != hitsetrange.second;
109  ++hitsetitr)
110  {
111  auto clusRange = _cluster_map->getClusters(hitsetitr->first);
112  for (auto clusiter = clusRange.first;
113  clusiter != clusRange.second; ++clusiter)
114  {
115  TrkrDefs::cluskey cluskey = clusiter->first;
116  auto trkrid = TrkrDefs::getTrkrId(cluskey);
117  if(trkrid != TrkrDefs::tpcId) continue;
118 
119  TrkrCluster *cluster = clusiter->second;
120  ActsTransformations transformer;
121  auto glob = transformer.getGlobalPosition(cluster,surfmaps,tgeometry);
122 
123  float x = glob(0);
124  float y = glob(1);
125  float z = glob(2);
126 
127  if(Verbosity() > 0)
128  {
129  unsigned int lyr = TrkrDefs::getLayer(cluskey);
130  unsigned short side = TpcDefs::getSide(cluskey);
131  std::cout << " z " << z << " side " << side << " layer " << lyr << " Adc " << cluster->getAdc() << " x " << x << " y " << y << std::endl;
132  }
133 
134  if(cluster->getAdc() < _min_adc_value) continue;
135 
136  i_pair.push_back(-1);
137  energy.push_back(cluster->getAdc());
138  nTpcClust++;
139  pos.push_back(TVector3(x,y,z));
140  layer.push_back((int)(TrkrDefs::getLayer(cluskey)));
141  if(Verbosity() > 0)
142  std::cout << ":\t" << x << "\t" << y << "\t" << z <<std::endl;
143  }
144  }
145 
146  TVector3 delta;
147  float dphi;
148  TVector2 delta2;
149 
150  for (int i=0; i<nTpcClust; ++i)
151  {
152  for (int j=i+1; j<nTpcClust; j++)
153  { //redundant to the 'adjacent row' check: if (i==j) continue; //can't match yourself.
154  if (abs(layer[i]-layer[j])==0)
155  {
156  delta=pos[i]-pos[j];
157  dphi=abs(pos[i].DeltaPhi(pos[j]));
158  if(_histos)
159  {
160  hDist->Fill(delta.Mag());
161  hDist2->Fill(dphi);
162  hDistRow->Fill(dphi,layer[i]);
163  }
164  }
165  if (abs(layer[i]-layer[j])==1)
166  { //match those centers to the known/expected stripe positions
167 
168  delta=pos[i]-pos[j];
169  dphi=abs(pos[i].DeltaPhi(pos[j]));
170  if(_histos)
171  {
172  hDist2Adj->Fill(dphi);
173  hDistRowAdj->Fill(dphi,layer[i]);
174  }
175  }
176  }
177  }
178 
179  //now for each cluster, find its nearest partner on an adjacent row:
180  const float maxphidist=0.003;//as read off the plots.
181  for (int i=0;i<nTpcClust;i++)
182  {
183  float bestphidist=maxphidist;
184  for (int j=0;j<nTpcClust;j++)
185  {
186  //redundant to the 'adjacent row' check: if (i==j) continue; //can't match yourself.
187  if (abs(layer[i]-layer[j])!=1)
188  {
189  continue; //must match to an ADJACENT row.
190  }
191  float newphidist=abs(pos[i].DeltaPhi(pos[j]));
192  if (newphidist<bestphidist)
193  {
194  i_pair[i]=j;
195  bestphidist=newphidist;
196  }
197  }
198  }
199 
200  // check to see if the cluster pairs each match each other
201  vector<bool>goodPair;
202  bool allGood=true;
203  int nGood=0;
204 
205  for (int i=0;i<nTpcClust;i++)
206  {
207  int myPair=i_pair[i];
208  int itsPair=myPair<0?-1:i_pair[myPair];
209  if (i!=itsPair){
210  goodPair.push_back(false);
211  allGood=false;
212  } else {
213  if (i<myPair) nGood++;
214  {
215  goodPair.push_back(true);
216  }
217  }
218  }
219  if (allGood)
220  {
221  printf("All Good!\n");
222  }
223  else
224  {
225  printf("nGood=%d out of %d\n",nGood,nTpcClust/2);
226  }
227 
228  //build the weighted cluster centers
229  //==========================
230  vector<float>aveenergy;
231  vector<TVector3> avepos;
232 
233  for (int i=0;i<nTpcClust;++i)
234  {
235  if(_histos) hClustE[0]->Fill(energy[i]);
236 
237  if (goodPair[i])
238  {
239  if (i_pair[i]>i)
240  {
241  if(_histos) hClustE[1]->Fill(energy[i]+energy[i_pair[i]]);
242 
243  aveenergy.push_back(energy[i]+energy[i_pair[i]]);
244  TVector3 temppos=energy[i]*pos[i];
245  temppos=temppos+(energy[i_pair[i]]*pos[i_pair[i]]);
246  temppos=temppos*(1./(energy[i]+energy[i_pair[i]]));
247  avepos.push_back(temppos);
248  }
249  }
250  else
251  {
252  if(_histos) hClustE[2]->Fill(energy[i]);
253 
254  aveenergy.push_back(energy[i]);
255  avepos.push_back(pos[i]);
256  }
257  }
258 
259  // Loop over the vectors and put the clusters on the node tree
260  //==============================================
261  if(Verbosity() > 1) std::cout << " vector size is " << avepos.size() << std::endl;
262 
263  for(unsigned int iv = 0; iv <avepos.size(); ++iv)
264  {
265  auto cmfc = new CMFlashClusterv1();
266 
267  cmfc->setX(avepos[iv].X());
268  cmfc->setY(avepos[iv].Y());
269  cmfc->setZ(avepos[iv].Z());
270  cmfc->setAdc(aveenergy[iv]);
271 
273 
274  }
275 
276  // read back the clusters and make some histograms
277 
278  auto clusrange = _corrected_CMcluster_map->getClusters();
279  for (auto cmitr = clusrange.first;
280  cmitr !=clusrange.second;
281  ++cmitr)
282  {
283  auto cmkey = cmitr->first;
284  auto cmclus = cmitr->second;
285 
286  if(Verbosity() > 0)
287  std::cout << "found cluster " << cmkey << " with adc " << cmclus->getAdc()
288  << " x " << cmclus->getX() << " y " << cmclus->getY() << " z " << cmclus->getZ()
289  << std::endl;
290 
291  if(_histos)
292  {
293  henergy->Fill(cmclus->getAdc());
294  hxy->Fill(cmclus->getX(), cmclus->getY());
295  hz->Fill(cmclus->getZ());
296  }
297  }
298 
300 }
301 
302 //____________________________________________________________________________..
304 {
305 
306  if(_histos)
307  {
308  fout->cd();
309 
310  henergy->Write();
311  hxy->Write();
312  hz->Write();
313 
314  hClustE[0]->Write();
315  hClustE[1]->Write();
316  hClustE[2]->Write();
317  hDist->Write();
318  hDistRow->Write();
319  hDist2->Write();
320  hDistRowAdj->Write();
321  hDist2Adj->Write();
322 
323  fout->Close();
324  }
325 
327 }
328 
329 //____________________________________________________________________________..
330 
332 {
333  //---------------------------------
334  // Get Objects off of the Node Tree
335  //---------------------------------
336 
337  _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
338  if (!_cluster_map)
339  {
340  cerr << PHWHERE << " ERROR: Can't find node TRKR_CLUSTER" << endl;
342  }
343 
344 
345  _cluster_hit_map = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
346  if (!_cluster_hit_map)
347  {
348  cerr << PHWHERE << " ERROR: Can't find node TRKR_CLUSTERHITASSOC" << endl;
350  }
351 
352 
353  _hitset_map = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
354  if (!_hitset_map)
355  {
356  cerr << PHWHERE << " ERROR: Can't find node TRKR_HITSET" << endl;
358  }
359 
360  _corrected_CMcluster_map = findNode::getClass<CMFlashClusterContainer>(topNode, "CORRECTED_CM_CLUSTER");
362  {
363  std::cout << "Creating node CORRECTED_CM_CLUSTER" << std::endl;
364 
365  PHNodeIterator iter(topNode);
366 
367  // Looking for the DST node
368  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
369  if (!dstNode)
370  {
371  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
373  }
374  PHNodeIterator dstiter(dstNode);
375  PHCompositeNode *DetNode =
376  dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
377  if (!DetNode)
378  {
379  DetNode = new PHCompositeNode("TRKR");
380  dstNode->addNode(DetNode);
381  }
382 
384  PHIODataNode<PHObject> *TrkrClusterContainerNode =
385  new PHIODataNode<PHObject>(_corrected_CMcluster_map, "CORRECTED_CM_CLUSTER", "PHObject");
386  DetNode->addNode(TrkrClusterContainerNode);
387  }
388 
390 }
391 
392