Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-16 09:20:22 

0001 #include "LaserClusterizer.h"
0002 
0003 #include "LaserEventInfo.h"
0004 
0005 #include <trackbase/LaserCluster.h>
0006 #include <trackbase/LaserClusterContainer.h>
0007 #include <trackbase/LaserClusterContainerv1.h>
0008 #include <trackbase/LaserClusterv2.h>
0009 #include <trackbase/TpcDefs.h>
0010 #include <trackbase/TrkrDefs.h>  // for hitkey, getLayer
0011 #include <trackbase/TrkrHit.h>
0012 #include <trackbase/TrkrHitSet.h>
0013 #include <trackbase/TrkrHitSetContainer.h>
0014 
0015 #include <ffaobjects/EventHeader.h>
0016 #include <fun4all/Fun4AllReturnCodes.h>
0017 #include <fun4all/SubsysReco.h>  // for SubsysReco
0018 
0019 #include <g4detectors/PHG4TpcGeom.h>
0020 #include <g4detectors/PHG4TpcGeomContainer.h>
0021 
0022 #include <phool/PHCompositeNode.h>
0023 #include <phool/PHIODataNode.h>  // for PHIODataNode
0024 #include <phool/PHNode.h>        // for PHNode
0025 #include <phool/PHNodeIterator.h>
0026 #include <phool/PHObject.h>  // for PHObject
0027 #include <phool/PHTimer.h>
0028 #include <phool/getClass.h>
0029 #include <phool/phool.h>  // for PHWHERE
0030 
0031 #include <Acts/Definitions/Units.hpp>
0032 #include <Acts/Surfaces/Surface.hpp>
0033 
0034 #include <TF1.h>
0035 #include <TFile.h>
0036 #include <TH3.h>
0037 //#include <TF3.h>
0038 //#include <ROOT/Fit/Fitter.h>
0039 #include <Fit/Fitter.h>
0040 //#include <Math/Functor.h>
0041 
0042 #include <boost/geometry.hpp>
0043 #include <boost/geometry/geometries/box.hpp>
0044 #include <boost/geometry/geometries/point.hpp>
0045 #include <boost/geometry/index/rtree.hpp>
0046 
0047 #include <array>
0048 #include <cmath>  // for sqrt, cos, sin
0049 #include <iostream>
0050 #include <limits>
0051 #include <map>  // for _Rb_tree_cons...
0052 #include <string>
0053 #include <utility>  // for pair
0054 #include <vector>
0055 #include <set>
0056 #include <queue>
0057 #include <tuple>
0058 
0059 #include <pthread.h>
0060 
0061 namespace bg = boost::geometry;
0062 namespace bgi = boost::geometry::index;
0063 
0064 using point = bg::model::point<float, 3, bg::cs::cartesian>;
0065 using box = bg::model::box<point>;
0066 using specHitKey = std::pair<TrkrDefs::hitkey, TrkrDefs::hitsetkey>;
0067 using adcKey = std::pair<unsigned int, specHitKey>;
0068 using pointKeyLaser = std::pair<point, specHitKey>;
0069 using hitData = std::pair<point, adcKey>;
0070 
0071 
0072 int layerMins[3] = {7,23,39};
0073 int layerMaxes[3] = {22, 38, 54};
0074 
0075 
0076 
0077 namespace
0078 {
0079   struct thread_data
0080   {
0081     PHG4TpcGeomContainer *geom_container = nullptr;
0082     ActsGeometry *tGeometry = nullptr;
0083     std::vector<TrkrHitSet *> hitsets;
0084     std::vector<unsigned int> layers;
0085     bool side = false;
0086     unsigned int sector = 0;
0087     unsigned int module = 0;
0088     std::vector<LaserCluster *> cluster_vector;
0089     std::vector<TrkrDefs::cluskey> cluster_key_vector;
0090     double adc_threshold = 74.4;
0091     int peakTimeBin = 325;
0092     int layerMin = 1;
0093     int layerMax = 1;
0094     double tdriftmax = 0;
0095     int eventNum = 0;
0096     int Verbosity = 0;
0097     TH3D *hitHist = nullptr;
0098     bool doFitting = false;
0099   };
0100 
0101   pthread_mutex_t mythreadlock;
0102 
0103   const std::vector<point> neighborOffsets = {
0104     point(1, 0, 0), point(-1, 0, 0),
0105     point(0, 1, 0), point(0, -1, 0),
0106     point(0, 0, 1), point(0, 0, -1),
0107     point(0, 0, 2), point(0, 0, -2)
0108   };
0109 
0110 
0111   double layerFunction(double *x, double *par)
0112   {
0113     double A = par[0];
0114     double mu = par[1];
0115 
0116 
0117     double binCenter = round(x[0]);
0118     double overlapLow = std::max(binCenter - 0.5, mu - 0.5);
0119     double overlapHigh = std::min(binCenter + 0.5, mu + 0.5);
0120     double overlap = overlapHigh - overlapLow;
0121     if(overlap <= 0.0)
0122     {
0123       return 0.0;
0124     }
0125     return A*overlap;
0126     /*
0127     if(fabs(x[0] - mu) < 1)
0128     {
0129       double frac = 1-fabs(mu - round(x[0]));
0130       return A*(frac);
0131     }
0132     return 0.0;
0133     */
0134 
0135 
0136   }
0137 
0138   double phiFunction(double *x, double *par)
0139   {
0140     if(par[2] < 0.0)
0141     {
0142       return 0.0;
0143     }
0144     return par[0] * TMath::Gaus(x[0],par[1],par[2],false);
0145   }
0146 
0147   double timeFunction(double *x, double *par)
0148   {
0149     if(par[2] < 0.0)
0150     {
0151       return 0.0;
0152     }
0153     double g = TMath::Gaus(x[0],par[1],par[2],true);
0154     double cdf = 1 + TMath::Erfc(par[3]*(x[0]-par[1])/(sqrt(2.0)*par[2]));
0155     return par[0]*g*cdf;
0156   }
0157 
0158   void findConnectedRegions3(std::vector<hitData> &clusHits, std::pair<TrkrDefs::hitkey, TrkrDefs::hitsetkey> &maxKey)
0159   {
0160     std::vector<std::vector<hitData>> regions;
0161 
0162     std::vector<hitData> unvisited;
0163     for(auto &clusHit : clusHits)
0164     {
0165       unvisited.push_back(clusHit);
0166     }
0167 
0168     while(!unvisited.empty())
0169     {
0170       std::vector<hitData> region;
0171       std::queue<hitData> q;
0172 
0173       unsigned int mIndex = 0;
0174       int i=0;
0175       for(auto hit : unvisited)
0176       {
0177         if(hit.second.second.first == maxKey.first && hit.second.second.second == maxKey.second)
0178         {
0179           mIndex = i;
0180           break;
0181         }
0182         i++;
0183       }
0184 
0185       auto seed = unvisited[mIndex];
0186       unvisited.erase(unvisited.begin()+mIndex);
0187       q.push(seed);
0188       region.push_back(seed);
0189 
0190       while(!q.empty())
0191       {
0192         float ix = q.front().first.get<0>();
0193         float iy = q.front().first.get<1>();
0194         float iz = q.front().first.get<2>();
0195         q.pop();
0196 
0197         for (auto neigh : neighborOffsets)
0198         {
0199           float nx = ix + neigh.get<0>();
0200           float ny = iy + neigh.get<1>();
0201           float nz = iz + neigh.get<2>();
0202 
0203           for(unsigned int v=0; v<unvisited.size(); v++)
0204           {
0205 
0206             if(fabs(unvisited[v].first.get<0>() - nx) < 0.01 && fabs(unvisited[v].first.get<1>() - ny) < 0.01 && fabs(unvisited[v].first.get<2>() - nz) < 0.01)
0207             {
0208               auto newSeed = unvisited[v];
0209               unvisited.erase(unvisited.begin()+v);
0210               q.push(newSeed);
0211               region.push_back(newSeed);
0212               break;
0213             }
0214           }
0215         }
0216       }
0217       regions.push_back(region);
0218 
0219     }
0220 
0221     clusHits.clear();
0222     for(auto hit : regions[0])
0223     {
0224       clusHits.push_back(hit);
0225     }
0226 
0227   }
0228 
0229    
0230   void remove_hits(std::vector<hitData> &clusHits, bgi::rtree<hitData, bgi::quadratic<16>> &rtree, std::multimap<unsigned int, pointKeyLaser> &adcMap)
0231   {
0232     for (auto &clusHit : clusHits)
0233     {
0234       auto spechitkey = clusHit.second.second;
0235 
0236       rtree.remove(clusHit);
0237 
0238       for (auto iterAdc = adcMap.begin(); iterAdc != adcMap.end();)
0239       {
0240           if(iterAdc->second.second == spechitkey)
0241           {
0242             iterAdc = adcMap.erase(iterAdc);
0243             break;
0244           }
0245           else
0246           {
0247             ++iterAdc;
0248           }
0249       }
0250     }
0251 
0252   }
0253 
0254   void calc_cluster_parameter(std::vector<hitData> &clusHits, thread_data &my_data, std::pair<TrkrDefs::hitkey, TrkrDefs::hitsetkey> maxADCKey)
0255   {
0256 
0257 
0258 
0259     findConnectedRegions3(clusHits, maxADCKey);
0260 
0261 
0262     double rSum = 0.0;
0263     double phiSum = 0.0;
0264     double tSum = 0.0;
0265     
0266     double layerSum = 0.0;
0267     double iphiSum = 0.0;
0268     double itSum = 0.0;
0269     
0270     double adcSum = 0.0;
0271     
0272     double maxAdc = 0.0;
0273     TrkrDefs::hitsetkey maxKey = 0;
0274     
0275     unsigned int nHits = clusHits.size();
0276     
0277     auto *clus = new LaserClusterv2;
0278     
0279     int meanSide = 0;
0280     
0281     std::vector<float> usedLayer;
0282     std::vector<float> usedIPhi;
0283     std::vector<float> usedIT;
0284     
0285     double meanLayer = 0.0;
0286     double meanIPhi = 0.0;
0287     double meanIT = 0.0;
0288 
0289     for (auto &clusHit : clusHits)
0290     {
0291       float coords[3] = {clusHit.first.get<0>(), clusHit.first.get<1>(), clusHit.first.get<2>()};
0292       std::pair<TrkrDefs::hitkey, TrkrDefs::hitsetkey> spechitkey = clusHit.second.second;
0293       unsigned int adc = clusHit.second.first;
0294 
0295       int side = TpcDefs::getSide(spechitkey.second);
0296           
0297       if (side)
0298       {
0299           meanSide++;
0300       }
0301       else
0302       {
0303           meanSide--;
0304       }
0305       
0306       PHG4TpcGeom *layergeom = my_data.geom_container->GetLayerCellGeom((int) coords[0]);
0307       
0308       double r = layergeom->get_radius();
0309       double phi = layergeom->get_phi(coords[1], side);
0310       double t = layergeom->get_zcenter(fabs(coords[2]));
0311       
0312       double hitzdriftlength = t * my_data.tGeometry->get_drift_velocity();
0313       double hitZ = my_data.tdriftmax * my_data.tGeometry->get_drift_velocity() - hitzdriftlength;
0314       
0315       
0316         bool foundLayer = false;
0317         for (float i : usedLayer)
0318         {
0319           if (coords[0] == i)
0320           {
0321             foundLayer = true;
0322             break;
0323           }
0324         }
0325       
0326         if (!foundLayer)
0327         {
0328           usedLayer.push_back(coords[0]);
0329         }
0330       
0331         bool foundIPhi = false;
0332         for (float i : usedIPhi)
0333         {
0334           if (coords[1] == i)
0335           {
0336             foundIPhi = true;
0337             break;
0338           }
0339         }
0340       
0341         if (!foundIPhi)
0342         {
0343           usedIPhi.push_back(coords[1]);
0344         }
0345       
0346         bool foundIT = false;
0347         for (float i : usedIT)
0348         {
0349           if (coords[2] == i)
0350           {
0351             foundIT = true;
0352             break;
0353           }
0354         }
0355       
0356         if (!foundIT)
0357         {
0358           usedIT.push_back(coords[2]);
0359         }
0360       
0361         clus->addHit();
0362         clus->setHitLayer(clus->getNhits() - 1, coords[0]);
0363         clus->setHitIPhi(clus->getNhits() - 1, coords[1]);
0364         clus->setHitIT(clus->getNhits() - 1, coords[2]);
0365         clus->setHitX(clus->getNhits() - 1, r * cos(phi));
0366         clus->setHitY(clus->getNhits() - 1, r * sin(phi));
0367         clus->setHitZ(clus->getNhits() - 1, hitZ);
0368         clus->setHitAdc(clus->getNhits() - 1, (float) adc);
0369       
0370         rSum += r * adc;
0371         phiSum += phi * adc;
0372         tSum += t * adc;
0373       
0374         layerSum += coords[0] * adc;
0375         iphiSum += coords[1] * adc;
0376         itSum += coords[2] * adc;
0377       
0378         meanLayer += coords[0];
0379         meanIPhi += coords[1];
0380         meanIT += coords[2];
0381       
0382         adcSum += adc;
0383       
0384         if (adc > maxAdc)
0385         {
0386           maxAdc = adc;
0387           maxKey = spechitkey.second;
0388         }
0389       
0390       }
0391     
0392     if (nHits == 0)
0393     {
0394       return;
0395     }
0396 
0397     
0398     double clusR = rSum / adcSum;
0399     double clusPhi = phiSum / adcSum;
0400     double clusT = tSum / adcSum;
0401     double zdriftlength = clusT * my_data.tGeometry->get_drift_velocity();
0402     
0403     double clusX = clusR * cos(clusPhi);
0404     double clusY = clusR * sin(clusPhi);
0405     double clusZ = my_data.tdriftmax * my_data.tGeometry->get_drift_velocity() - zdriftlength;
0406     if (meanSide < 0)
0407     {
0408       clusZ = -clusZ;
0409       for (int i = 0; i < (int) clus->getNhits(); i++)
0410       {
0411           clus->setHitZ(i, -1 * clus->getHitZ(i));
0412       }
0413     }
0414     
0415     std::sort(usedLayer.begin(), usedLayer.end());
0416     std::sort(usedIPhi.begin(), usedIPhi.end());
0417     std::sort(usedIT.begin(), usedIT.end());
0418 
0419     meanLayer = meanLayer / nHits;
0420     meanIPhi = meanIPhi / nHits;
0421     meanIT = meanIT / nHits;
0422     
0423     double sigmaLayer = 0.0;
0424     double sigmaIPhi = 0.0;
0425     double sigmaIT = 0.0;
0426     
0427     double sigmaWeightedLayer = 0.0;
0428     double sigmaWeightedIPhi = 0.0;
0429     double sigmaWeightedIT = 0.0;
0430     
0431     pthread_mutex_lock(&mythreadlock);
0432     my_data.hitHist = new TH3D(Form("hitHist_event%d_side%d_sector%d_module%d_cluster%d",my_data.eventNum,(int)my_data.side,(int)my_data.sector,(int)my_data.module,(int)my_data.cluster_vector.size()),";layer;iphi;it",usedLayer.size()+2,usedLayer[0]-1.5,*usedLayer.rbegin()+1.5,usedIPhi.size()+2,usedIPhi[0]-1.5,*usedIPhi.rbegin()+1.5,usedIT.size()+2,usedIT[0]-1.5,*usedIT.rbegin()+1.5);
0433 
0434 
0435     //TH3D *hitHist = new TH3D(Form("hitHist_event%d_side%d_sector%d_module%d_cluster%d",my_data.eventNum,(int)my_data.side,(int)my_data.sector,(int)my_data.module,(int)my_data.cluster_vector.size()),";layer;iphi;it",usedLayer.size()+2,usedLayer[0]-1.5,*usedLayer.rbegin()+1.5,usedIPhi.size()+2,usedIPhi[0]-1.5,*usedIPhi.rbegin()+1.5,usedIT.size()+2,usedIT[0]-1.5,*usedIT.rbegin()+1.5);
0436 
0437     for (int i = 0; i < (int) clus->getNhits(); i++)
0438     {
0439 
0440       my_data.hitHist->Fill(clus->getHitLayer(i), clus->getHitIPhi(i), clus->getHitIT(i), clus->getHitAdc(i));
0441 
0442       sigmaLayer += pow(clus->getHitLayer(i) - meanLayer, 2);
0443       sigmaIPhi += pow(clus->getHitIPhi(i) - meanIPhi, 2);
0444       sigmaIT += pow(clus->getHitIT(i) - meanIT, 2);
0445       
0446       sigmaWeightedLayer += clus->getHitAdc(i) * pow(clus->getHitLayer(i) - (layerSum / adcSum), 2);
0447       sigmaWeightedIPhi += clus->getHitAdc(i) * pow(clus->getHitIPhi(i) - (iphiSum / adcSum), 2);
0448       sigmaWeightedIT += clus->getHitAdc(i) * pow(clus->getHitIT(i) - (itSum / adcSum), 2);
0449     }
0450 
0451     bool fitSuccess = false;
0452     ROOT::Fit::Fitter *fit3D = new ROOT::Fit::Fitter;
0453 
0454     if(my_data.doFitting)
0455     {
0456 
0457       double par_init[7] = {
0458         maxAdc,
0459         meanLayer,
0460         meanIPhi, 0.75,
0461         meanIT, 0.5, 1
0462       };
0463 
0464       double satThreshold = 900.0;
0465       double sigma_ADC = 20.0;
0466 
0467       auto nll = [&](const double* par)
0468       {
0469         double nll_val = 0.0;
0470 
0471         int nx = my_data.hitHist->GetNbinsX();
0472         int ny = my_data.hitHist->GetNbinsY();
0473         int nz = my_data.hitHist->GetNbinsZ();
0474 
0475         double parLayer[2] = {1.0,par[1]};
0476         double parPhi[4] = {1.0,par[2],par[3]};
0477         double parTime[4] = {1.0,par[4],par[5],par[6]};
0478 
0479         double xyz[3];
0480 
0481         for (int i = 1; i <= nx; ++i)
0482         {
0483           xyz[0] = my_data.hitHist->GetXaxis()->GetBinCenter(i);  
0484           for (int j = 1; j <= ny; ++j)
0485           {
0486             xyz[1] = my_data.hitHist->GetYaxis()->GetBinCenter(j);    
0487             for (int k = 1; k <= nz; ++k)
0488             {
0489               xyz[2] = my_data.hitHist->GetZaxis()->GetBinCenter(k);
0490               double observed = my_data.hitHist->GetBinContent(i, j, k);
0491 
0492               double expected = par[0]*layerFunction(&xyz[0], parLayer)*phiFunction(&xyz[1], parPhi)*timeFunction(&xyz[2], parTime);
0493 
0494               if(observed <= my_data.adc_threshold)
0495               {
0496                 double arg = (expected - my_data.adc_threshold) / (sqrt(2.0) * sigma_ADC);
0497                 double tail_prob = 0.5 * TMath::Erfc(arg);
0498                 nll_val -= log(tail_prob + 1e-12);
0499               }
0500               else if(observed < satThreshold)
0501               {
0502                 double resid = (observed - expected) / sigma_ADC;
0503                 nll_val += 0.5 * (resid * resid + log(2 * TMath::Pi() * sigma_ADC * sigma_ADC));
0504               }
0505               else if(observed >= satThreshold)
0506               {
0507                 double arg = (satThreshold - expected) / (sqrt(2.0) * sigma_ADC);
0508                 double tail_prob = 0.5 * TMath::Erfc(arg);
0509                 nll_val -= log(tail_prob + 1e-12);
0510               }
0511             }
0512           }
0513         }
0514         return nll_val;
0515       };
0516 
0517       fit3D->SetFCN(7, nll, par_init);
0518 
0519       fit3D->Config().ParSettings(0).SetName("amp");
0520       fit3D->Config().ParSettings(0).SetStepSize(10);
0521       fit3D->Config().ParSettings(0).SetLimits(0,5000);
0522       fit3D->Config().ParSettings(1).SetName("mu_layer");
0523       fit3D->Config().ParSettings(1).SetStepSize(0.1);
0524       fit3D->Config().ParSettings(1).SetLimits(usedLayer[0],*usedLayer.rbegin());
0525       fit3D->Config().ParSettings(2).SetName("mu_phi");
0526       fit3D->Config().ParSettings(2).SetStepSize(0.1);
0527       fit3D->Config().ParSettings(2).SetLimits(usedIPhi[0],*usedIPhi.rbegin());
0528       fit3D->Config().ParSettings(3).SetName("sig_phi");
0529       fit3D->Config().ParSettings(3).SetStepSize(0.1);
0530       fit3D->Config().ParSettings(3).SetLimits(0.01,2);
0531       fit3D->Config().ParSettings(4).SetName("mu_t");
0532       fit3D->Config().ParSettings(4).SetStepSize(0.1);
0533       fit3D->Config().ParSettings(4).SetLimits(usedIT[0],*usedIT.rbegin());
0534       fit3D->Config().ParSettings(5).SetName("sig_t");
0535       fit3D->Config().ParSettings(5).SetStepSize(0.1);
0536       fit3D->Config().ParSettings(5).SetLimits(0.01,10);
0537       fit3D->Config().ParSettings(6).SetName("lambda_t");
0538       fit3D->Config().ParSettings(6).SetStepSize(0.01);
0539       fit3D->Config().ParSettings(6).SetLimits(0,5);
0540 
0541  
0542       if(usedLayer.size() == 1)
0543       {
0544         fit3D->Config().ParSettings(1).Fix();
0545       }
0546 
0547       fitSuccess = fit3D->FitFCN();
0548 
0549 
0550       if (my_data.Verbosity > 2)
0551       {
0552         std::cout << "fit success: " << fitSuccess << std::endl;
0553       }
0554     }
0555     pthread_mutex_unlock(&mythreadlock);
0556 
0557 
0558 
0559     if(my_data.doFitting && fitSuccess)
0560     {
0561 
0562       const ROOT::Fit::FitResult& result = fit3D->Result();
0563 
0564 
0565       PHG4TpcGeom *layergeomLow = my_data.geom_container->GetLayerCellGeom((int) floor(result.Parameter(1)));
0566       PHG4TpcGeom *layergeomHigh = my_data.geom_container->GetLayerCellGeom((int) ceil(result.Parameter(1)));
0567 
0568       double RLow = layergeomLow->get_radius();
0569       double RHigh = layergeomHigh->get_radius();
0570 
0571       double phiHigh_RLow = -999.0;
0572       if(ceil(result.Parameter(2)) < layergeomLow->get_phibins())
0573       {
0574         phiHigh_RLow = layergeomLow->get_phi(ceil(result.Parameter(2)), (meanSide < 0 ? 0 : 1));
0575       }
0576       double phiHigh_RHigh = -999.0;
0577       if(ceil(result.Parameter(2)) < layergeomHigh->get_phibins())
0578       {
0579         phiHigh_RHigh = layergeomHigh->get_phi(ceil(result.Parameter(2)), (meanSide < 0 ? 0 : 1));
0580       }
0581 
0582       double phiLow_RLow = layergeomLow->get_phi(floor(result.Parameter(2)), (meanSide < 0 ? 0 : 1));
0583       double phiLow_RHigh = layergeomHigh->get_phi(floor(result.Parameter(2)), (meanSide < 0 ? 0 : 1));
0584 
0585       double meanR = (result.Parameter(1) - floor(result.Parameter(1))) * (RHigh - RLow) + RLow;
0586 
0587       double meanPhi_RLow = ((result.Parameter(2) - floor(result.Parameter(2)))) * (phiHigh_RLow - phiLow_RLow) + phiLow_RLow;
0588       double meanPhi_RHigh = ((result.Parameter(2) - floor(result.Parameter(2)))) * (phiHigh_RHigh - phiLow_RHigh) + phiLow_RHigh;
0589 
0590       double meanPhi = 0.5*(meanPhi_RLow + meanPhi_RHigh);
0591       if(phiHigh_RLow == -999.0 && phiHigh_RHigh != -999.0)
0592       {
0593         meanPhi = meanPhi_RHigh;
0594       }
0595       else if(phiHigh_RLow != -999.0 && phiHigh_RHigh == -999.0)
0596       {
0597         meanPhi = meanPhi_RLow;
0598       }
0599       
0600 
0601       if(phiHigh_RLow == -999.0 && phiHigh_RHigh == -999.0)
0602       {
0603         clus->setAdc(adcSum);
0604         clus->setX(clusX);
0605         clus->setY(clusY);
0606         clus->setZ(clusZ);
0607         clus->setFitMode(false);
0608         clus->setLayer(layerSum / adcSum);
0609         clus->setIPhi(iphiSum / adcSum);
0610         clus->setIT(itSum / adcSum);
0611         clus->setNLayers(usedLayer.size());
0612         clus->setNIPhi(usedIPhi.size());
0613         clus->setNIT(usedIT.size());
0614         clus->setSDLayer(sqrt(sigmaLayer / nHits));
0615         clus->setSDIPhi(sqrt(sigmaIPhi / nHits));
0616         clus->setSDIT(sqrt(sigmaIT / nHits));
0617         clus->setSDWeightedLayer(sqrt(sigmaWeightedLayer / adcSum));
0618         clus->setSDWeightedIPhi(sqrt(sigmaWeightedIPhi / adcSum));
0619         clus->setSDWeightedIT(sqrt(sigmaWeightedIT / adcSum));
0620       }
0621       else
0622       { 
0623         clus->setAdc(adcSum);
0624         clus->setX(meanR*cos(meanPhi));
0625         clus->setY(meanR*sin(meanPhi));
0626         clus->setZ(clusZ);
0627         clus->setFitMode(true);
0628         clus->setLayer(result.Parameter(1));
0629         clus->setIPhi(result.Parameter(2));
0630         clus->setIT(result.Parameter(4));
0631         clus->setNLayers(usedLayer.size());
0632         clus->setNIPhi(usedIPhi.size());
0633         clus->setNIT(usedIT.size());
0634         clus->setSDLayer(sqrt(sigmaLayer / nHits));
0635         clus->setSDIPhi(sqrt(sigmaIPhi / nHits));
0636         clus->setSDIT(sqrt(sigmaIT / nHits));
0637         clus->setSDWeightedLayer(sqrt(sigmaWeightedLayer / adcSum));
0638         clus->setSDWeightedIPhi(result.Parameter(3));
0639         clus->setSDWeightedIT(result.Parameter(5));
0640       }
0641     }
0642     else
0643     {
0644       clus->setAdc(adcSum);
0645       clus->setX(clusX);
0646       clus->setY(clusY);
0647       clus->setZ(clusZ);
0648       clus->setFitMode(false);
0649       clus->setLayer(layerSum / adcSum);
0650       clus->setIPhi(iphiSum / adcSum);
0651       clus->setIT(itSum / adcSum);
0652       clus->setNLayers(usedLayer.size());
0653       clus->setNIPhi(usedIPhi.size());
0654       clus->setNIT(usedIT.size());
0655       clus->setSDLayer(sqrt(sigmaLayer / nHits));
0656       clus->setSDIPhi(sqrt(sigmaIPhi / nHits));
0657       clus->setSDIT(sqrt(sigmaIT / nHits));
0658       clus->setSDWeightedLayer(sqrt(sigmaWeightedLayer / adcSum));
0659       clus->setSDWeightedIPhi(sqrt(sigmaWeightedIPhi / adcSum));
0660       clus->setSDWeightedIT(sqrt(sigmaWeightedIT / adcSum));
0661     }
0662 
0663     const auto ckey = TrkrDefs::genClusKey(maxKey, my_data.cluster_vector.size());
0664     my_data.cluster_vector.push_back(clus);
0665     my_data.cluster_key_vector.push_back(ckey);
0666     
0667     if(fit3D)
0668     {
0669       delete fit3D;
0670     }
0671 
0672     if(my_data.hitHist)
0673     {
0674       delete my_data.hitHist;
0675       my_data.hitHist = nullptr;
0676     }
0677     
0678 
0679   }
0680   
0681 
0682   void ProcessModuleData(thread_data *my_data)
0683   {
0684     
0685     if (my_data->Verbosity > 2)
0686     {
0687       pthread_mutex_lock(&mythreadlock);
0688       std::cout << "working on side: " << my_data->side << "   sector: " << my_data->sector << "   module: " << my_data->module << std::endl;
0689       pthread_mutex_unlock(&mythreadlock);
0690     }
0691 
0692     bgi::rtree<hitData, bgi::quadratic<16>> rtree;
0693 
0694     std::multimap<unsigned int, pointKeyLaser> adcMap;
0695 
0696     if (my_data->hitsets.size() == 0)
0697     {
0698       return;
0699     }
0700     for(int i=0; i<(int)my_data->hitsets.size(); i++)
0701     {
0702       auto *hitset = my_data->hitsets[i];
0703       unsigned int layer = my_data->layers[i];
0704       bool side = my_data->side;
0705       unsigned int sector = my_data->sector;
0706 
0707       TrkrDefs::hitsetkey hitsetKey = TpcDefs::genHitSetKey(layer, sector, (int)side);
0708 
0709       TrkrHitSet::ConstRange hitrangei = hitset->getHits();
0710 
0711       for (TrkrHitSet::ConstIterator hitr = hitrangei.first; hitr != hitrangei.second; ++hitr)
0712       {
0713           float_t fadc = hitr->second->getAdc();
0714           unsigned short adc = 0;
0715           if (fadc > my_data->adc_threshold)
0716           {
0717             adc = (unsigned short) fadc;
0718           }
0719           else
0720           {
0721             continue;
0722           }
0723     
0724           int iphi = TpcDefs::getPad(hitr->first);
0725           int it = TpcDefs::getTBin(hitr->first);
0726         
0727           if(fabs(it - my_data->peakTimeBin) > 5)
0728           {
0729             continue;
0730           }
0731 
0732           point coords = point((int) layer, iphi, it);
0733 
0734           std::vector<hitData> testduplicate;
0735           rtree.query(bgi::intersects(box(point(layer - 0.001, iphi - 0.001, it - 0.001),
0736                     point(layer + 0.001, iphi + 0.001, it + 0.001))),
0737             std::back_inserter(testduplicate));
0738           if (!testduplicate.empty())
0739           {
0740             testduplicate.clear();
0741             continue;
0742           }
0743 
0744           TrkrDefs::hitkey hitKey = TpcDefs::genHitKey(iphi, it);
0745     
0746           auto spechitkey = std::make_pair(hitKey, hitsetKey);
0747           pointKeyLaser coordsKey = std::make_pair(coords, spechitkey);
0748           adcMap.insert(std::make_pair(adc, coordsKey));
0749         auto adckey = std::make_pair(adc, spechitkey);
0750           rtree.insert(std::make_pair(point(1.0*layer, 1.0*iphi, 1.0*it), adckey));
0751       }
0752     }
0753     //finished filling rtree
0754 
0755     while (adcMap.size() > 0)
0756     {
0757       auto iterKey = adcMap.rbegin();
0758       if(iterKey == adcMap.rend())
0759       {
0760           break;
0761       }
0762       
0763 
0764       auto coords = iterKey->second.first;
0765       int layer = coords.get<0>();
0766       int iphi = coords.get<1>();
0767       int it = coords.get<2>();
0768    
0769       if (my_data->Verbosity > 2)
0770       {
0771         pthread_mutex_lock(&mythreadlock);
0772         std::cout << "working on cluster " << my_data->cluster_vector.size() << "   side: " << my_data->side << "   sector: " << my_data->sector << "   module: " << (layer<23 ? 1 : (layer<39 ? 2 : 3) ) << std::endl;
0773         pthread_mutex_unlock(&mythreadlock);
0774 
0775       }
0776 
0777       std::vector<hitData> clusHits;
0778 
0779       rtree.query(bgi::intersects(box(point(layer - my_data->layerMin, iphi - 6, it - 5), point(layer + my_data->layerMax, iphi + 6, it + 5))), std::back_inserter(clusHits));
0780 
0781       calc_cluster_parameter(clusHits, *my_data, iterKey->second.second);
0782 
0783       remove_hits(clusHits, rtree, adcMap);
0784 
0785     }
0786   }
0787 
0788   void *ProcessModule(void *threadarg)
0789   {
0790     auto my_data = static_cast<thread_data *>(threadarg);
0791     ProcessModuleData(my_data);
0792     pthread_exit(nullptr);
0793   }
0794 } //namespace
0795 
0796 LaserClusterizer::LaserClusterizer(const std::string &name)
0797   : SubsysReco(name)
0798 {
0799 }
0800 
0801 int LaserClusterizer::InitRun(PHCompositeNode *topNode)
0802 {
0803   TH1::AddDirectory(kFALSE);
0804 
0805   PHNodeIterator iter(topNode);
0806 
0807   // Looking for the DST node
0808   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0809   if (!dstNode)
0810   {
0811     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0812     return Fun4AllReturnCodes::ABORTRUN;
0813   }
0814 
0815   // Create the Cluster node if required
0816   std::string laserClusterNodeName = "LASER_CLUSTER";
0817   if (m_lamination)
0818   {
0819     laserClusterNodeName = "LAMINATION_CLUSTER";
0820   }
0821   auto laserclusters = findNode::getClass<LaserClusterContainer>(dstNode, laserClusterNodeName);
0822   if (!laserclusters)
0823   {
0824     PHNodeIterator dstiter(dstNode);
0825     PHCompositeNode *DetNode =
0826         dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0827     if (!DetNode)
0828     {
0829       DetNode = new PHCompositeNode("TRKR");
0830       dstNode->addNode(DetNode);
0831     }
0832 
0833     laserclusters = new LaserClusterContainerv1;
0834     PHIODataNode<PHObject> *LaserClusterContainerNode =
0835         new PHIODataNode<PHObject>(laserclusters, laserClusterNodeName, "PHObject");
0836     DetNode->addNode(LaserClusterContainerNode);
0837   }
0838   
0839   m_geom_container =
0840       findNode::getClass<PHG4TpcGeomContainer>(topNode, "TPCGEOMCONTAINER");
0841   if (!m_geom_container)
0842   {
0843     std::cout << PHWHERE << "ERROR: Can't find node TPCGEOMCONTAINER" << std::endl;
0844     return Fun4AllReturnCodes::ABORTRUN;
0845   }
0846   // get the first layer to get the clock freq
0847   AdcClockPeriod = m_geom_container->GetFirstLayerCellGeom()->get_zstep();
0848   m_tdriftmax = AdcClockPeriod * NZBinsSide;
0849 
0850   return Fun4AllReturnCodes::EVENT_OK;
0851 }
0852 
0853 int LaserClusterizer::process_event(PHCompositeNode *topNode)
0854 {
0855   eventHeader = findNode::getClass<EventHeader>(topNode, "EventHeader");
0856   if (!eventHeader)
0857   {
0858     std::cout << PHWHERE << " EventHeader Node missing, doing nothing." << std::endl;
0859     return Fun4AllReturnCodes::ABORTRUN;
0860   }
0861 
0862   m_event = eventHeader->get_EvtSequence();
0863 
0864   if (Verbosity() > 1)
0865   {
0866     std::cout << "LaserClusterizer::process_event working on event " << m_event << std::endl;
0867   }
0868 
0869   m_laserEventInfo = findNode::getClass<LaserEventInfo>(topNode, "LaserEventInfo");
0870   if (!m_laserEventInfo)
0871   {
0872     std::cout << PHWHERE << "ERROR: Can't find node LaserEventInfo" << std::endl;
0873     return Fun4AllReturnCodes::ABORTRUN;
0874   }
0875 
0876   if((eventHeader->get_RunNumber() > 66153 && !m_laserEventInfo->isGl1LaserEvent()) || (eventHeader->get_RunNumber() <= 66153 && !m_laserEventInfo->isLaserEvent()))
0877   {
0878     return Fun4AllReturnCodes::EVENT_OK;
0879   }
0880 
0881   if (Verbosity() > 1)
0882   {
0883     std::cout << "LaserClusterizer::process_event laser event found" << std::endl;
0884   }
0885 
0886   PHNodeIterator iter(topNode);
0887   PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0888   if (!dstNode)
0889   {
0890     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0891     return Fun4AllReturnCodes::ABORTRUN;
0892   }
0893   // get node containing the digitized hits
0894   m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0895   if (!m_hits)
0896   {
0897     std::cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << std::endl;
0898     return Fun4AllReturnCodes::ABORTRUN;
0899   }
0900   
0901   // get node for clusters
0902   std::string laserClusterNodeName = "LASER_CLUSTER";
0903   if (m_lamination)
0904   {
0905     laserClusterNodeName = "LAMINATION_CLUSTER";
0906   }
0907   m_clusterlist = findNode::getClass<LaserClusterContainer>(topNode, laserClusterNodeName);
0908   if (!m_clusterlist)
0909   {
0910     std::cout << PHWHERE << " ERROR: Can't find " << laserClusterNodeName << "." << std::endl;
0911     return Fun4AllReturnCodes::ABORTRUN;
0912   }
0913 
0914   m_tGeometry = findNode::getClass<ActsGeometry>(topNode,
0915                                                  "ActsGeometry");
0916   if (!m_tGeometry)
0917   {
0918     std::cout << PHWHERE
0919               << "ActsGeometry not found on node tree. Exiting"
0920               << std::endl;
0921     return Fun4AllReturnCodes::ABORTRUN;
0922   }
0923 
0924   TrkrHitSetContainer::ConstRange hitsetrange = m_hits->getHitSets(TrkrDefs::TrkrId::tpcId);;
0925   
0926   struct thread_pair_t
0927   {
0928     pthread_t thread{};
0929     thread_data data;
0930   };
0931 
0932   std::vector<thread_pair_t> threads;
0933   threads.reserve(72);
0934 
0935   pthread_attr_t attr;
0936   pthread_attr_init(&attr);
0937   pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
0938 
0939   if (pthread_mutex_init(&mythreadlock, nullptr) != 0)
0940   {
0941     std::cout << std::endl << " mutex init failed" << std::endl;
0942     return 1;
0943   }
0944   
0945   for (unsigned int sec=0; sec<12; sec++)
0946   {
0947     for (int s=0; s<2; s++)
0948     {
0949       for (unsigned int mod=0; mod<3; mod++)
0950       {
0951     
0952         if(Verbosity() > 2)
0953         {
0954           std::cout << "making thread for side: " << s << "   sector: " << sec << "   module: " << mod << std::endl;
0955         }
0956 
0957           thread_pair_t &thread_pair = threads.emplace_back();
0958 
0959           std::vector<TrkrHitSet *> hitsets;
0960           std::vector<unsigned int> layers;
0961 
0962           std::vector<LaserCluster *> cluster_vector;
0963           std::vector<TrkrDefs::cluskey> cluster_key_vector;
0964 
0965           for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
0966               hitsetitr != hitsetrange.second;
0967               ++hitsetitr)
0968           {
0969             unsigned int layer = TrkrDefs::getLayer(hitsetitr->first);
0970             int side = TpcDefs::getSide(hitsetitr->first);
0971             unsigned int sector = TpcDefs::getSectorId(hitsetitr->first);
0972             if (sector != sec || side != s)
0973             {
0974               continue;
0975             }
0976             if ((mod==0 && (layer<7 || layer>22)) || (mod==1 && (layer<=22 || layer>38) ) || (mod==2 && (layer<=38 || layer>54)))
0977             {
0978               continue;
0979             }
0980 
0981             TrkrHitSet *hitset = hitsetitr->second;
0982       
0983             hitsets.push_back(hitset);
0984             layers.push_back(layer);
0985       
0986         }
0987 
0988           thread_pair.data.geom_container = m_geom_container;
0989           thread_pair.data.tGeometry = m_tGeometry;
0990           thread_pair.data.hitsets = hitsets;
0991           thread_pair.data.layers = layers;
0992           thread_pair.data.side = (bool)s;
0993           thread_pair.data.sector = sec;
0994         thread_pair.data.module = mod;
0995           thread_pair.data.cluster_vector = cluster_vector;
0996           thread_pair.data.cluster_key_vector = cluster_key_vector;
0997           thread_pair.data.adc_threshold = m_adc_threshold;
0998           thread_pair.data.peakTimeBin = m_laserEventInfo->getPeakSample(s);
0999           thread_pair.data.layerMin = 3;
1000           thread_pair.data.layerMax = 3;
1001           thread_pair.data.tdriftmax = m_tdriftmax;
1002         thread_pair.data.eventNum = m_event;
1003         thread_pair.data.Verbosity = Verbosity();
1004         thread_pair.data.hitHist = nullptr;
1005         thread_pair.data.doFitting = m_do_fitting;
1006 
1007           int rc;
1008           rc = pthread_create(&thread_pair.thread, &attr, ProcessModule, (void *) &thread_pair.data);
1009 
1010           if (rc)
1011           {
1012             std::cout << "Error:unable to create thread," << rc << std::endl;
1013           }
1014 
1015           if (m_do_sequential)
1016           {
1017             //wait for termination of thread
1018             int rc2 = pthread_join(thread_pair.thread, nullptr);
1019             if (rc2)
1020             {
1021               std::cout << "Error:unable to join," << rc2 << std::endl;
1022             }
1023 
1024           //add clusters from thread to laserClusterContainer
1025             const auto &data(thread_pair.data);
1026             for(int index = 0; index < (int) data.cluster_vector.size(); ++index)
1027             {
1028               auto cluster = data.cluster_vector[index];
1029               const auto ckey = data.cluster_key_vector[index];
1030         
1031             m_clusterlist->addClusterSpecifyKey(ckey, cluster);
1032             }
1033           }
1034       }
1035     }
1036   }
1037   
1038   pthread_attr_destroy(&attr);
1039 
1040   if (!m_do_sequential)
1041   {
1042     for (const auto & thread_pair : threads)
1043     {
1044       int rc2 = pthread_join(thread_pair.thread, nullptr);
1045       if (rc2)
1046       {
1047           std::cout << "Error:unable to join," << rc2 << std::endl;
1048       }
1049       
1050       //const auto &data(thread_pair.data);
1051       
1052       for(int index = 0; index < (int) thread_pair.data.cluster_vector.size(); ++index)
1053       {
1054           auto cluster = thread_pair.data.cluster_vector[index];
1055           const auto ckey = thread_pair.data.cluster_key_vector[index];
1056     
1057           m_clusterlist->addClusterSpecifyKey(ckey, cluster);
1058       }
1059     }
1060   }
1061 
1062   threads.clear();
1063 
1064   if (Verbosity() > 1)
1065   {
1066     std::cout << "LaserClusterizer::process_event " << m_clusterlist->size() << " clusters found" << std::endl;
1067   }
1068 
1069   return Fun4AllReturnCodes::EVENT_OK;
1070 
1071 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!