packages/tpc2019/TpcPrototypeClusterizer.cc

0001 #include "TpcPrototypeClusterizer.h"
0002
0003 #include <tpc/TpcDefs.h>
0004
0005 #include <trackbase/TrkrClusterContainer.h>
0006 #include <trackbase/TrkrClusterHitAssoc.h>
0007 #include <trackbase/TrkrClusterv1.h>
0008 #include <trackbase/TrkrDefs.h>  // for hitkey, getLayer
0009 #include <trackbase/TrkrHit.h>
0010 #include <trackbase/TrkrHitSet.h>
0011 #include <trackbase/TrkrHitSetContainer.h>
0012
0013 #include <fun4all/Fun4AllReturnCodes.h>
0014 #include <fun4all/SubsysReco.h>
0015
0016 #include <g4detectors/PHG4CylinderCellGeom.h>
0017 #include <g4detectors/PHG4CylinderCellGeomContainer.h>
0018
0019 #include <phool/PHCompositeNode.h>
0020 #include <phool/PHIODataNode.h>
0021 #include <phool/PHNode.h>
0022 #include <phool/PHNodeIterator.h>
0023 #include <phool/PHObject.h>
0024 #include <phool/getClass.h>
0025 #include <phool/phool.h>  // for PHWHERE
0026
0027 #include <TMatrixFfwd.h>    // for TMatrixF
0028 #include <TMatrixT.h>       // for TMatrixT, ope...
0029 #include <TMatrixTUtils.h>  // for TMatrixTRow
0030
0031 #include <cmath>  // for sqrt, cos, sin
0032 #include <iostream>
0033 #include <map>  // for _Rb_tree_cons...
0034 #include <string>
0035 #include <utility>  // for pair
0036 #include <vector>
0037
0038 using namespace std;
0039
0040 TpcPrototypeClusterizer::TpcPrototypeClusterizer(const string &name)
0041   : SubsysReco(name)
0042   , m_hits(nullptr)
0043   , m_clusterlist(nullptr)
0044   , m_clusterhitassoc(nullptr)
0045   , zz_shaping_correction(0.0754)
0046   , pedestal(74.4)
0047   , NPhiBinsMax(0)
0048   , NPhiBinsMin(0)
0049   , NZBinsMax(0)
0050   , NZBinsMin(0)
0051 {
0052 }
0053
0054 //===================
0055 bool TpcPrototypeClusterizer::is_local_maximum(int phibin, int zbin, std::vector<std::vector<double>> &adcval)
0056 {
0057   bool retval = true;
0058   double centval = adcval[phibin][zbin];
0059   //cout << "enter is_local_maximum for phibin " << phibin << " zbin " << zbin << " adcval " << centval <<  endl;
0060
0061   // search contiguous adc values for a larger signal
0062   for (int iz = zbin - 4; iz <= zbin + 4; iz++)
0063     for (int iphi = phibin - 2; iphi <= phibin + 2; iphi++)
0064     {
0065       //cout << " is_local_maximum: iphi " <<  iphi << " iz " << iz << " adcval " << adcval[iphi][iz] << endl;
0066       if (iz >= NZBinsMax) continue;
0067       if (iz < NZBinsMin) continue;
0068
0069       if (iphi >= NPhiBinsMax) continue;
0070       if (iphi < NPhiBinsMin) continue;
0071
0072       if (adcval[iphi][iz] > centval)
0073       {
0074         retval = false;
0075       }
0076     }
0077
0078   //if(retval)  cout << "**********************   success: returning " << retval << endl;
0079
0080   return retval;
0081 }
0082
0083 void TpcPrototypeClusterizer::get_cluster(int phibin, int zbin, int &phiup, int &phidown, int &zup, int &zdown, std::vector<std::vector<double>> &adcval)
0084 {
0085   // search along phi at the peak in z
0086
0087   for (int iphi = phibin + 1; iphi < phibin + 4; iphi++)
0088   {
0089     if (iphi >= NPhiBinsMax) continue;
0090
0091     if (adcval[iphi][zbin] <= 0)
0092       break;
0093
0094     if (adcval[iphi][zbin] <= adcval[iphi - 1][zbin])
0095       phiup++;
0096     else
0097       break;
0098   }
0099
0100   for (int iphi = phibin - 1; iphi > phibin - 4; iphi--)
0101   {
0102     if (iphi < NPhiBinsMin) continue;
0103
0104     if (adcval[iphi][zbin] <= 0)
0105       break;
0106
0107     if (adcval[iphi][zbin] <= adcval[iphi + 1][zbin])
0108       phidown++;
0109     else
0110       break;
0111   }
0112
0113   // search along z at the peak in phi
0114
0115   for (int iz = zbin + 1; iz < zbin + 10; iz++)
0116   {
0117     if (iz >= NZBinsMax) continue;
0118
0119     if (adcval[phibin][iz] <= 0)
0120       break;
0121
0122     if (adcval[phibin][iz] <= adcval[phibin][iz - 1])
0123       zup++;
0124     else
0125       break;
0126   }
0127
0128   for (int iz = zbin - 1; iz > zbin - 10; iz--)
0129   {
0130     if (iz < NZBinsMin) continue;
0131
0132     if (adcval[phibin][iz] <= 0)
0133       break;
0134
0135     if (adcval[phibin][iz] <= adcval[phibin][iz + 1])
0136       zdown++;
0137     else
0138       break;
0139   }
0140 }
0141
0142 int TpcPrototypeClusterizer::InitRun(PHCompositeNode *topNode)
0143 {
0144   PHNodeIterator iter(topNode);
0145
0146   // Looking for the DST node
0147   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0148   if (!dstNode)
0149   {
0150     cout << PHWHERE << "DST Node missing, doing nothing." << endl;
0151     return Fun4AllReturnCodes::ABORTRUN;
0152   }
0153
0154   // Create the Cluster node if required
0155   TrkrClusterContainer *trkrclusters = findNode::getClass<TrkrClusterContainer>(dstNode, "TRKR_CLUSTER");
0156   if (!trkrclusters)
0157   {
0158     PHNodeIterator dstiter(dstNode);
0159     PHCompositeNode *DetNode =
0160         dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0161     if (!DetNode)
0162     {
0163       DetNode = new PHCompositeNode("TRKR");
0164       dstNode->addNode(DetNode);
0165     }
0166
0167     trkrclusters = new TrkrClusterContainer();
0168     PHIODataNode<PHObject> *TrkrClusterContainerNode =
0169         new PHIODataNode<PHObject>(trkrclusters, "TRKR_CLUSTER", "PHObject");
0170     DetNode->addNode(TrkrClusterContainerNode);
0171   }
0172
0173   TrkrClusterHitAssoc *clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
0174   if (!clusterhitassoc)
0175   {
0176     PHNodeIterator dstiter(dstNode);
0177     PHCompositeNode *DetNode =
0178         dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0179     if (!DetNode)
0180     {
0181       DetNode = new PHCompositeNode("TRKR");
0182       dstNode->addNode(DetNode);
0183     }
0184
0185     clusterhitassoc = new TrkrClusterHitAssoc();
0186     PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(clusterhitassoc, "TRKR_CLUSTERHITASSOC", "PHObject");
0187     DetNode->addNode(newNode);
0188   }
0189
0190   return Fun4AllReturnCodes::EVENT_OK;
0191 }
0192
0193 int TpcPrototypeClusterizer::process_event(PHCompositeNode *topNode)
0194 {
0195   int print_layer = 47;
0196
0197   if (Verbosity() > 1000)
0198     std::cout << "TpcPrototypeClusterizer::Process_Event" << std::endl;
0199
0200   PHNodeIterator iter(topNode);
0201   PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0202   if (!dstNode)
0203   {
0204     cout << PHWHERE << "DST Node missing, doing nothing." << endl;
0205     return Fun4AllReturnCodes::ABORTRUN;
0206   }
0207
0208   // get node containing the digitized hits
0209   m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0210   if (!m_hits)
0211   {
0212     cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << endl;
0213     return Fun4AllReturnCodes::ABORTRUN;
0214   }
0215
0216   // get node for clusters
0217   m_clusterlist = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0218   if (!m_clusterlist)
0219   {
0220     cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTER." << endl;
0221     return Fun4AllReturnCodes::ABORTRUN;
0222   }
0223
0224   // get node for cluster hit associations
0225   m_clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
0226   if (!m_clusterhitassoc)
0227   {
0228     cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTERHITASSOC" << endl;
0229     return Fun4AllReturnCodes::ABORTRUN;
0230   }
0231
0232   PHG4CylinderCellGeomContainer *geom_container =
0233       findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0234   if (!geom_container)
0235   {
0236     std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
0237     return Fun4AllReturnCodes::ABORTRUN;
0238   }
0239
0240   // The hits are stored in hitsets, where each hitset contains all hits in a given TPC readout (layer, sector, side), so clusters are confined to a hitset
0241   // The TPC clustering is more complicated than for the silicon, because we have to deal with overlapping clusters
0242
0243   // loop over the TPC HitSet objects
0244   TrkrHitSetContainer::ConstRange hitsetrange = m_hits->getHitSets(TrkrDefs::TrkrId::tpcId);
0245   for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
0246        hitsetitr != hitsetrange.second;
0247        ++hitsetitr)
0248   {
0249     TrkrHitSet *hitset = hitsetitr->second;
0250     if (Verbosity() > 1)
0251       cout << "TpcPrototypeClusterizer process hitsetkey " << hitsetitr->first
0252            << " layer " << (int) TrkrDefs::getLayer(hitsetitr->first)
0253            << " side " << (int) TpcDefs::getSide(hitsetitr->first)
0254            << " sector " << (int) TpcDefs::getSectorId(hitsetitr->first)
0255            << endl;
0256     if (Verbosity() > 2) hitset->identify();
0257
0258     // we have a single hitset, get the info that identifies the module
0259     int layer = TrkrDefs::getLayer(hitsetitr->first);
0260     // int sector = TpcDefs::getSectorId(hitsetitr->first);
0261     int side = TpcDefs::getSide(hitsetitr->first);
0262
0263     // we will need the geometry object for this layer to get the global position
0264     PHG4CylinderCellGeom *layergeom = geom_container->GetLayerCellGeom(layer);
0265     int NPhiBins = layergeom->get_phibins();
0266     NPhiBinsMin = 0;
0267     NPhiBinsMax = NPhiBins;
0268
0269     int NZBins = layergeom->get_zbins();
0270     if (side == 0)
0271     {
0272       NZBinsMin = 0;
0273       NZBinsMax = NZBins / 2;
0274     }
0275     else
0276     {
0277       NZBinsMin = NZBins / 2 + 1;
0278       NZBinsMax = NZBins;
0279     }
0280
0281     // for convenience, create a 2D vector to store adc values in and initialize to zero
0282     std::vector<std::vector<double>> adcval(NPhiBins, std::vector<double>(NZBins, 0));
0283
0284     TrkrHitSet::ConstRange hitrangei = hitset->getHits();
0285     for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0286          hitr != hitrangei.second;
0287          ++hitr)
0288     {
0289       int phibin = TpcDefs::getPad(hitr->first);
0290       int zbin = TpcDefs::getTBin(hitr->first);
0291       if (hitr->second->getAdc() > 0)
0292         adcval[phibin][zbin] = (double) hitr->second->getAdc() - pedestal;
0293
0294       if (Verbosity() > 2)
0295         //        if (layer == print_layer)
0296         cout << " add hit in layer " << layer << " with phibin " << phibin << " zbin " << zbin << " adcval " << hitr->second->getAdc() << "->" << adcval[phibin][zbin] << endl;
0297     }
0298
0299     vector<int> phibinlo;
0300     vector<int> phibinhi;
0301     vector<int> zbinlo;
0302     vector<int> zbinhi;
0303     // we want to search the hit list for local maxima in phi-z space and cluster around them
0304     for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0305          hitr != hitrangei.second;
0306          ++hitr)
0307     {
0308       int phibin = TpcDefs::getPad(hitr->first);
0309       int zbin = TpcDefs::getTBin(hitr->first);
0310       if (!is_local_maximum(phibin, zbin, adcval)) continue;
0311
0312       int phiup = 0;
0313       int phidown = 0;
0314       int zup = 0;
0315       int zdown = 0;
0316       // cluster the hits around this local maximum
0317       get_cluster(phibin, zbin, phiup, phidown, zup, zdown, adcval);
0318
0319       if (phiup == 0 && phidown == 0 && zup == 0 and zdown == 0) continue;  // ignore isolated noise hit
0320
0321       // Add this cluster to a vector of clusters for later analysis
0322       phibinlo.push_back(phibin - phidown);
0323       phibinhi.push_back(phibin + phiup);
0324       zbinlo.push_back(zbin - zdown);
0325       zbinhi.push_back(zbin + zup);
0326
0327       if (Verbosity() > 2)
0328         //        if (layer == print_layer)
0329         cout << " cluster found in layer " << layer << " around hitkey " << hitr->first << " with zbin " << zbin << " zup " << zup << " zdown " << zdown
0330              << " phibin " << phibin << " phiup " << phiup << " phidown " << phidown << endl;
0331     }  // end loop over hits in this hitset
0332
0333     // Now we analyze the clusters to get their parameters
0334     for (unsigned int iclus = 0; iclus < phibinlo.size(); iclus++)
0335     {
0336       //cout << "TpcPrototypeClusterizer: process cluster iclus = " << iclus <<  " in layer " << layer << endl;
0337       // loop over the hits in this cluster
0338       double zsum = 0.0;
0339       double phi_sum = 0.0;
0340       double adc_sum = 0.0;
0341       double radius = layergeom->get_radius();  // returns center of layer
0342       if (Verbosity() > 2)
0343         if (layer == print_layer)
0344         {
0345           cout << "iclus " << iclus << " layer " << layer << " radius " << radius << endl;
0346           cout << "    z bin range " << zbinlo[iclus] << " to " << zbinhi[iclus] << " phibin range " << phibinlo[iclus] << " to " << phibinhi[iclus] << endl;
0347         }
0348
0349       std::vector<TrkrDefs::hitkey> hitkeyvec;
0350       for (int iphi = phibinlo[iclus]; iphi <= phibinhi[iclus]; iphi++)
0351       {
0352         for (int iz = zbinlo[iclus]; iz <= zbinhi[iclus]; iz++)
0353         {
0354           double phi_center = layergeom->get_phicenter(iphi);
0355           double z = layergeom->get_zcenter(iz);
0356
0357           zsum += z * adcval[iphi][iz];
0358           phi_sum += phi_center * adcval[iphi][iz];
0359           adc_sum += adcval[iphi][iz];
0360
0361           // capture the hitkeys for all non-zero adc values
0362           if (adcval[iphi][iz] != 0)
0363           {
0364             TrkrDefs::hitkey hitkey = TpcDefs::genHitKey(iphi, iz);
0365             hitkeyvec.push_back(hitkey);
0366           }
0367         }
0368       }
0369
0370       if (adc_sum < 10) continue;  // skip obvious noise "clusters"
0371
0372       // This is the global position
0373       double clusphi = phi_sum / adc_sum;
0374       double clusz = zsum / adc_sum;
0375
0376       // create the cluster entry directly in the node tree
0377       TrkrDefs::cluskey ckey = TpcDefs::genClusKey(hitset->getHitSetKey(), iclus);
0378       TrkrClusterv1 *clus = static_cast<TrkrClusterv1 *>((m_clusterlist->findOrAddCluster(ckey))->second);
0379
0380       double phi_size = (double) (phibinhi[iclus] - phibinlo[iclus] + 1) * radius * layergeom->get_phistep();
0381       double z_size = (double) (zbinhi[iclus] - zbinlo[iclus] + 1) * layergeom->get_zstep();
0382
0383       // Estimate the errors
0384       double dphi2_adc = 0.0;
0385       double dphi_adc = 0.0;
0386       double dz2_adc = 0.0;
0387       double dz_adc = 0.0;
0388       for (int iz = zbinlo[iclus]; iz <= zbinhi[iclus]; iz++)
0389       {
0390         for (int iphi = phibinlo[iclus]; iphi <= phibinhi[iclus]; iphi++)
0391         {
0392           double dphi = layergeom->get_phicenter(iphi) - clusphi;
0393           dphi2_adc += dphi * dphi * adcval[iphi][iz];
0394           dphi_adc += dphi * adcval[iphi][iz];
0395
0396           double dz = layergeom->get_zcenter(iz) - clusz;
0397           dz2_adc += dz * dz * adcval[iphi][iz];
0398           dz_adc += dz * adcval[iphi][iz];
0399         }
0400       }
0401       double phi_cov = (dphi2_adc / adc_sum - dphi_adc * dphi_adc / (adc_sum * adc_sum));
0402       double z_cov = dz2_adc / adc_sum - dz_adc * dz_adc / (adc_sum * adc_sum);
0403
0404       //cout << " layer " << layer << " z_cov " << z_cov << " dz2_adc " << dz2_adc << " adc_sum " <<  adc_sum << " dz_adc " << dz_adc << endl;
0405
0406       // phi_cov = (weighted mean of dphi^2) - (weighted mean of dphi)^2,  which is essentially the weighted mean of dphi^2. The error is then:
0407       // e_phi = sigma_dphi/sqrt(N) = sqrt( sigma_dphi^2 / N )  -- where N is the number of samples of the distribution with standard deviation sigma_dphi
0408       //    - N is the number of electrons that drift to the readout plane
0409       // We have to convert (sum of adc units for all bins in the cluster) to number of ionization electrons N
0410       // Conversion gain is 20 mV/fC - relates total charge collected on pad to PEAK voltage out of ADC. The GEM gain is assumed to be 2000
0411       // To get equivalent charge per Z bin, so that summing ADC input voltage over all Z bins returns total input charge, divide voltages by 2.4 for 80 ns SAMPA
0412       // Equivalent charge per Z bin is then  (ADU x 2200 mV / 1024) / 2.4 x (1/20) fC/mV x (1/1.6e-04) electrons/fC x (1/2000) = ADU x 0.14
0413
0414       double phi_err = radius * sqrt(phi_cov / (adc_sum * 0.14));
0415       if (phi_err == 0.0)  // a single phi bin will cause this
0416         phi_err = radius * layergeom->get_phistep() / sqrt(12.0);
0417
0418       double z_err = sqrt(z_cov / (adc_sum * 0.14));
0419       if (z_err == 0.0)
0420         z_err = layergeom->get_zstep() / sqrt(12.0);
0421
0422       // This corrects the bias introduced by the asymmetric SAMPA chip shaping - assumes 80 ns shaping time
0423       if (clusz < 0)
0424         clusz -= zz_shaping_correction;
0425       else
0426         clusz += zz_shaping_correction;
0427
0428       // Fill in the cluster details
0429       //================
0430       clus->setAdc(adc_sum);
0431       clus->setPosition(0, radius * cos(clusphi));
0432       clus->setPosition(1, radius * sin(clusphi));
0433       clus->setPosition(2, clusz);
0434       clus->setGlobal();
0435
0436       TMatrixF DIM(3, 3);
0437       DIM[0][0] = 0.0;
0438       DIM[0][1] = 0.0;
0439       DIM[0][2] = 0.0;
0440       DIM[1][0] = 0.0;
0441       DIM[1][1] = phi_size / radius * phi_size / radius;  //cluster_v1 expects polar coordinates covariance
0442       DIM[1][2] = 0.0;
0443       DIM[2][0] = 0.0;
0444       DIM[2][1] = 0.0;
0445       DIM[2][2] = z_size * z_size;
0446
0447       TMatrixF ERR(3, 3);
0448       ERR[0][0] = 0.0;
0449       ERR[0][1] = 0.0;
0450       ERR[0][2] = 0.0;
0451       ERR[1][0] = 0.0;
0452       ERR[1][1] = phi_err * phi_err;  //cluster_v1 expects rad, arc, z as elementsof covariance
0453       ERR[1][2] = 0.0;
0454       ERR[2][0] = 0.0;
0455       ERR[2][1] = 0.0;
0456       ERR[2][2] = z_err * z_err;
0457
0458       TMatrixF ROT(3, 3);
0459       ROT[0][0] = cos(clusphi);
0460       ROT[0][1] = -sin(clusphi);
0461       ROT[0][2] = 0.0;
0462       ROT[1][0] = sin(clusphi);
0463       ROT[1][1] = cos(clusphi);
0464       ROT[1][2] = 0.0;
0465       ROT[2][0] = 0.0;
0466       ROT[2][1] = 0.0;
0467       ROT[2][2] = 1.0;
0468
0469       TMatrixF ROT_T(3, 3);
0470       ROT_T.Transpose(ROT);
0471
0472       TMatrixF COVAR_DIM(3, 3);
0473       COVAR_DIM = ROT * DIM * ROT_T;
0474
0475       clus->setSize(0, 0, COVAR_DIM[0][0]);
0476       clus->setSize(0, 1, COVAR_DIM[0][1]);
0477       clus->setSize(0, 2, COVAR_DIM[0][2]);
0478       clus->setSize(1, 0, COVAR_DIM[1][0]);
0479       clus->setSize(1, 1, COVAR_DIM[1][1]);
0480       clus->setSize(1, 2, COVAR_DIM[1][2]);
0481       clus->setSize(2, 0, COVAR_DIM[2][0]);
0482       clus->setSize(2, 1, COVAR_DIM[2][1]);
0483       clus->setSize(2, 2, COVAR_DIM[2][2]);
0484       //cout << " covar_dim[2][2] = " <<  COVAR_DIM[2][2] << endl;
0485
0486       TMatrixF COVAR_ERR(3, 3);
0487       COVAR_ERR = ROT * ERR * ROT_T;
0488
0489       clus->setError(0, 0, COVAR_ERR[0][0]);
0490       clus->setError(0, 1, COVAR_ERR[0][1]);
0491       clus->setError(0, 2, COVAR_ERR[0][2]);
0492       clus->setError(1, 0, COVAR_ERR[1][0]);
0493       clus->setError(1, 1, COVAR_ERR[1][1]);
0494       clus->setError(1, 2, COVAR_ERR[1][2]);
0495       clus->setError(2, 0, COVAR_ERR[2][0]);
0496       clus->setError(2, 1, COVAR_ERR[2][1]);
0497       clus->setError(2, 2, COVAR_ERR[2][2]);
0498
0499       /*
0500       for(int i=0;i<3;i++)
0501         for(int j=0;j<3;j++)
0502           cout << "    i " << i << " j " << j << " clusphi " << clusphi << " clus error " << clus->getError(i,j)
0503            << " clus size " << clus->getSize(i,j) << " ROT " << ROT[i][j] << " ROT_T " << ROT_T[i][j] << " COVAR_ERR " << COVAR_ERR[i][j] << endl;
0504       */
0505
0506       // Add the hit associations to the TrkrClusterHitAssoc node
0507       // we need the cluster key and all associated hit keys (note: the cluster key includes the hitset key)
0508       for (unsigned int i = 0; i < hitkeyvec.size(); i++)
0509       {
0510         m_clusterhitassoc->addAssoc(ckey, hitkeyvec[i]);
0511       }
0512
0513     }  // end loop over clusters for this hitset
0514   }    // end loop over hitsets
0515
0516   if (Verbosity() > 2)
0517   {
0518     cout << "Dump clusters after TpcPrototypeClusterizer" << endl;
0519     m_clusterlist->identify();
0520   }
0521
0522   if (Verbosity() > 2)
0523   {
0524     cout << "Dump cluster hit associations after TpcPrototypeClusterizer" << endl;
0525     m_clusterhitassoc->identify();
0526   }
0527
0528   return Fun4AllReturnCodes::EVENT_OK;
0529 }