51 "./eval_output/Energy_Histograms_%i.root",
_process);
52 fout =
new TFile(temp,
"RECREATE");
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);
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);
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);
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)
94 std::cout <<
PHWHERE <<
"No Acts geometry on node tree. Can't continue."
101 std::vector<TVector3>
pos;
102 std::vector<int>
layer;
103 std::vector<int>i_pair;
108 hitsetitr != hitsetrange.second;
112 for (
auto clusiter = clusRange.first;
113 clusiter != clusRange.second; ++clusiter)
131 std::cout <<
" z " << z <<
" side " << side <<
" layer " << lyr <<
" Adc " << cluster->
getAdc() <<
" x " << x <<
" y " << y << std::endl;
136 i_pair.push_back(-1);
137 energy.push_back(cluster->
getAdc());
139 pos.push_back(TVector3(x,y,z));
142 std::cout <<
":\t" << x <<
"\t" << y <<
"\t" << z <<std::endl;
150 for (
int i=0; i<nTpcClust; ++i)
152 for (
int j=i+1; j<nTpcClust; j++)
154 if (
abs(layer[i]-layer[j])==0)
157 dphi=
abs(pos[i].DeltaPhi(pos[j]));
160 hDist->Fill(delta.Mag());
165 if (
abs(layer[i]-layer[j])==1)
169 dphi=
abs(pos[i].DeltaPhi(pos[j]));
180 const float maxphidist=0.003;
181 for (
int i=0;i<nTpcClust;i++)
183 float bestphidist=maxphidist;
184 for (
int j=0;j<nTpcClust;j++)
187 if (
abs(layer[i]-layer[j])!=1)
191 float newphidist=
abs(pos[i].DeltaPhi(pos[j]));
192 if (newphidist<bestphidist)
195 bestphidist=newphidist;
201 vector<bool>goodPair;
205 for (
int i=0;i<nTpcClust;i++)
207 int myPair=i_pair[i];
208 int itsPair=myPair<0?-1:i_pair[myPair];
210 goodPair.push_back(
false);
213 if (i<myPair) nGood++;
215 goodPair.push_back(
true);
225 printf(
"nGood=%d out of %d\n",nGood,nTpcClust/2);
230 vector<float>aveenergy;
231 vector<TVector3> avepos;
233 for (
int i=0;i<nTpcClust;++i)
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);
254 aveenergy.push_back(energy[i]);
255 avepos.push_back(pos[i]);
261 if(
Verbosity() > 1) std::cout <<
" vector size is " << avepos.size() << std::endl;
263 for(
unsigned int iv = 0; iv <avepos.size(); ++iv)
267 cmfc->setX(avepos[iv].
X());
268 cmfc->setY(avepos[iv].
Y());
269 cmfc->setZ(avepos[iv].
Z());
270 cmfc->setAdc(aveenergy[iv]);
279 for (
auto cmitr = clusrange.first;
280 cmitr !=clusrange.second;
283 auto cmkey = cmitr->first;
284 auto cmclus = cmitr->second;
287 std::cout <<
"found cluster " << cmkey <<
" with adc " << cmclus->getAdc()
288 <<
" x " << cmclus->getX() <<
" y " << cmclus->getY() <<
" z " << cmclus->getZ()
293 henergy->Fill(cmclus->getAdc());
294 hxy->Fill(cmclus->getX(), cmclus->getY());
295 hz->Fill(cmclus->getZ());
337 _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode,
"TRKR_CLUSTER");
340 cerr <<
PHWHERE <<
" ERROR: Can't find node TRKR_CLUSTER" << endl;
345 _cluster_hit_map = findNode::getClass<TrkrClusterHitAssoc>(topNode,
"TRKR_CLUSTERHITASSOC");
348 cerr <<
PHWHERE <<
" ERROR: Can't find node TRKR_CLUSTERHITASSOC" << endl;
353 _hitset_map = findNode::getClass<TrkrHitSetContainer>(topNode,
"TRKR_HITSET");
356 cerr <<
PHWHERE <<
" ERROR: Can't find node TRKR_HITSET" << endl;
363 std::cout <<
"Creating node CORRECTED_CM_CLUSTER" << std::endl;
371 std::cout <<
PHWHERE <<
"DST Node missing, doing nothing." << std::endl;
386 DetNode->
addNode(TrkrClusterContainerNode);