Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHTpcCentralMembraneMatcher.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHTpcCentralMembraneMatcher.cc
1 
7 #include "PHTpcCentralMembraneMatcher.h"
8 
9 #include <fun4all/Fun4AllReturnCodes.h>
10 
11 #include <phool/PHCompositeNode.h>
12 #include <phool/getClass.h>
13 #include <phool/phool.h>
14 #include <trackbase/CMFlashClusterv3.h>
15 #include <trackbase/CMFlashClusterContainerv1.h>
16 #include <trackbase/CMFlashDifferencev1.h>
17 #include <trackbase/CMFlashDifferenceContainerv1.h>
18 
19 #include <TFile.h>
20 #include <TGraph.h>
21 #include <TH1F.h>
22 #include <TH1D.h>
23 #include <TH2F.h>
24 #include <TF1.h>
25 #include <TNtuple.h>
26 #include <TString.h>
27 #include <TVector3.h>
28 #include <TStyle.h>
29 
30 #include <cmath>
31 #include <set>
32 #include <string>
33 
34 namespace
35 {
36  template<class T> inline constexpr T delta_phi(const T& phi)
37  {
38  if (phi > M_PI) return phi - 2. * M_PI;
39  else if (phi <= -M_PI) return phi + 2.* M_PI;
40  else return phi;
41  }
42 
43  template<class T> inline constexpr T square( const T& x ) { return x*x; }
44 
45  template<class T> inline T get_r( const T& x, const T& y ) { return std::sqrt( square(x) + square(y) ); }
46 
47  // stream acts vector3
48  [[maybe_unused]] std::ostream& operator << (std::ostream& out, const Acts::Vector3& v )
49  {
50  out << "(" << v.x() << ", " << v.y() << ", " << v.z() << ")";
51  return out;
52  }
53 
55  [[maybe_unused]] void normalize_distortions( TpcDistortionCorrectionContainer* dcc )
56  {
57  // loop over side
58  for( unsigned int i = 0; i<2; ++i )
59  {
60 
61  // loop over bins in entries
62  for( int ip = 0; ip < dcc->m_hentries[i]->GetNbinsX(); ++ip )
63  for( int ir = 0; ir < dcc->m_hentries[i]->GetNbinsY(); ++ir )
64  {
65  // count number of times a given cell was filled
66  const auto entries = dcc->m_hentries[i]->GetBinContent( ip+1, ir+1 );
67  if( entries <= 1 ) continue;
68 
69  // normalize histograms
70  for( const auto& h:{dcc->m_hDRint[i], dcc->m_hDPint[i], dcc->m_hDZint[i]} )
71  {
72  h->SetBinContent( ip+1, ir+1, h->GetBinContent( ip+1, ir+1 )/entries );
73  h->SetBinError( ip+1, ir+1, h->GetBinError( ip+1, ir+1 )/entries );
74  }
75  }
76  }
77  }
78 
80  [[maybe_unused]] void fill_guarding_bins( TpcDistortionCorrectionContainer* dcc )
81  {
82 
83  // loop over side
84  for( unsigned int i = 0; i<2; ++i )
85  {
86  for( const auto& h:{dcc->m_hDRint[i], dcc->m_hDPint[i], dcc->m_hDZint[i], dcc->m_hentries[i]} )
87  {
88 
89  // fill guarding phi bins
90  /*
91  * we use 2pi periodicity to do that:
92  * - last valid bin is copied to first guarding bin;
93  * - first valid bin is copied to last guarding bin
94  */
95  const auto phibins = h->GetNbinsX();
96  const auto rbins = h->GetNbinsY();
97  for( int ir = 0; ir < rbins; ++ir )
98  {
99  // copy last valid bin to first guarding bin
100  h->SetBinContent( 1, ir+1, h->GetBinContent( phibins-1, ir+1 ) );
101  h->SetBinError( 1, ir+1, h->GetBinError( phibins-1, ir+1 ) );
102 
103  // copy first valid bin to last guarding bin
104  h->SetBinContent( phibins, ir+1, h->GetBinContent( 2, ir+1 ) );
105  h->SetBinError( phibins, ir+1, h->GetBinError( 2, ir+1 ) );
106  }
107 
108  // fill guarding r bins
109  for( int iphi = 0; iphi < phibins; ++iphi )
110  {
111  // copy first valid bin to first guarding bin
112  h->SetBinContent( iphi+1, 1, h->GetBinContent( iphi+1, 2 ) );
113  h->SetBinError( iphi+1, 1, h->GetBinError( iphi+1, 2 ) );
114 
115  // copy last valid bin to last guarding bin
116  h->SetBinContent( iphi+1, rbins, h->GetBinContent( iphi+1, rbins-1 ) );
117  h->SetBinError( iphi+1, rbins, h->GetBinError( iphi+1, rbins-1 ) );
118  }
119  }
120  }
121  }
122 
123 }
124 
125 //____________________________________________________________________________..
126 PHTpcCentralMembraneMatcher::PHTpcCentralMembraneMatcher(const std::string &name):
127  SubsysReco(name)
128 {
129  // calculate stripes center positions
130  CalculateCenters(nPads_R1, R1_e, nGoodStripes_R1_e, keepUntil_R1_e, nStripesIn_R1_e, nStripesBefore_R1_e, cx1_e, cy1_e);
131  CalculateCenters(nPads_R1, R1, nGoodStripes_R1, keepUntil_R1, nStripesIn_R1, nStripesBefore_R1, cx1, cy1);
132  CalculateCenters(nPads_R2, R2, nGoodStripes_R2, keepUntil_R2, nStripesIn_R2, nStripesBefore_R2, cx2, cy2);
133  CalculateCenters(nPads_R3, R3, nGoodStripes_R3, keepUntil_R3, nStripesIn_R3, nStripesBefore_R3, cx3, cy3);
134 }
135 
136 
137 //___________________________________________________________
138 void PHTpcCentralMembraneMatcher::set_grid_dimensions( int phibins, int rbins )
139 {
140  m_phibins = phibins;
141  m_rbins = rbins;
142 }
143 
144 /*
145 std::vector<double> PHTpcCentralMembraneMatcher::getRGaps( TH2F *r_phi ){
146 
147  TH1D *proj = r_phi->ProjectionY("R_proj",1,360);
148 
149  std::vector<double> pass1;
150 
151  for(int i=2; i<proj->GetNbinsX(); i++){
152  if(proj->GetBinContent(i) > 0.15*proj->GetMaximum() && proj->GetBinContent(i) >= proj->GetBinContent(i-1) && proj->GetBinContent(i) >= proj->GetBinContent(i+1)) pass1.push_back(proj->GetBinCenter(i));
153  }
154 
155  for(int i=0; i<(int)pass1.size()-1; i++){
156  if(pass1[i+1]-pass1[i] > 0.75) continue;
157 
158  if(proj->GetBinContent(proj->FindBin(pass1[i])) > proj->GetBinContent(proj->FindBin(pass1[i+1]))) pass1.erase(std::next(pass1.begin(), i+1));
159  else pass1.erase(std::next(pass1.begin(), i));
160 
161  i--;
162  }
163 
164  return pass1;
165 
166 }
167 */
168 
169 //get the average phi rotation using smoothed histograms
170 double PHTpcCentralMembraneMatcher::getPhiRotation_smoothed(TH1D *hitHist, TH1D *clustHist){
171 
172  //smooth the truth and cluster histograms
173  hitHist->Smooth();
174  clustHist->Smooth();
175 
176  //make a TF1 with a lambda function to make a function out of the truth histogram and shift it by a constant
177  TF1 *f1 = new TF1("f1",[&](double *x, double *p){ return p[0] * hitHist->GetBinContent(hitHist->FindBin((x[0] - p[1])>M_PI?x[0] - p[1] - 2*M_PI: x[0] - p[1])); }, -M_PI, M_PI, 2);
178  f1->SetParNames("A","shift");
179  f1->SetParameters(1.0,0.0);
180  // f1->SetParLimits(1,-M_PI/18,M_PI/18);
181 
182  f1->SetNpx(500);
183 
184  clustHist->Fit("f1","IQ");
185 
186  // clustHist->Draw();
187  //f1->Draw("same");
188 
189  return f1->GetParameter(1);
190 }
191 
192 //get vector with peak positions in Y (R) of histogram
193 std::vector<double> PHTpcCentralMembraneMatcher::getRPeaks(TH2F *r_phi){
194 
195  TH1D *proj = r_phi->ProjectionY("R_proj");
196 
197  std::vector<double> rPeaks;
198 
199  for(int i=2; i<proj->GetNbinsX(); i++){
200  //peak is when content is higher than 0.15* maximum value and content is greater than or equal to both adjacent bins
201  if(proj->GetBinContent(i) > 0.15*proj->GetMaximum() && proj->GetBinContent(i) >= proj->GetBinContent(i-1) && proj->GetBinContent(i) >= proj->GetBinContent(i+1)) rPeaks.push_back(proj->GetBinCenter(i));
202  }
203 
204  //if two peaks are within 0.75 cm of eachother, remove one with fewer counts
205  for(int i=0; i<(int)rPeaks.size()-1; i++){
206  if(rPeaks[i+1]-rPeaks[i] > 0.75) continue;
207  if(proj->GetBinContent(proj->FindBin(rPeaks[i])) > proj->GetBinContent(proj->FindBin(rPeaks[i+1]))) rPeaks.erase(rPeaks.begin()+i+1);
208  else rPeaks.erase(rPeaks.begin()+i);
209  i--;
210  }
211  return rPeaks;
212 }
213 
214 int PHTpcCentralMembraneMatcher::getClusterRMatch( std::vector<int> hitMatches, std::vector<double> clusterPeaks, double clusterR){
215 
216  double closestDist = 100.;
217  int closestPeak = -1;
218  //find cluster peak closest to position of passed cluster
219  for(int j=0; j<(int)clusterPeaks.size(); j++){
220  if(std::abs(clusterR - clusterPeaks[j]) < closestDist){
221  closestDist = std::abs(clusterR - clusterPeaks[j]);
222  closestPeak = j;
223  }
224  }
225 
226  //return hit match to cluster peak or -1 if closest peak failed (shouldn't be possible)
227  if(closestPeak != -1) return hitMatches[closestPeak];
228  else return -1;
229 
230 }
231 
232 //____________________________________________________________________________..
233 int PHTpcCentralMembraneMatcher::InitRun(PHCompositeNode *topNode)
234 {
235  if( m_savehistograms )
236  {
237 
238  static constexpr float max_dr = 5.0;
239  static constexpr float max_dphi = 0.05;
240 
241  fout.reset( new TFile(m_histogramfilename.c_str(),"RECREATE") );
242  hxy_reco = new TH2F("hxy_reco","reco cluster x:y",800,-100,100,800,-80,80);
243  hxy_truth = new TH2F("hxy_truth","truth cluster x:y",800,-100,100,800,-80,80);
244 
245  hdrdphi = new TH2F("hdrdphi","dr vs dphi",800,-max_dr,max_dr,800,-max_dphi,max_dphi);
246  hdrdphi->GetXaxis()->SetTitle("dr");
247  hdrdphi->GetYaxis()->SetTitle("dphi");
248 
249  hrdr = new TH2F("hrdr","dr vs r",800,0.0,80.0,800,-max_dr, max_dr);
250  hrdr->GetXaxis()->SetTitle("r");
251  hrdr->GetYaxis()->SetTitle("dr");
252 
253  hrdphi = new TH2F("hrdphi","dphi vs r",800,0.0,80.0,800,-max_dphi,max_dphi);
254  hrdphi->GetXaxis()->SetTitle("r");
255  hrdphi->GetYaxis()->SetTitle("dphi");
256 
257  hdphi = new TH1F("hdphi","dph",800,-max_dphi,max_dphi);
258  hdphi->GetXaxis()->SetTitle("dphi");
259 
260  hdr1_single = new TH1F("hdr1_single", "innner dr single", 200,-max_dr, max_dr);
261  hdr2_single = new TH1F("hdr2_single", "mid dr single", 200,-max_dr, max_dr);
262  hdr3_single = new TH1F("hdr3_single", "outer dr single", 200,-max_dr, max_dr);
263  hdr1_double = new TH1F("hdr1_double", "innner dr double", 200,-max_dr, max_dr);
264  hdr2_double = new TH1F("hdr2_double", "mid dr double", 200,-max_dr, max_dr);
265  hdr3_double = new TH1F("hdr3_double", "outer dr double", 200,-max_dr, max_dr);
266  hdrphi = new TH1F("hdrphi","r * dphi", 200, -0.05, 0.05);
267  hnclus = new TH1F("hnclus", " nclusters ", 3, 0., 3.);
268 
269  fout2.reset ( new TFile(m_histogramfilename2.c_str(),"RECREATE") );
270  match_ntup = new TNtuple("match_ntup","Match NTuple","event:truthR:truthPhi:recoR:recoPhi:recoZ:nclus:r1:phi1:e1:layer1:r2:phi2:e2:layer2");
271  }
272 
273  hit_r_phi = new TH2F("hit_r_phi","hit r vs #phi;#phi (rad); r (cm)",360,-M_PI,M_PI,500,0,100);
274 
275  clust_r_phi_pos = new TH2F("clust_r_phi_pos","clust R vs #phi Z>0;#phi (rad); r (cm)",360,-M_PI,M_PI,500,0,100);
276  clust_r_phi_neg = new TH2F("clust_r_phi_neg","clust R vs #phi Z<0;#phi (rad); r (cm)",360,-M_PI,M_PI,500,0,100);
277 
278 
279  // Get truth cluster positions
280  //=====================
281 
282  const double phi_petal = M_PI / 9.0; // angle span of one petal
283 
284 
285  /*
286  * utility function to
287  * - duplicate generated truth position to cover both sides of the central membrane
288  * - assign proper z,
289  * - insert in container
290  */
291  auto save_truth_position = [&](TVector3 source)
292  {
293  source.SetZ( +1 );
294  m_truth_pos.push_back( source );
295 
296  hit_r_phi->Fill(source.Phi(), source.Perp());
297 
298  source.SetZ( -1 );
299  m_truth_pos.push_back( source );
300 
301  hit_r_phi->Fill(source.Phi(), source.Perp());
302  };
303 
304  // inner region extended is the 8 layers inside 30 cm
305  for(int j = 0; j < nRadii; ++j)
306  for(int i=0; i < nGoodStripes_R1_e[j]; ++i)
307  for(int k =0; k<18; ++k)
308  {
309  TVector3 dummyPos(cx1_e[i][j], cy1_e[i][j], 0.0);
310  dummyPos.RotateZ(k * phi_petal);
311  save_truth_position(dummyPos);
312 
313  if(Verbosity() > 2)
314  std::cout << " i " << i << " j " << j << " k " << k << " x1 " << dummyPos.X() << " y1 " << dummyPos.Y()
315  << " theta " << std::atan2(dummyPos.Y(), dummyPos.X())
316  << " radius " << get_r( dummyPos.X(), dummyPos.y()) << std::endl;
317  if(m_savehistograms) hxy_truth->Fill(dummyPos.X(),dummyPos.Y());
318  }
319 
320  // inner region is the 8 layers outside 30 cm
321  for(int j = 0; j < nRadii; ++j)
322  for(int i=0; i < nGoodStripes_R1[j]; ++i)
323  for(int k =0; k<18; ++k)
324  {
325  TVector3 dummyPos(cx1[i][j], cy1[i][j], 0.0);
326  dummyPos.RotateZ(k * phi_petal);
327  save_truth_position(dummyPos);
328 
329  if(Verbosity() > 2)
330  std::cout << " i " << i << " j " << j << " k " << k << " x1 " << dummyPos.X() << " y1 " << dummyPos.Y()
331  << " theta " << std::atan2(dummyPos.Y(), dummyPos.X())
332  << " radius " << get_r( dummyPos.X(), dummyPos.y()) << std::endl;
333  if(m_savehistograms) hxy_truth->Fill(dummyPos.X(),dummyPos.Y());
334  }
335 
336  for(int j = 0; j < nRadii; ++j)
337  for(int i=0; i < nGoodStripes_R2[j]; ++i)
338  for(int k =0; k<18; ++k)
339  {
340  TVector3 dummyPos(cx2[i][j], cy2[i][j], 0.0);
341  dummyPos.RotateZ(k * phi_petal);
342  save_truth_position(dummyPos);
343 
344  if(Verbosity() > 2)
345  std::cout << " i " << i << " j " << j << " k " << k << " x1 " << dummyPos.X() << " y1 " << dummyPos.Y()
346  << " theta " << std::atan2(dummyPos.Y(), dummyPos.X())
347  << " radius " << get_r( dummyPos.X(), dummyPos.y()) << std::endl;
348  if(m_savehistograms) hxy_truth->Fill(dummyPos.X(),dummyPos.Y());
349  }
350 
351  for(int j = 0; j < nRadii; ++j)
352  for(int i=0; i < nGoodStripes_R3[j]; ++i)
353  for(int k =0; k<18; ++k)
354  {
355  TVector3 dummyPos(cx3[i][j], cy3[i][j], 0.0);
356  dummyPos.RotateZ(k * phi_petal);
357  save_truth_position(dummyPos);
358 
359  if(Verbosity() > 2)
360  std::cout << " i " << i << " j " << j << " k " << k << " x1 " << dummyPos.X() << " y1 " << dummyPos.Y()
361  << " theta " << std::atan2(dummyPos.Y(), dummyPos.X())
362  << " radius " << get_r( dummyPos.X(), dummyPos.y()) << std::endl;
363  if(m_savehistograms) hxy_truth->Fill(dummyPos.X(),dummyPos.Y());
364  }
365 
366  int ret = GetNodes(topNode);
367  return ret;
368 }
369 
370 //____________________________________________________________________________..
371 int PHTpcCentralMembraneMatcher::process_event(PHCompositeNode * /*topNode*/)
372 {
373  std::vector<TVector3> reco_pos;
374  std::vector<TVector3> pos1;
375  std::vector<TVector3> pos2;
376  std::vector<unsigned int> reco_nclusters;
377  std::vector<unsigned int> reco_adc;
378  std::vector<unsigned int> adc1;
379  std::vector<unsigned int> adc2;
380  std::vector<unsigned int> layer1;
381  std::vector<unsigned int> layer2;
382 
383  // reset output distortion correction container histograms
384  for( const auto& harray:{m_dcc_out->m_hDRint, m_dcc_out->m_hDPint, m_dcc_out->m_hDZint, m_dcc_out->m_hentries} )
385  {
386 
387  clust_r_phi_pos->Reset();
388  clust_r_phi_neg->Reset();
389 
390  if(!m_corrected_CMcluster_map || m_corrected_CMcluster_map->size() < 100){
391  m_event_index++;
392  return Fun4AllReturnCodes::EVENT_OK;
393  }
394 
395  for( const auto& h:harray ) { h->Reset(); } }
396 
397  // read the reconstructed CM clusters
398  auto clusrange = m_corrected_CMcluster_map->getClusters();
399  for (auto cmitr = clusrange.first;
400  cmitr !=clusrange.second;
401  ++cmitr)
402  {
403  const auto& [cmkey, cmclus_orig] = *cmitr;
404  CMFlashCluster *cmclus = cmclus_orig;
405  const unsigned int nclus = cmclus->getNclusters();
406  const unsigned int adc = cmclus->getAdc();
407 
408  if(m_useOnly_nClus2 && nclus != 2) continue;
409 
410  const bool isRGap = cmclus->getIsRGap();
411 
412 
413  // Do the static + average distortion corrections if the container was found
414  Acts::Vector3 pos(cmclus->getX(), cmclus->getY(), cmclus->getZ());
415  Acts::Vector3 apos1(cmclus->getX1(), cmclus->getY1(), cmclus->getZ1());
416  Acts::Vector3 apos2(cmclus->getX2(), cmclus->getY2(), cmclus->getZ2());
417  if( m_dcc_in_static){
418  pos = m_distortionCorrection.get_corrected_position( pos, m_dcc_in_static );
419  apos1 = m_distortionCorrection.get_corrected_position( apos1, m_dcc_in_static );
420  apos2 = m_distortionCorrection.get_corrected_position( apos2, m_dcc_in_static );
421  }
422  if( m_dcc_in_average){
423  pos = m_distortionCorrection.get_corrected_position( pos, m_dcc_in_average );
424  apos1 = m_distortionCorrection.get_corrected_position( apos1, m_dcc_in_average );
425  apos2 = m_distortionCorrection.get_corrected_position( apos2, m_dcc_in_average );
426  }
427 
428 
429 
430  TVector3 tmp_pos(pos[0], pos[1], pos[2]);
431  TVector3 tmp_pos1(apos1[0], apos1[1], apos1[2]);
432  TVector3 tmp_pos2(apos2[0], apos2[1], apos2[2]);
433 
434 
435  if(nclus == 1 && isRGap) continue;
436 
437  reco_pos.push_back(tmp_pos);
438  pos1.push_back(tmp_pos1);
439  pos2.push_back(tmp_pos2);
440  reco_nclusters.push_back(nclus);
441  reco_adc.push_back(adc);
442  adc1.push_back(cmclus->getAdc1());
443  adc2.push_back(cmclus->getAdc2());
444  layer1.push_back(cmclus->getLayer1());
445  layer2.push_back(cmclus->getLayer2());
446 
447  if(tmp_pos.Z() < 0) clust_r_phi_neg->Fill(tmp_pos.Phi(),tmp_pos.Perp());
448  else clust_r_phi_pos->Fill(tmp_pos.Phi(),tmp_pos.Perp());
449 
450 
451 
452  if(Verbosity())
453  {
454  double raw_rad = sqrt( cmclus->getX()*cmclus->getX() + cmclus->getY()*cmclus->getY() );
455  double corr_rad = sqrt( tmp_pos.X()*tmp_pos.X() + tmp_pos.Y()*tmp_pos.Y() );
456  std::cout << "found raw cluster " << cmkey << " with x " << cmclus->getX() << " y " << cmclus->getY() << " z " << cmclus->getZ() << " radius " << raw_rad << std::endl;
457  std::cout << " --- corrected positions: " << tmp_pos.X() << " " << tmp_pos.Y() << " " << tmp_pos.Z() << " radius " << corr_rad << std::endl;
458  }
459 
460  if(m_savehistograms)
461  {
462  hxy_reco->Fill(tmp_pos.X(), tmp_pos.Y());
463  }
464  }
465 
466 
467  //get global phi rotation for each module
468  m_clustRotation_pos[0] = getPhiRotation_smoothed(hit_r_phi->ProjectionX("hR1",151,206),clust_r_phi_pos->ProjectionX("cR1_pos",151,206));
469  m_clustRotation_pos[1] = getPhiRotation_smoothed(hit_r_phi->ProjectionX("hR2",206,290),clust_r_phi_pos->ProjectionX("cR2_pos",206,290));
470  m_clustRotation_pos[2] = getPhiRotation_smoothed(hit_r_phi->ProjectionX("hR3",290,499),clust_r_phi_pos->ProjectionX("cR3_pos",290,499));
471 
472  m_clustRotation_neg[0] = getPhiRotation_smoothed(hit_r_phi->ProjectionX("hR1",151,206),clust_r_phi_neg->ProjectionX("cR1_neg",151,206));
473  m_clustRotation_neg[1] = getPhiRotation_smoothed(hit_r_phi->ProjectionX("hR2",206,290),clust_r_phi_neg->ProjectionX("cR2_neg",206,290));
474  m_clustRotation_neg[2] = getPhiRotation_smoothed(hit_r_phi->ProjectionX("hR3",290,499),clust_r_phi_neg->ProjectionX("cR3_neg",290,499));
475 
476 
477  //get hit and cluster peaks
478  std::vector<double> hit_RPeaks = getRPeaks(hit_r_phi);
479  std::vector<double> clust_RPeaks_pos = getRPeaks(clust_r_phi_pos);
480  std::vector<double> clust_RPeaks_neg = getRPeaks(clust_r_phi_neg);
481 
482  //identify where gaps between module 1&2 and 2&3 are by finding largest distances between peaks
483  std::vector<double> clust_RGaps_pos;
484  int R12Gap_pos = -1;
485  int R23Gap_pos = -1;
486  for(int i=0; i<(int)clust_RPeaks_pos.size()-1; i++) clust_RGaps_pos.push_back(clust_RPeaks_pos[i+1] - clust_RPeaks_pos[i]);
487  int tmpGap_pos = std::distance(clust_RGaps_pos.begin(),std::max_element(clust_RGaps_pos.begin(),clust_RGaps_pos.end()));
488  if(tmpGap_pos > (int)clust_RGaps_pos.size()/2){
489  R23Gap_pos = tmpGap_pos;
490  R12Gap_pos = std::distance(clust_RGaps_pos.begin(),std::max_element(clust_RGaps_pos.begin(),clust_RGaps_pos.begin()+R23Gap_pos));
491  }else{
492  R12Gap_pos = tmpGap_pos;
493  R23Gap_pos = std::distance(clust_RGaps_pos.begin(),std::max_element(clust_RGaps_pos.begin()+R12Gap_pos+1,clust_RGaps_pos.end()));
494  }
495 
496  std::vector<double> clust_RGaps_neg;
497  int R12Gap_neg = -1;
498  int R23Gap_neg = -1;
499  for(int i=0; i<(int)clust_RPeaks_neg.size()-1; i++) clust_RGaps_neg.push_back(clust_RPeaks_neg[i+1] - clust_RPeaks_neg[i]);
500  int tmpGap_neg = std::distance(clust_RGaps_neg.begin(),std::max_element(clust_RGaps_neg.begin(),clust_RGaps_neg.end()));
501  if(tmpGap_neg > (int)clust_RGaps_neg.size()/2){
502  R23Gap_neg = tmpGap_neg;
503  R12Gap_neg = std::distance(clust_RGaps_neg.begin(),std::max_element(clust_RGaps_neg.begin(),clust_RGaps_neg.begin()+R23Gap_neg));
504  }else{
505  R12Gap_neg = tmpGap_neg;
506  R23Gap_neg = std::distance(clust_RGaps_neg.begin(),std::max_element(clust_RGaps_neg.begin()+R12Gap_neg+1,clust_RGaps_neg.end()));
507  }
508 
509  std::vector<int> hitMatches_pos;
510  //match cluster peaks to hit peaks using gap positions
511  for(int i=0; i<(int)clust_RPeaks_pos.size(); i++){
512 
513  //Module 1
514  if(i < (R12Gap_pos+1)){
515  //module 1-2 gap is between 15 & 16 in hitPeaks
516  hitMatches_pos.push_back(15 + i - R12Gap_pos );
517  //if multiple rows missing, offset for each one
518  if(clust_RGaps_pos[R12Gap_pos] > 3.6) hitMatches_pos[i] -= 1;
519  if(clust_RGaps_pos[R12Gap_pos] > 4.6) hitMatches_pos[i] -= 1;
520  if(clust_RGaps_pos[R12Gap_pos] > 5.8) hitMatches_pos[i] -= 1;
521  }
522  //module 1-2 gap is between 15 & 16
523  else if(i < (R23Gap_pos+1) && i >= (R12Gap_pos+1)) hitMatches_pos.push_back(15+1 + i - (R12Gap_pos+1));
524  //module 2-3 gap is between 22 & 23
525  else if(i >= R23Gap_pos+1) hitMatches_pos.push_back(23+1 + i - (R23Gap_pos+1));
526  }
527 
528  std::vector<int> hitMatches_neg;
529  for(int i=0; i<(int)clust_RPeaks_neg.size(); i++){
530 
531  if(i < (R12Gap_neg+1)){
532 
533  hitMatches_neg.push_back(15 + i - R12Gap_neg );
534  if(clust_RGaps_neg[R12Gap_neg] > 3.6) hitMatches_neg[i] -= 1;
535  if(clust_RGaps_neg[R12Gap_neg] > 4.6) hitMatches_neg[i] -= 1;
536  if(clust_RGaps_neg[R12Gap_neg] > 5.8) hitMatches_neg[i] -= 1;
537  }
538  else if(i < (R23Gap_neg+1) && i >= (R12Gap_neg+1)) hitMatches_neg.push_back(15+1 + i - (R12Gap_neg+1));
539  else if(i >= R23Gap_neg+1) hitMatches_neg.push_back(23+1 + i - (R23Gap_neg+1));
540  }
541 
542 
543  // Match reco and truth positions
544  //std::map<unsigned int, unsigned int> matched_pair;
545  std::vector<std::pair<unsigned int, unsigned int>> matched_pair;
546  std::vector<unsigned int> matched_nclus;
547 
548  std::vector<bool> hits_matched(m_truth_pos.size());
549  std::vector<bool> clusts_matched(reco_pos.size());
550 
551  //do iterative matching m_nMatchIter times
552  for(int matchIt=0; matchIt<m_nMatchIter; matchIt++){
553  // loop over truth positions
554  for(unsigned int i=0; i<m_truth_pos.size(); ++i)
555  {
556 
557  if(hits_matched[i]) continue;
558 
559  const double z1 = m_truth_pos[i].Z();
560  const double rad1= get_r( m_truth_pos[i].X(),m_truth_pos[i].Y());
561  const double phi1 = m_truth_pos[i].Phi();
562 
563  int hitRadIndex = -1;
564 
565  //get which hit radial index this it
566  for(int k=0; k<(int)hit_RPeaks.size(); k++){
567  if(std::abs(rad1 - hit_RPeaks[k]) < 0.5){
568  hitRadIndex = k;
569  break;
570  }
571  }
572 
573  //for iterative looping: identify closest phi
574  double prev_dphi = 1.1*m_phi_cut;
575  int matchJ = -1;
576 
577  // loop over cluster positions
578  for(unsigned int j = 0; j < reco_pos.size(); ++j)
579  {
580  if(clusts_matched[j]) continue;
581 
582  int angleR = -1;
583 
584  if(reco_pos[j].Perp() < 41) angleR = 0;
585  else if(reco_pos[j].Perp() >= 41 && reco_pos[j].Perp() < 58) angleR = 1;
586  if(reco_pos[j].Perp() >= 58) angleR = 2;
587 
588 
589  //const auto& nclus = reco_nclusters[j];
590  double phi2 = reco_pos[j].Phi();
591  const double z2 = reco_pos[j].Z();
592  const double rad2=get_r(reco_pos[j].X(), reco_pos[j].Y());
593  if(angleR != -1){
594  if(z2 > 0) phi2 -= m_clustRotation_pos[angleR];
595  else phi2 -= m_clustRotation_neg[angleR];
596  }
597 
598 
599  int clustRMatchIndex = -1;
600  if(z2 > 0) clustRMatchIndex = getClusterRMatch(hitMatches_pos, clust_RPeaks_pos, rad2);
601  else clustRMatchIndex = getClusterRMatch(hitMatches_neg, clust_RPeaks_neg, rad2);
602 
603 
604  if(clustRMatchIndex == -1) continue;
605 
606  //const double phi2 = reco_pos[j].Phi();
607 
608  // only match pairs that are on the same side of the TPC
609  const bool accepted_z = ((z1>0)==(z2>0));
610  if( !accepted_z ) continue;
611 
612 
613 
614  const bool accepted_r = (hitRadIndex == clustRMatchIndex);
615 
616  const auto dphi = delta_phi(phi1-phi2);
617  const bool accepted_phi = std::abs(dphi) < m_phi_cut;
618 
619  if(!accepted_r || !accepted_phi) continue;
620 
621 
622 
623  if(fabs(dphi)<fabs(prev_dphi)){
624  prev_dphi = dphi;
625  matchJ = j;
626  hits_matched[i] = true;
627  }
628  }//end loop over reco_pos
629 
630  if(matchJ != -1){
631  clusts_matched[matchJ] = true;
632  matched_pair.emplace_back(i,matchJ);
633  matched_nclus.push_back(reco_nclusters[matchJ]);
634 
635  if(m_savehistograms)
636  {
637 
638  const auto& nclus = reco_nclusters[matchJ];
639  const double rad2=get_r(reco_pos[matchJ].X(), reco_pos[matchJ].Y());
640  const double phi2 = reco_pos[matchJ].Phi();
641 
642  const auto dr = rad1-rad2;
643  const auto dphi = delta_phi(phi1-phi2);
644 
645  hnclus->Fill( (float) nclus);
646 
647  double r = rad2;
648 
649  hdrphi->Fill(r * dphi);
650  hdphi->Fill(dphi);
651  hrdphi->Fill(r,dphi);
652  hdrdphi->Fill(dr, dphi);
653  hrdr->Fill(r,dr);
654  if(nclus==1)
655  {
656  if(r < 40.0) hdr1_single->Fill(dr);
657  if(r >= 40.0 && r < 58.0) hdr2_single->Fill(dr);
658  if(r >= 58.0) hdr3_single->Fill(dr);
659  }
660  else
661  {
662  if(r < 40.0) hdr1_double->Fill(dr);
663  if(r >= 40.0 && r < 58.0) hdr2_double->Fill(dr);
664  if(r >= 58.0) hdr3_double->Fill(dr);
665  }
666  }//end save histos
667  }//end if match was found
668  }//end loop over truth pos
669  }//end loop over matching iterations
670 
671  // print some statistics:
672  if( Verbosity() )
673  {
674  const auto n_valid_truth = std::count_if( m_truth_pos.begin(), m_truth_pos.end(), []( const TVector3& pos ) { return get_r( pos.x(), pos.y() ) > 30; } );
675  const auto n_reco_size1 = std::count_if( reco_nclusters.begin(), reco_nclusters.end(), []( const unsigned int& value ) { return value==1; } );
676  const auto n_reco_size2 = std::count_if( reco_nclusters.begin(), reco_nclusters.end(), []( const unsigned int& value ) { return value==2; } );
677  std::cout << "PHTpcCentralMembraneMatcher::process_event - m_truth_pos size: " << m_truth_pos.size() << std::endl;
678  std::cout << "PHTpcCentralMembraneMatcher::process_event - m_truth_pos size, r>30cm: " << n_valid_truth << std::endl;
679  std::cout << "PHTpcCentralMembraneMatcher::process_event - reco_pos size: " << reco_pos.size() << std::endl;
680  std::cout << "PHTpcCentralMembraneMatcher::process_event - reco_pos size (nclus==1): " << n_reco_size1 << std::endl;
681  std::cout << "PHTpcCentralMembraneMatcher::process_event - reco_pos size (nclus==2): " << n_reco_size2 << std::endl;
682  std::cout << "PHTpcCentralMembraneMatcher::process_event - matched_pair size: " << matched_pair.size() << std::endl;
683  }
684 
685  for(unsigned int ip = 0; ip < matched_pair.size(); ++ip)
686  {
687  const std::pair<unsigned int, unsigned int>& p = matched_pair[ip];
688  const unsigned int& nclus = matched_nclus[ip];
689 
690  // add to node tree
691  unsigned int key = p.first;
692  auto cmdiff = new CMFlashDifferencev1();
693  cmdiff->setTruthPhi(m_truth_pos[p.first].Phi());
694  cmdiff->setTruthR(m_truth_pos[p.first].Perp());
695  cmdiff->setTruthZ(m_truth_pos[p.first].Z());
696 
697  cmdiff->setRecoPhi(reco_pos[p.second].Phi());
698  cmdiff->setRecoR(reco_pos[p.second].Perp());
699  cmdiff->setRecoZ(reco_pos[p.second].Z());
700 
701  cmdiff->setNclusters(nclus);
702 
703  m_cm_flash_diffs->addDifferenceSpecifyKey(key, cmdiff);
704 
705  if(m_savehistograms) match_ntup->Fill(m_event_index,m_truth_pos[p.first].Perp(),m_truth_pos[p.first].Phi(),reco_pos[p.second].Perp(),reco_pos[p.second].Phi(),reco_pos[p.second].Z(),nclus,pos1[p.second].Perp(),pos1[p.second].Phi(),adc1[p.second],layer1[p.second],pos2[p.second].Perp(),pos2[p.second].Phi(),adc2[p.second],layer2[p.second]);
706 
707  // store cluster position
708  const double clus_r = reco_pos[p.second].Perp();
709  double clus_phi = reco_pos[p.second].Phi();
710  if ( clus_phi < 0 ) clus_phi += 2*M_PI;
711 
712  const double clus_z = reco_pos[p.second].z();
713  const unsigned int side = (clus_z<0) ? 0:1;
714 
715  // calculate residuals (cluster - truth)
716  const double dr = reco_pos[p.second].Perp() - m_truth_pos[p.first].Perp();
717  const double dphi = delta_phi( reco_pos[p.second].Phi() - m_truth_pos[p.first].Phi() );
718  const double rdphi = reco_pos[p.second].Perp() * dphi;
719  const double dz = reco_pos[p.second].z() - m_truth_pos[p.first].z();
720 
721  // fill distortion correction histograms
722  /*
723  * TODO:
724  * - we might need to only fill the histograms for cm clusters that have 2 clusters only
725  * - we might need a smoothing procedure to fill the bins that have no entries using neighbors
726  */
727  for( const auto& dcc:{m_dcc_out, m_dcc_out_aggregated.get()} )
728  {
729  static_cast<TH2*>(dcc->m_hDRint[side])->Fill( clus_phi, clus_r, dr );
730  static_cast<TH2*>(dcc->m_hDPint[side])->Fill( clus_phi, clus_r, rdphi );
731  static_cast<TH2*>(dcc->m_hDZint[side])->Fill( clus_phi, clus_r, dz );
732  static_cast<TH2*>(dcc->m_hentries[side])->Fill( clus_phi, clus_r );
733  }
734 
735  }
736 
737  if(Verbosity())
738  {
739  std::cout << "PHTpcCentralMembraneMatcher::process_events - cmclusters: " << m_corrected_CMcluster_map->size() << std::endl;
740  std::cout << "PHTpcCentralMembraneMatcher::process_events - matched pairs: " << matched_pair.size() << std::endl;
741  std::cout << "PHTpcCentralMembraneMatcher::process_events - differences: " << m_cm_flash_diffs->size() << std::endl;
742  std::cout << "PHTpcCentralMembraneMatcher::process_events - entries: " << m_dcc_out->m_hentries[0]->GetEntries() << ", " << m_dcc_out->m_hentries[1]->GetEntries() << std::endl;
743  }
744 
745  // normalize per-event distortion correction histograms and fill guarding bins
746  normalize_distortions( m_dcc_out );
747  fill_guarding_bins( m_dcc_out );
748 
749  if(Verbosity())
750  {
751  // read back differences from node tree as a check
752  auto diffrange = m_cm_flash_diffs->getDifferences();
753  for (auto cmitr = diffrange.first;
754  cmitr !=diffrange.second;
755  ++cmitr)
756  {
757  auto key = cmitr->first;
758  auto cmreco = cmitr->second;
759 
760  std::cout << " key " << key
761  << " nclus " << cmreco->getNclusters()
762  << " truth Phi " << cmreco->getTruthPhi() << " reco Phi " << cmreco->getRecoPhi()
763  << " truth R " << cmreco->getTruthR() << " reco R " << cmreco->getRecoR()
764  << " truth Z " << cmreco->getTruthZ() << " reco Z " << cmreco->getRecoZ()
765  << std::endl;
766  }
767  }
768 
769  m_event_index++;
770 
771  return Fun4AllReturnCodes::EVENT_OK;
772 }
773 
774 //____________________________________________________________________________..
775 int PHTpcCentralMembraneMatcher::End(PHCompositeNode * /*topNode*/ )
776 {
777 
778  // write distortion corrections
779  if( m_dcc_out_aggregated )
780  {
781 
782  // normalize aggregated distortion correction histograms and fill guarding bins
783  normalize_distortions( m_dcc_out_aggregated.get() );
784  fill_guarding_bins( m_dcc_out_aggregated.get() );
785 
786  // create TFile and write all histograms
787  std::unique_ptr<TFile> outputfile( TFile::Open( m_outputfile.c_str(), "RECREATE" ) );
788  outputfile->cd();
789 
790  // loop over side
791  for( unsigned int i = 0; i<2; ++i )
792  {
793  for( const auto& h:{m_dcc_out_aggregated->m_hDRint[i], m_dcc_out_aggregated->m_hDPint[i], m_dcc_out_aggregated->m_hDZint[i], m_dcc_out_aggregated->m_hentries[i]} )
794  { if( h ) h->Write(); }
795  }
796  }
797 
798  // write evaluation histograms
799  if(m_savehistograms && fout)
800  {
801  fout->cd();
802 
803  hxy_reco->Write();
804  hxy_truth->Write();
805  hdrdphi->Write();
806  hrdr->Write();
807  hrdphi->Write();
808  hdphi->Write();
809  hdrphi->Write();
810  hdr1_single->Write();
811  hdr2_single->Write();
812  hdr3_single->Write();
813  hdr1_double->Write();
814  hdr2_double->Write();
815  hdr3_double->Write();
816  hnclus->Write();
817 
818  fout->Close();
819  }
820 
821  if(m_savehistograms && fout2)
822  {
823  fout2->cd();
824 
825  match_ntup->Write();
826  hit_r_phi->Write();
827  clust_r_phi_pos->Write();
828  clust_r_phi_neg->Write();
829  }
830 
831  return Fun4AllReturnCodes::EVENT_OK;
832 }
833 
834 //____________________________________________________________________________..
835 
836 int PHTpcCentralMembraneMatcher::GetNodes(PHCompositeNode* topNode)
837 {
838  //---------------------------------
839  // Get Objects off of the Node Tree
840  //---------------------------------
841 
842  m_corrected_CMcluster_map = findNode::getClass<CMFlashClusterContainer>(topNode, "CORRECTED_CM_CLUSTER");
843  if(!m_corrected_CMcluster_map)
844  {
845  std::cout << PHWHERE << "CORRECTED_CM_CLUSTER Node missing, abort." << std::endl;
846  return Fun4AllReturnCodes::ABORTRUN;
847  }
848 
849  // input tpc distortion correction static
850  m_dcc_in_static = findNode::getClass<TpcDistortionCorrectionContainer>(topNode,"TpcDistortionCorrectionContainerStatic");
851  if( m_dcc_in_static )
852  {
853  std::cout << "PHTpcCentralMembraneMatcher: found TPC distortion correction container static" << std::endl;
854  }
855 
856  // input tpc distortion correction average
857  m_dcc_in_average = findNode::getClass<TpcDistortionCorrectionContainer>(topNode,"TpcDistortionCorrectionContainerAverage");
858  if( m_dcc_in_average )
859  {
860  std::cout << "PHTpcCentralMembraneMatcher: found TPC distortion correction container average" << std::endl;
861  }
862 
863  // create node for results of matching
864  std::cout << "Creating node CM_FLASH_DIFFERENCES" << std::endl;
865  PHNodeIterator iter(topNode);
866 
867  // Looking for the DST node
868  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
869  if (!dstNode)
870  {
871  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
872  return Fun4AllReturnCodes::ABORTRUN;
873  }
874  PHNodeIterator dstiter(dstNode);
875  PHCompositeNode *DetNode =
876  dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
877  if (!DetNode)
878  {
879  DetNode = new PHCompositeNode("TRKR");
880  dstNode->addNode(DetNode);
881  }
882 
883  m_cm_flash_diffs = new CMFlashDifferenceContainerv1;
884  PHIODataNode<PHObject> *CMFlashDifferenceNode =
885  new PHIODataNode<PHObject>(m_cm_flash_diffs, "CM_FLASH_DIFFERENCES", "PHObject");
886  DetNode->addNode(CMFlashDifferenceNode);
887 
888 
891  const std::string dcc_out_node_name = "TpcDistortionCorrectionContainerFluctuation";
892  m_dcc_out = findNode::getClass<TpcDistortionCorrectionContainer>(topNode,dcc_out_node_name);
893  if( !m_dcc_out )
894  {
895 
897  auto runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
898  if (!runNode)
899  {
900  std::cout << "PHTpcCentralMembraneMatcher::InitRun - RUN Node missing, quitting" << std::endl;
901  return Fun4AllReturnCodes::ABORTRUN;
902  }
903 
904  std::cout << "PHTpcCentralMembraneMatcher::GetNodes - creating TpcDistortionCorrectionContainer in node " << dcc_out_node_name << std::endl;
905  m_dcc_out = new TpcDistortionCorrectionContainer;
906  auto node = new PHDataNode<TpcDistortionCorrectionContainer>(m_dcc_out, dcc_out_node_name);
907  runNode->addNode(node);
908  }
909 
910 
911  // create per event distortions. Do not put on the node tree
912  //m_dcc_out = new TpcDistortionCorrectionContainer;
913 
914  // also prepare the local distortion container, used to aggregate multple events
915  m_dcc_out_aggregated.reset( new TpcDistortionCorrectionContainer );
916 
917  // compute axis limits to include guarding bins, needed for TH2::Interpolate to work
918  const float phiMin = m_phiMin - (m_phiMax-m_phiMin)/m_phibins;
919  const float phiMax = m_phiMax + (m_phiMax-m_phiMin)/m_phibins;
920 
921  const float rMin = m_rMin - (m_rMax-m_rMin)/m_rbins;
922  const float rMax = m_rMax + (m_rMax-m_rMin)/m_rbins;
923 
924  /*
925  double r_bins_mm[69] = {217.83-2,217.83,
926  223.83, 229.83, 235.83, 241.83, 247.83, 253.83, 259.83, 265.83, 271.83, 277.83, 283.83, 289.83, 295.83, 301.83, 306.83,
927  311.05, 317.92, 323.31, 329.27, 334.63, 340.59, 345.95, 351.91, 357.27, 363.23, 368.59, 374.55, 379.91, 385.87, 391.23, 397.19, 402.49,
928  411.53, 421.70, 431.90, 442.11, 452.32, 462.52, 472.73, 482.94, 493.14, 503.35, 513.56, 523.76, 533.97, 544.18, 554.39, 564.59, 574.76,
929  583.67, 594.59, 605.57, 616.54, 627.51, 638.48, 649.45, 660.42, 671.39, 682.36, 693.33, 704.30, 715.27, 726.24, 737.21, 748.18, 759.11, 759.11+2};
930 
931  double r_bins[69];
932 
933  for(int i=0; i<69; i++){
934  r_bins[i] = r_bins_mm[i]/10.0;
935  }
936 
937 
938 
939  double phiBins[206] = { 0., 6.3083-2 * M_PI, 6.3401-2 * M_PI, 6.372-2 * M_PI, 6.4039-2 * M_PI, 6.4358-2 * M_PI, 6.4676-2 * M_PI, 6.4995-2 * M_PI, 6.5314-2 * M_PI,
940  0.2618, 0.2937, 0.3256, 0.3574, 0.3893, 0.4212, 0.453, 0.4849, 0.5168, 0.5487, 0.5805, 0.6124, 0.6443, 0.6762, 0.7081,
941  0.7399, 0.7718, 0.7854, 0.8173, 0.8491, 0.881, 0.9129, 0.9448, 0.9767, 1.0085, 1.0404, 1.0723, 1.1041, 1.136, 1.1679,
942  1.1998, 1.2317, 1.2635, 1.2954, 1.309, 1.3409, 1.3727, 1.4046, 1.4365, 1.4684, 1.5002, 1.5321, 1.564, 1.5959, 1.6277,
943  1.6596, 1.6915, 1.7234, 1.7552, 1.7871, 1.819, 1.8326, 1.8645, 1.8963, 1.9282, 1.9601, 1.992, 2.0238, 2.0557, 2.0876,
944  2.1195, 2.1513, 2.1832, 2.2151, 2.247, 2.2788, 2.3107, 2.3426, 2.3562, 2.3881, 2.42, 2.4518, 2.4837, 2.5156, 2.5474,
945  2.5793, 2.6112, 2.6431, 2.6749, 2.7068, 2.7387, 2.7706, 2.8024, 2.8343, 2.8662, 2.8798, 2.9117, 2.9436, 2.9754, 3.0073,
946  3.0392, 3.0711, 3.1029, 3.1348, 3.1667, 3.1986, 3.2304, 3.2623, 3.2942, 3.326, 3.3579, 3.3898, 3.4034, 3.4353, 3.4671,
947  3.499, 3.5309, 3.5628, 3.5946, 3.6265, 3.6584, 3.6903, 3.7221, 3.754, 3.7859, 3.8178, 3.8496, 3.8815, 3.9134, 3.927,
948  3.9589, 3.9907, 4.0226, 4.0545, 4.0864, 4.1182, 4.1501, 4.182, 4.2139, 4.2457, 4.2776, 4.3095, 4.3414, 4.3732, 4.4051,
949  4.437, 4.4506, 4.4825, 4.5143, 4.5462, 4.5781, 4.61, 4.6418, 4.6737, 4.7056, 4.7375, 4.7693, 4.8012, 4.8331, 4.865,
950  4.8968, 4.9287, 4.9606, 4.9742, 5.0061, 5.0379, 5.0698, 5.1017, 5.1336, 5.1654, 5.1973, 5.2292, 5.2611, 5.2929, 5.3248,
951  5.3567, 5.3886, 5.4204, 5.4523, 5.4842, 5.4978, 5.5297, 5.5615, 5.5934, 5.6253, 5.6572, 5.689, 5.7209, 5.7528, 5.7847,
952  5.8165, 5.8484, 5.8803, 5.9122, 5.944, 5.9759, 6.0078, 6.0214, 6.0533, 6.0851, 6.117, 6.1489, 6.1808, 6.2127, 6.2445,
953  6.2764, 2 * M_PI};
954  */
955 
956  // reset all output distortion container so that they match the requested grid size
957  const std::array<const std::string,2> extension = {{ "_negz", "_posz" }};
958  for( const auto& dcc:{m_dcc_out, m_dcc_out_aggregated.get()} )
959  {
960  // set dimensions to 2, since central membrane flashes only provide distortions at z = 0
961  dcc->dimensions = 2;
962 
963  // create all histograms
964  for( int i =0; i < 2; ++i )
965  {
966  delete dcc->m_hDPint[i]; dcc->m_hDPint[i] = new TH2F( Form("hIntDistortionP%s", extension[i].c_str()), Form("hIntDistortionP%s", extension[i].c_str()), m_phibins+2, phiMin, phiMax, m_rbins+2, rMin, rMax );
967  delete dcc->m_hDRint[i]; dcc->m_hDRint[i] = new TH2F( Form("hIntDistortionR%s", extension[i].c_str()), Form("hIntDistortionR%s", extension[i].c_str()), m_phibins+2, phiMin, phiMax, m_rbins+2, rMin, rMax );
968  delete dcc->m_hDZint[i]; dcc->m_hDZint[i] = new TH2F( Form("hIntDistortionZ%s", extension[i].c_str()), Form("hIntDistortionZ%s", extension[i].c_str()), m_phibins+2, phiMin, phiMax, m_rbins+2, rMin, rMax );
969  delete dcc->m_hentries[i]; dcc->m_hentries[i] = new TH2I( Form("hEntries%s", extension[i].c_str()), Form("hEntries%s", extension[i].c_str()), m_phibins+2, phiMin, phiMax, m_rbins+2, rMin, rMax );
970 
971 
972  //delete dcc->m_hDPint[i]; dcc->m_hDPint[i] = new TH2F( Form("hIntDistortionP%s", extension[i].c_str()), Form("hIntDistortionP%s", extension[i].c_str()), 205, phiBins, 68, r_bins );
973  //delete dcc->m_hDRint[i]; dcc->m_hDRint[i] = new TH2F( Form("hIntDistortionR%s", extension[i].c_str()), Form("hIntDistortionR%s", extension[i].c_str()), 205, phiBins, 68, r_bins );
974  //delete dcc->m_hDZint[i]; dcc->m_hDZint[i] = new TH2F( Form("hIntDistortionZ%s", extension[i].c_str()), Form("hIntDistortionZ%s", extension[i].c_str()), 205, phiBins, 68, r_bins );
975  //delete dcc->m_hentries[i]; dcc->m_hentries[i] = new TH2I( Form("hEntries%s", extension[i].c_str()), Form("hEntries%s", extension[i].c_str()), 205, phiBins, 68, r_bins);
976  }
977  }
978 
979  return Fun4AllReturnCodes::EVENT_OK;
980 }
981 
982 //_____________________________________________________________
983 void PHTpcCentralMembraneMatcher::CalculateCenters(
984  int nPads,
985  const std::array<double, nRadii>& R,
986  std::array<int, nRadii>& nGoodStripes,
987  const std::array<int, nRadii>& keepUntil,
988  std::array<int, nRadii>& nStripesIn,
989  std::array<int, nRadii>& nStripesBefore,
990  double cx[][nRadii], double cy[][nRadii])
991 {
992  const double phi_module = M_PI / 6.0; // angle span of a module
993  const int pr_mult = 3; // multiples of intrinsic resolution of pads
994  const int dw_mult = 8; // multiples of diffusion width
995  const double diffwidth = 0.6 * mm; // diffusion width
996  const double adjust = 0.015; //arbitrary angle to center the pattern in a petal
997 
998  double theta = 0.0;
999 
1000  //center coords
1001 
1002  //calculate spacing first:
1003  std::array<double, nRadii> spacing;
1004  for (int i = 0; i < nRadii; i++)
1005  {
1006  spacing[i] = 2.0 * ((dw_mult * diffwidth / R[i]) + (pr_mult * phi_module / nPads));
1007  }
1008 
1009  //center calculation
1010  for (int j = 0; j < nRadii; j++)
1011  {
1012  int i_out = 0;
1013  for (int i = keepThisAndAfter[j]; i < keepUntil[j]; i++)
1014  {
1015  if (j % 2 == 0)
1016  {
1017  theta = i * spacing[j] + (spacing[j] / 2) - adjust;
1018  cx[i_out][j] = R[j] * cos(theta) / cm;
1019  cy[i_out][j] = R[j] * sin(theta) / cm;
1020  }
1021  else
1022  {
1023  theta = (i + 1) * spacing[j] - adjust;
1024  cx[i_out][j] = R[j] * cos(theta) / cm;
1025  cy[i_out][j] = R[j] * sin(theta) / cm;
1026  }
1027 
1028  if( Verbosity() > 2 )
1029  std::cout << " j " << j << " i " << i << " i_out " << i_out << " theta " << theta << " cx " << cx[i_out][j] << " cy " << cy[i_out][j]
1030  << " radius " << sqrt(pow(cx[i_out][j],2)+pow(cy[i_out][j],2)) << std::endl;
1031 
1032  i_out++;
1033 
1034  nStripesBefore_R1_e[0] = 0;
1035 
1036  nStripesIn[j] = keepUntil[j] - keepThisAndAfter[j];
1037  if (j==0)
1038  {
1039  nStripesBefore[j] = 0;
1040  }
1041  else
1042  {
1043  nStripesBefore[j] = nStripesIn[j - 1] + nStripesBefore[j - 1];
1044  }
1045  nStripesBefore_R1_e[0] = 0;
1046  }
1047  nGoodStripes[j] = i_out;
1048  }
1049 }
1050 
1051