packages/tpccalib/TpcLaminationFitting.cc

0001 #include "TpcLaminationFitting.h"
0002
0003 #include <trackbase/CMFlashDifferenceContainerv1.h>
0004 #include <trackbase/CMFlashDifferencev1.h>
0005 #include <trackbase/LaserClusterContainerv1.h>
0006 #include <trackbase/LaserClusterv1.h>
0007 #include <trackbase/TpcDefs.h>
0008
0009 #include <ffaobjects/EventHeader.h>
0010
0011 #include <fun4all/Fun4AllReturnCodes.h>
0012
0013 #include <phool/PHCompositeNode.h>
0014 #include <phool/getClass.h>
0015 #include <phool/phool.h>
0016
0017 #include <odbc++/connection.h>
0018 #include <odbc++/drivermanager.h>
0019 #include <odbc++/resultset.h>
0020 #include <odbc++/statement.h>
0021 #include <odbc++/types.h>
0022
0023 #include <TCanvas.h>
0024 #include <TF1.h>
0025 #include <TFile.h>
0026 #include <TGraph.h>
0027 #include <TH2.h>
0028 #include <TH3.h>
0029 #include <TLine.h>
0030 #include <TPaveText.h>
0031 #include <TString.h>
0032 #include <TStyle.h>
0033 #include <TTree.h>
0034 #include <TVector3.h>
0035
0036 #include <boost/format.hpp>
0037
0038 #include <cmath>
0039 #include <iomanip>
0040 #include <set>
0041 #include <string>
0042 #include <vector>
0043
0044 //________________________________
0045 TpcLaminationFitting::TpcLaminationFitting(const std::string &name)
0046   : SubsysReco(name)
0047 {
0048 }
0049
0050 //___________________________________________________________
0051 void TpcLaminationFitting::set_grid_dimensions(int phibins, int rbins)
0052 {
0053   m_phibins = phibins;
0054   m_rbins = rbins;
0055 }
0056
0057 //____________________________________________
0058 int TpcLaminationFitting::InitRun(PHCompositeNode *topNode)
0059 {
0060   for (int s = 0; s < 2; s++)
0061   {
0062     for (int l = 0; l < 18; l++)
0063     {
0064       double shift = l * M_PI / 9;
0065       if (s == 0)
0066       {
0067         shift += M_PI / 18;
0068       }
0069
0070       m_hLamination[l][s] = new TH2D((boost::format("hLamination%d_%s") %l %(s == 1 ? "North" : "South")).str().c_str(), (boost::format("Lamination %d %s, #phi_{expected}=%.2f;R [cm];#phi") %l %(s == 1 ? "North" : "South") %shift).str().c_str(), 200, 30, 80, 200, shift - 0.2, shift + 0.2);
0071       //m_fLamination[l][s] = new TF1((boost::format("fLamination%d_%s") %l %(s == 1 ? "North" : "South")).str().c_str(), "[0]+[1]*exp(-[2]*x)", 30, 80);
0072       //m_fLamination[l][s] = new TF1((boost::format("fLamination%d_%s") %l %(s == 1 ? "North" : "South")).str().c_str(), "[0]*(1+exp(-[2]*(x-[1])))", 30, 80);
0073       m_fLamination[l][s] = new TF1((boost::format("fLamination%d_%s") %l %(s == 1 ? "North" : "South")).str().c_str(), "[3]+[0]*(1-exp(-[2]*(x-[1])))", 30, 80);
0074       //m_fLamination[l][s]->SetParameters(-0.022 + shift, log(3.0/(-0.22 + shift)), 0.12);
0075       m_fLamination[l][s]->SetParameters(-0.011, 30, 0.16, 0.0);
0076       m_fLamination[l][s]->SetParLimits(0, -0.22, 0.0);
0077       m_fLamination[l][s]->SetParLimits(1, 0, 80);
0078       m_fLamination[l][s]->SetParLimits(2, 0.0, 3);
0079       m_fLamination[l][s]->FixParameter(3, shift);
0080       m_laminationCenter[l][s] = shift;
0081     }
0082   }
0083
0084   /*
0085   //Make map for run and ZDC rate for pp mode
0086   m_run_ZDC_map_pp.insert(std::pair<int, float>(49709, 555.0));
0087   m_run_ZDC_map_pp.insert(std::pair<int, float>(52077, 0.0));
0088   m_run_ZDC_map_pp.insert(std::pair<int, float>(52078, 0.0));
0089   m_run_ZDC_map_pp.insert(std::pair<int, float>(53534, 3013.5));
0090   m_run_ZDC_map_pp.insert(std::pair<int, float>(53630, 6849.3));
0091   m_run_ZDC_map_pp.insert(std::pair<int, float>(53631, 5577.8));
0092   m_run_ZDC_map_pp.insert(std::pair<int, float>(53632, 5151.2));
0093   m_run_ZDC_map_pp.insert(std::pair<int, float>(53652, 4600.0));
0094   m_run_ZDC_map_pp.insert(std::pair<int, float>(53687, 3967.2));
0095   m_run_ZDC_map_pp.insert(std::pair<int, float>(53716, 3070.1));
0096   m_run_ZDC_map_pp.insert(std::pair<int, float>(53738, 4510.7));
0097   m_run_ZDC_map_pp.insert(std::pair<int, float>(53739, 4165.0));
0098   m_run_ZDC_map_pp.insert(std::pair<int, float>(53741, 3738.1));
0099   m_run_ZDC_map_pp.insert(std::pair<int, float>(53742, 3721.4));
0100   m_run_ZDC_map_pp.insert(std::pair<int, float>(53743, 3693.4));
0101   m_run_ZDC_map_pp.insert(std::pair<int, float>(53744, 3581.9));
0102   m_run_ZDC_map_pp.insert(std::pair<int, float>(53756, 4471.4));
0103   m_run_ZDC_map_pp.insert(std::pair<int, float>(53783, 4825.7));
0104   m_run_ZDC_map_pp.insert(std::pair<int, float>(53871, 6871.5));
0105   m_run_ZDC_map_pp.insert(std::pair<int, float>(53876, 5082.3));
0106   m_run_ZDC_map_pp.insert(std::pair<int, float>(53877, 4758.5));
0107   m_run_ZDC_map_pp.insert(std::pair<int, float>(53879, 4315.0));
0108
0109   //beam off go into pp
0110   m_run_ZDC_map_pp.insert(std::pair<int, float>(53098, 0.0));
0111   m_run_ZDC_map_pp.insert(std::pair<int, float>(53271, 0.0));
0112
0113   m_run_ZDC_map_auau.insert(std::pair<int, float>(54966, 12400.));
0114   m_run_ZDC_map_auau.insert(std::pair<int, float>(54967, 11600.));
0115   m_run_ZDC_map_auau.insert(std::pair<int, float>(54968, 10500.));
0116   m_run_ZDC_map_auau.insert(std::pair<int, float>(54969, 9680.));
0117   */
0118
0119   int ret = GetNodes(topNode);
0120   return ret;
0121 }
0122
0123 //______________________________________
0124 int TpcLaminationFitting::GetNodes(PHCompositeNode *topNode)
0125 {
0126   //m_correctedCMcluster_map = findNode::getClass<LaserClusterContainer>(topNode, "LAMINATION_CLUSTER");
0127   m_correctedCMcluster_map = findNode::getClass<LaserClusterContainer>(topNode, "LASER_CLUSTER");
0128   if (!m_correctedCMcluster_map)
0129   {
0130     std::cout << PHWHERE << "CORRECTED_CM_CLUSTER Node missing, abort." << std::endl;
0131     return Fun4AllReturnCodes::ABORTRUN;
0132   }
0133
0134   m_dcc_in_module_edge = findNode::getClass<TpcDistortionCorrectionContainer>(topNode, "TpcDistortionCorrectionContainerModuleEdge");
0135   if (m_dcc_in_module_edge)
0136   {
0137     std::cout << "TpcLaminationFitting::GetNodes - found TPC distortion correction container module edge" << std::endl;
0138   }
0139
0140   m_dcc_in_static = findNode::getClass<TpcDistortionCorrectionContainer>(topNode, "TpcDistortionCorrectionContainerStatic");
0141   if (m_dcc_in_static)
0142   {
0143     std::cout << "TpcLaminationFitting::GetNodes - found TPC distortion correction container static" << std::endl;
0144   }
0145
0146   PHNodeIterator iter(topNode);
0147
0148   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0149   if (!dstNode)
0150   {
0151     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0152     return Fun4AllReturnCodes::ABORTRUN;
0153   }
0154
0155   auto flashDiffContainer = findNode::getClass<CMFlashDifferenceContainerv1>(topNode, "CM_FLASH_DIFFERENCES");
0156   if (!flashDiffContainer)
0157   {
0158     PHNodeIterator dstIter(dstNode);
0159     PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstIter.findFirst("PHCompositeNode", "TRKR"));
0160     if (!DetNode)
0161     {
0162       DetNode = new PHCompositeNode("TRKR");
0163       dstNode->addNode(DetNode);
0164     }
0165
0166     flashDiffContainer = new CMFlashDifferenceContainerv1;
0167     PHIODataNode<PHObject> *CMFlashDifferenceNode = new PHIODataNode<PHObject>(flashDiffContainer, "CM_FLASH_DIFFERENCES", "PHObject");
0168     DetNode->addNode(CMFlashDifferenceNode);
0169   }
0170
0171   const std::string dcc_out_node_name = "TpcDistortionCorrectionContainerAverage";
0172   m_dcc_out = findNode::getClass<TpcDistortionCorrectionContainer>(topNode, dcc_out_node_name);
0173   if (!m_dcc_out)
0174   {
0175     auto runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0176     if (!runNode)
0177     {
0178       std::cout << "TpcLaminationFitting::InitRun = RUN Node missing, quitting" << std::endl;
0179       return Fun4AllReturnCodes::ABORTRUN;
0180     }
0181
0182     std::cout << "TpcLaminationFitting::GetNodes - creating TpcDistortionCorrectionContainer in node " << dcc_out_node_name << std::endl;
0183     m_dcc_out = new TpcDistortionCorrectionContainer;
0184     auto node = new PHDataNode<TpcDistortionCorrectionContainer>(m_dcc_out, dcc_out_node_name);
0185     runNode->addNode(node);
0186   }
0187
0188   const float phiMin = m_phiMin - (m_phiMax - m_phiMin) / m_phibins;
0189   const float phiMax = m_phiMax + (m_phiMax - m_phiMin) / m_phibins;
0190
0191   const float rMin = m_rMin - (m_rMax - m_rMin) / m_rbins;
0192   const float rMax = m_rMax + (m_rMax - m_rMin) / m_rbins;
0193
0194   const std::array<const std::string, 2> extension = {{"_negz", "_posz"}};
0195
0196   m_dcc_out->m_dimensions = 2;
0197
0198   for (int i = 0; i < 2; ++i)
0199   {
0200     delete m_dcc_out->m_hDPint[i];
0201     m_dcc_out->m_hDPint[i] = new TH2F((boost::format("hIntDistortionP%s") % extension[i]).str().c_str(), (boost::format("hIntDistortionP%s") % extension[i]).str().c_str(), m_phibins + 2, phiMin, phiMax, m_rbins + 2, rMin, rMax);
0202     delete m_dcc_out->m_hDRint[i];
0203     m_dcc_out->m_hDRint[i] = new TH2F((boost::format("hIntDistortionR%s") % extension[i]).str().c_str(), (boost::format("hIntDistortionR%s") % extension[i]).str().c_str(), m_phibins + 2, phiMin, phiMax, m_rbins + 2, rMin, rMax);
0204     delete m_dcc_out->m_hDZint[i];
0205     m_dcc_out->m_hDZint[i] = new TH2F((boost::format("hIntDistortionZ%s") % extension[i]).str().c_str(), (boost::format("hIntDistortionZ%s") % extension[i]).str().c_str(), m_phibins + 2, phiMin, phiMax, m_rbins + 2, rMin, rMax);
0206     delete m_dcc_out->m_hentries[i];
0207     m_dcc_out->m_hentries[i] = new TH2I((boost::format("hEntries%s") % extension[i]).str().c_str(), (boost::format("hEntries%s") % extension[i]).str().c_str(), m_phibins + 2, phiMin, phiMax, m_rbins + 2, rMin, rMax);
0208   }
0209
0210   m_laminationTree = new TTree("laminationTree","laminationTree");
0211   m_laminationTree->Branch("side",&m_side);
0212   m_laminationTree->Branch("lamIndex",&m_lamIndex);
0213   m_laminationTree->Branch("lamPhi",&m_lamPhi);
0214   m_laminationTree->Branch("goodFit",&m_goodFit);
0215   m_laminationTree->Branch("A",&m_A);
0216   m_laminationTree->Branch("B",&m_B);
0217   m_laminationTree->Branch("C",&m_C);
0218   m_laminationTree->Branch("A_err",&m_A_err);
0219   m_laminationTree->Branch("B_err",&m_B_err);
0220   m_laminationTree->Branch("C_err",&m_C_err);
0221   m_laminationTree->Branch("distanceToFit",&m_dist);
0222   m_laminationTree->Branch("nBinsFit",&m_nBins);
0223
0224
0225   return Fun4AllReturnCodes::EVENT_OK;
0226 }
0227
0228 //__________________________________________
0229 int TpcLaminationFitting::process_event(PHCompositeNode *topNode)
0230 {
0231   eventHeader = findNode::getClass<EventHeader>(topNode, "EventHeader");
0232   if (!eventHeader)
0233   {
0234     std::cout << PHWHERE << " EventHeader Node missing, abort" << std::endl;
0235     return Fun4AllReturnCodes::ABORTRUN;
0236   }
0237
0238   m_runnumber = eventHeader->get_RunNumber();
0239
0240   if (m_useHeader && eventHeader->get_EvtSequence() == 0)
0241   {
0242     m_useHeader = false;
0243   }
0244
0245   if (m_useHeader)
0246   {
0247     m_event_index = eventHeader->get_EvtSequence();
0248   }
0249
0250   if (!m_correctedCMcluster_map || m_correctedCMcluster_map->size() < 1000)
0251   {
0252     if (!m_useHeader)
0253     {
0254       m_event_index++;
0255     }
0256     return Fun4AllReturnCodes::EVENT_OK;
0257   }
0258
0259   m_nClusters += m_correctedCMcluster_map->size();
0260   m_nEvents++;
0261
0262   auto clusrange = m_correctedCMcluster_map->getClusters();
0263   for (auto cmitr = clusrange.first; cmitr != clusrange.second; ++cmitr)
0264   {
0265     const auto &[cmkey, cmclus_orig] = *cmitr;
0266     LaserCluster *cmclus = cmclus_orig;
0267     // const unsigned int adc = cmclus->getAdc();
0268     bool side = (bool) TpcDefs::getSide(cmkey);
0269     if (cmclus->getNLayers() <= m_nLayerCut)
0270     {
0271       continue;
0272     }
0273
0274     //Acts::Vector3 pos(cmclus->getX(), cmclus->getY(), cmclus->getZ());
0275     Acts::Vector3 pos(cmclus->getX(), cmclus->getY(), (side ? 1.0 : -1.0));
0276     if (m_dcc_in_module_edge)
0277     {
0278       pos = m_distortionCorrection.get_corrected_position(pos, m_dcc_in_module_edge);
0279     }
0280     if (m_dcc_in_static)
0281     {
0282       pos = m_distortionCorrection.get_corrected_position(pos, m_dcc_in_static);
0283     }
0284
0285     TVector3 tmp_pos(pos[0], pos[1], pos[2]);
0286
0287     for (int l = 0; l < 18; l++)
0288     {
0289       double shift = m_laminationCenter[l][side];
0290
0291       double phi2pi = tmp_pos.Phi();
0292       if (side && phi2pi < -0.2)
0293       {
0294         phi2pi += 2 * M_PI;
0295       }
0296       if (!side && phi2pi < M_PI / 18 - 0.2)
0297       {
0298         phi2pi += 2 * M_PI;
0299       }
0300
0301       if (phi2pi > shift - 0.2 && phi2pi < shift + 0.2)
0302       {
0303         m_hLamination[l][side]->Fill(tmp_pos.Perp(), phi2pi);
0304       }
0305     }
0306   }
0307
0308   return Fun4AllReturnCodes::EVENT_OK;
0309 }
0310
0311 //____________________________________
0312 int TpcLaminationFitting::fitLaminations()
0313 {
0314   //int nBinAvg = 6;
0315   int nBinAvg = 4;
0316   // double contentCut = 0.0;
0317
0318   /*
0319   for(int s=0; s<2; s++)
0320   {
0321           for(int l=0; l<18; l++)
0322           {
0323                   if(m_hLamination[l][s]->GetMaximum()/5.0 > contentCut)
0324                   {
0325                           contentCut = m_hLamination[l][s]->GetMaximum()/5.0;
0326                   }
0327           }
0328   }
0329   */
0330   //double seedScale = (m_nClusters / m_nEvents) / 3718.8030;
0331
0332   //float ZDC = 4500.0;
0333   TF1 *Af[2] = {new TF1("AN","pol1",0,100000), new TF1("AS","pol1",0,100000)};
0334   TF1 *Bf[2] = {new TF1("BN","pol1",0,100000), new TF1("BS","pol1",0,100000)};
0335   double Cseed[2] = {0.16, 0.125};
0336
0337   if(ppMode)
0338   {
0339     /*
0340     auto it = m_run_ZDC_map_pp.find(m_runnumber);
0341     if( it != m_run_ZDC_map_pp.end() )
0342     {
0343       std::cout << "pp runnumber " << m_runnumber << " found. It has ZDC NS rate of " << it->second << std::endl;
0344       ZDC = it->second;
0345     }
0346     else
0347     {
0348       std::cout << "pp runnumber " << m_runnumber << " not found. Using default value of " << ZDC << std::endl;
0349     }
0350     */
0351
0352     std::cout << "in ppMode with runnumber " << m_runnumber << " which has a ZDC coincidence rate of " << m_ZDC_coincidence << std::endl;
0353
0354     Af[0]->SetParameters(-0.007999,-1.783e-6);
0355     Af[1]->SetParameters(-0.003288,-2.297e-6);
0356
0357     Bf[0]->SetParameters(31.55,0.0006141);
0358     Bf[1]->SetParameters(34.7,0.0005226);
0359   }
0360   else
0361   {
0362     /*
0363     ZDC = 10000.0;
0364
0365     auto it = m_run_ZDC_map_auau.find(m_runnumber);
0366     if( it != m_run_ZDC_map_auau.end() )
0367     {
0368       std::cout << "AuAu runnumber " << m_runnumber << " found. It has ZDC NS rate of " << it->second << std::endl;
0369       ZDC = it->second;
0370     }
0371     else
0372     {
0373       std::cout << "AuAu runnumber " << m_runnumber << " not found. Using default value of " << ZDC << std::endl;
0374     }
0375     */
0376
0377     std::cout << "in AuAuMode with runnumber " << m_runnumber << " which has a ZDC coincidence rate of " << m_ZDC_coincidence << std::endl;
0378
0379
0380     Af[0]->SetParameters(-0.003836,-1.025e-6);
0381     Af[1]->SetParameters(-0.003283,-8.176e-7);
0382
0383     Bf[0]->SetParameters(32.96,0.0002997);
0384     Bf[1]->SetParameters(31.19,0.0005622);
0385
0386     Cseed[0] = 0.125;
0387     Cseed[1] = 0.122;
0388   }
0389
0390
0391
0392   for (int s = 0; s < 2; s++)
0393   {
0394     for (int l = 0; l < 18; l++)
0395     {
0396       if (m_hLamination[l][s]->GetEntries() < 100)
0397       {
0398         m_laminationGoodFit[l][s] = false;
0399         continue;
0400       }
0401
0402       // double histMax = m_hLamination[l][s]->GetMaximum();
0403
0404       TGraph *gr = new TGraph();
0405       TGraph *proj = new TGraph();
0406
0407       //m_fLamination[l][s]->SetParameters(-0.022 + m_laminationCenter[l][s], 4.595 * seedScale, 0.138);
0408       //m_fLamination[l][s]->SetParameters(-0.022 + m_laminationCenter[l][s], log(4.595 * seedScale/(-0.022 + m_laminationCenter[l][s])), 0.138);
0409       //m_fLamination[l][s]->SetParameters(-0.011 + m_laminationCenter[l][s], 0.025, 0.16);
0410       //m_fLamination[l][s]->SetParameters(-0.011, 30, 0.16, m_laminationCenter[l][s]);
0411       m_fLamination[l][s]->SetParameters(Af[s]->Eval(m_ZDC_coincidence), Bf[s]->Eval(m_ZDC_coincidence), Cseed[s], m_laminationCenter[l][s]);
0412       m_fLamination[l][s]->FixParameter(3, m_laminationCenter[l][s]);
0413
0414       TF1 *fitSeed = (TF1 *) m_fLamination[l][s]->Clone();
0415       fitSeed->SetName((boost::format("fitSeed%d_%s") %l %(s == 1 ? "North" : "South")).str().c_str());
0416
0417       TF1 *localFit = (TF1 *) m_fLamination[l][s]->Clone();
0418       localFit->SetName((boost::format("localFit%d_%s") %l %(s == 1 ? "North" : "South")).str().c_str());
0419
0420       for (int i = m_hLamination[l][s]->GetNbinsX(); i >= 1; i--)
0421       {
0422         double R = m_hLamination[l][s]->GetXaxis()->GetBinCenter(i);
0423         double sum = 0.0;
0424         double c = 0.0;
0425         double noFitPred = fitSeed->Eval(R);
0426         int noFitPredBin = m_hLamination[l][s]->GetYaxis()->FindBin(noFitPred);
0427
0428         double pred2 = -999;
0429         int pred2Bin = -999;
0430
0431         if (gr->GetN() >= 5)
0432         {
0433           gr->Fit(localFit);
0434           pred2 = localFit->Eval(R);
0435           pred2Bin = m_hLamination[l][s]->GetYaxis()->FindBin(pred2);
0436         }
0437
0438         if (pred2 == -999)
0439         {
0440           proj->AddPoint(R, noFitPred);
0441         }
0442         else
0443         {
0444           proj->AddPoint(R, pred2);
0445         }
0446
0447         for (int j = 1; j <= m_hLamination[l][s]->GetNbinsY(); j++)
0448         {
0449           double phi = m_hLamination[l][s]->GetYaxis()->GetBinCenter(j);
0450
0451           if (fabs(phi - m_laminationCenter[l][s]) > 0.05)
0452           {
0453             continue;
0454           }
0455           if (pred2 == -999 && std::abs(j - noFitPredBin) > nBinAvg)
0456           {
0457             continue;
0458           }
0459           if (pred2 != -999 && std::abs(j - noFitPredBin) > nBinAvg && std::abs(j - pred2Bin) > nBinAvg)
0460           {
0461             continue;
0462           }
0463
0464           double content = m_hLamination[l][s]->GetBinContent(i, j);
0465           // if(content < histMax/20.0 || content > histMax/10.0)
0466           if (content <= 0)
0467           {
0468             continue;
0469           }
0470
0471           sum += content * phi;
0472           c += content;
0473         }
0474
0475         if (c > 0)
0476         {
0477           gr->AddPoint(R, sum / c);
0478         }
0479       }
0480
0481       gr->Fit(m_fLamination[l][s]);
0482
0483       // m_hLamination[l][s]->GetXaxis()->SetRangeUser(contentCut, 2*m_hLamination[l][s]->GetMaximum());
0484       // m_hLamination[l][s]->GetXaxis()->SetRangeUser(histMax/20.0,histMax/10.0);
0485
0486       double distToFunc = 0.0;
0487       int nBinsUsed = 0;
0488
0489       for (int i = 1; i <= m_hLamination[l][s]->GetNbinsX(); i++)
0490       {
0491         double R = m_hLamination[l][s]->GetXaxis()->GetBinCenter(i);
0492         double funcVal = m_fLamination[l][s]->Eval(R);
0493         int funcBin = m_hLamination[l][s]->GetYaxis()->FindBin(funcVal);
0494         for (int j = 0; j <= nBinAvg; j++)
0495         {
0496           double contentp = -999;
0497           double contentm = -999;
0498
0499           if (funcBin + j < m_hLamination[l][s]->GetNbinsY())
0500           {
0501             contentp = m_hLamination[l][s]->GetBinContent(i, funcBin + j);
0502           }
0503           if (funcBin - j >= 1)
0504           {
0505             contentm = m_hLamination[l][s]->GetBinContent(i, funcBin - j);
0506           }
0507           // if(contentp >= contentCut || contentm >= contentCut)
0508           // if(contentp >= histMax/20.0 || contentm >= histMax/20.0)
0509           if (contentp >= 1 || contentm >= 1)
0510           {
0511             distToFunc += j;
0512             nBinsUsed++;
0513             break;
0514           }
0515         }
0516       }
0517
0518       m_distanceToFit[l][s] = distToFunc / nBinsUsed;
0519       m_nBinsFit[l][s] = nBinsUsed;
0520       if (nBinsUsed < 10 || distToFunc / nBinsUsed > 1.0)
0521       {
0522         m_laminationGoodFit[l][s] = false;
0523       }
0524       else
0525       {
0526         m_laminationGoodFit[l][s] = true;
0527       }
0528     }
0529   }
0530
0531   return Fun4AllReturnCodes::EVENT_OK;
0532 }
0533
0534 int TpcLaminationFitting::InterpolatePhiDistortions(TH2 *simPhiDistortion[2])
0535 {
0536   phiDistortionLamination[0] = (TH2 *) simPhiDistortion[0]->Clone();
0537   phiDistortionLamination[0]->Reset();
0538   phiDistortionLamination[0]->SetName("phiDistortionLamination0");
0539   phiDistortionLamination[1] = (TH2 *) simPhiDistortion[1]->Clone();
0540   phiDistortionLamination[1]->Reset();
0541   phiDistortionLamination[1]->SetName("phiDistortionLamination1");
0542
0543   for (int s = 0; s < 2; s++)
0544   {
0545     for (int l = 0; l < 18; l++)
0546     {
0547       if (!m_laminationGoodFit[l][s])
0548       {
0549         continue;
0550       }
0551
0552       for (int i = 1; i <= phiDistortionLamination[s]->GetNbinsY(); i++)
0553       {
0554         double R = phiDistortionLamination[s]->GetYaxis()->GetBinCenter(i);
0555         if (R < 30 || R > 80)
0556         {
0557           continue;
0558         }
0559         double phi = m_fLamination[l][s]->Eval(R);
0560         if (phi < 0.0)
0561         {
0562           phi += 2 * M_PI;
0563         }
0564         if (phi > 2 * M_PI)
0565         {
0566           phi -= 2 * M_PI;
0567         }
0568         int phiBin = phiDistortionLamination[s]->GetXaxis()->FindBin(phi);
0569         //m_fLamination[l][s]->SetParameter(0, m_fLamination[l][s]->GetParameter(0) - m_laminationCenter[l][s]);
0570     m_fLamination[l][s]->SetParameter(3, m_laminationOffset[l][s]);
0571     /*
0572     if(s==0)
0573     {
0574       m_fLamination[l][s]->SetParameter(3, 0.0);
0575     }
0576         else
0577     {
0578       m_fLamination[l][s]->SetParameter(3, 0.0);
0579     }
0580     */
0581     //m_fLamination[l][s]->SetParameter(3, 0.0);
0582         double phiDistortion = R * m_fLamination[l][s]->Integral(phiDistortionLamination[s]->GetYaxis()->GetBinLowEdge(i), phiDistortionLamination[s]->GetYaxis()->GetBinLowEdge(i + 1)) / (phiDistortionLamination[s]->GetYaxis()->GetBinLowEdge(i + 1) - phiDistortionLamination[s]->GetYaxis()->GetBinLowEdge(i));
0583         //m_fLamination[l][s]->SetParameter(0, m_fLamination[l][s]->GetParameter(0) + m_laminationCenter[l][s]);
0584         m_fLamination[l][s]->SetParameter(3, m_laminationCenter[l][s]);
0585         phiDistortionLamination[s]->SetBinContent(phiBin, i, phiDistortion);
0586       }
0587     }
0588   }
0589
0590   for (int s = 0; s < 2; s++)
0591   {
0592     m_dcc_out->m_hDPint[s] = (TH2 *) phiDistortionLamination[s]->Clone();
0593     m_dcc_out->m_hDPint[s]->SetName((boost::format("hIntDistortionP%s") %(s == 0 ? "_negz" : "_posz")).str().c_str());
0594   }
0595
0596   /*
0597   for(int s=0; s<2; s++)
0598   {
0599           for(int i=1; i<=m_dcc_out->m_hDPint[s]->GetNbinsX(); i++)
0600           {
0601                   for(int j=1; j<=m_dcc_out->m_hDPint[s]->GetNbinsY(); j++)
0602                   {
0603                           if(phiDistortionLamination[s]->GetBinContent(i,j) != 0.0)
0604                           {
0605                                   m_dcc_out->m_hDPint[s]->SetBinContent(i,j, phiDistortionLamination[s]->GetBinContent(i,j));
0606                           }
0607                   }
0608           }
0609   }
0610   */
0611
0612   // m_dcc_out->m_hDPint[0] = (TH2*)phiDistortionLamination[0]->Clone();
0613   // m_dcc_out->m_hDPint[1] = (TH2*)phiDistortionLamination[1]->Clone();
0614
0615   for (int s = 0; s < 2; s++)
0616   {
0617     for (int i = 1; i <= m_dcc_out->m_hDPint[s]->GetNbinsY(); i++)
0618     {
0619       double R = m_dcc_out->m_hDPint[s]->GetYaxis()->GetBinCenter(i);
0620       if (R < 30 || R > 80)
0621       {
0622         continue;
0623       }
0624       std::vector<int> laminationPhiBins;
0625       for (int j = 1; j <= m_dcc_out->m_hDPint[s]->GetNbinsX(); j++)
0626       {
0627         if (m_dcc_out->m_hDPint[s]->GetBinContent(j, i) != 0.0)
0628         {
0629           laminationPhiBins.push_back(j);
0630         }
0631       }
0632
0633       if (laminationPhiBins.size() > 1)
0634       {
0635         laminationPhiBins.push_back(laminationPhiBins[0]);
0636       }
0637
0638       for (int lamPair = 0; lamPair < (int) laminationPhiBins.size() - 1; lamPair++)
0639       {
0640         double dist0 = m_dcc_out->m_hDPint[s]->GetBinContent(laminationPhiBins[lamPair], i);
0641         double dist1 = m_dcc_out->m_hDPint[s]->GetBinContent(laminationPhiBins[lamPair + 1], i);
0642
0643         int nEmptyBins = laminationPhiBins[lamPair + 1] - laminationPhiBins[lamPair] - 1;
0644         if (laminationPhiBins[lamPair] > laminationPhiBins[lamPair + 1])
0645         {
0646           double binW = m_dcc_out->m_hDPint[s]->GetXaxis()->GetBinWidth(1);
0647           nEmptyBins = (int) round((m_dcc_out->m_hDPint[s]->GetXaxis()->GetBinCenter(laminationPhiBins[lamPair + 1]) + 2 * M_PI - m_dcc_out->m_hDPint[s]->GetXaxis()->GetBinCenter(laminationPhiBins[lamPair])) / binW) - 1;
0648         }
0649
0650         bool wrap = false;
0651         int wrapBin = -1;
0652         for (int j = 1; j <= nEmptyBins; j++)
0653         // for(int j=+ 1; j<laminationPhiBins[lamPair+1]; j++)
0654         {
0655           if (!wrap && laminationPhiBins[lamPair] + j > m_dcc_out->m_hDPint[s]->GetNbinsX() - 1)
0656           {
0657             wrap = true;
0658             wrapBin = j;
0659           }
0660           double dist = dist0 + ((dist1 - dist0) * (1.0 * j / (nEmptyBins + 1)));
0661           if (wrap)
0662           {
0663             m_dcc_out->m_hDPint[s]->SetBinContent(2 + (j - wrapBin), i, dist);
0664           }
0665           else
0666           {
0667             m_dcc_out->m_hDPint[s]->SetBinContent(laminationPhiBins[lamPair] + j, i, dist);
0668           }
0669         }
0670       }
0671     }
0672   }
0673
0674   return Fun4AllReturnCodes::EVENT_OK;
0675 }
0676
0677 int TpcLaminationFitting::End(PHCompositeNode * /*topNode*/)
0678 {
0679
0680   odbc::Connection *dbConnection = odbc::DriverManager::getConnection("daq", "", "");
0681   std::string sql = "SELECT * FROM gl1_scalers WHERE runnumber = " + std::to_string(m_runnumber) + ";";
0682   odbc::Statement *stmt = dbConnection->createStatement();
0683   odbc::ResultSet *resultSet = stmt->executeQuery(sql);
0684   std::array<std::array<uint64_t, 3>, 64> scalers{};  // initialize to zero
0685   if (!resultSet)
0686   {
0687     std::cerr << "No db found for run number " << m_runnumber << ". Cannot get ZDC rate so aborting run" << std::endl;
0688     delete resultSet;
0689     delete stmt;
0690     delete dbConnection;
0691     return Fun4AllReturnCodes::ABORTRUN;
0692   }
0693
0694   while (resultSet->next())
0695   {
0696     int index = resultSet->getInt("index");
0697     // Iterate over the columns and fill the TriggerRunInfo object
0698     scalers[index][0] = resultSet->getLong("scaled");
0699     scalers[index][1] = resultSet->getLong("live");
0700     scalers[index][2] = resultSet->getLong("raw");
0701   }
0702
0703   delete resultSet;
0704   delete stmt;
0705   delete dbConnection;
0706
0707   m_ZDC_coincidence = (1.0*scalers[3][2]/scalers[0][2])/(106e-9);
0708
0709   std::cout << "Runnumber: " << m_runnumber << "   ppMode: " << ppMode << "   ZDC coindicence rate: " << m_ZDC_coincidence << std::endl;
0710
0711   int fitSuccess = fitLaminations();
0712   if (fitSuccess != Fun4AllReturnCodes::EVENT_OK)
0713   {
0714     std::cout << PHWHERE << " Return code for lamination fitting was " << fitSuccess << " and not successful" << std::endl;
0715     return Fun4AllReturnCodes::ABORTRUN;
0716   }
0717
0718   if(m_QAFileName != "")
0719   {
0720     TCanvas *c1 = new TCanvas();
0721     gStyle->SetPalette(56);
0722     c1->SetLogz();
0723     gStyle->SetOptStat(0);
0724     c1->SaveAs((boost::format("%s[") % m_QAFileName).str().c_str());
0725     for (int s = 0; s < 2; s++)
0726     {
0727       for (int l = 0; l < 18; l++)
0728       {
0729     c1->Clear();
0730     m_hLamination[l][s]->Draw("COLZ");
0731     m_fLamination[l][s]->SetLineColor(kRed);
0732     if (!m_laminationGoodFit[l][s])
0733     {
0734       m_fLamination[l][s]->SetLineColor(kMagenta);
0735     }
0736     m_fLamination[l][s]->Draw("same");
0737
0738     TLine *line = new TLine(30,m_laminationCenter[l][s],80,m_laminationCenter[l][s]);
0739     line->SetLineColor(kBlue);
0740     line->SetLineStyle(2);
0741     line->Draw("same");
0742
0743     TPaveText *pars = new TPaveText(0.6, 0.55, 0.85, 0.85, "NDC");
0744     pars->AddText("#phi = #phi_{ideal} + A#times (1 - e^{-C#times (R - B)})");
0745     pars->AddText((boost::format("A=%.3f#pm %.3f") %m_fLamination[l][s]->GetParameter(0) %m_fLamination[l][s]->GetParError(0)).str().c_str());
0746     pars->AddText((boost::format("#phi_{ideal}=%.3f#pm %.3f") %m_fLamination[l][s]->GetParameter(3) %m_fLamination[l][s]->GetParError(3)).str().c_str());
0747     pars->AddText((boost::format("B=%.3f#pm %.3f") %m_fLamination[l][s]->GetParameter(1) %m_fLamination[l][s]->GetParError(1)).str().c_str());
0748     pars->AddText((boost::format("C=%.3f#pm %.3f") %m_fLamination[l][s]->GetParameter(2) %m_fLamination[l][s]->GetParError(2)).str().c_str());
0749     pars->AddText((boost::format("Distance to line=%.2f") %m_distanceToFit[l][s]).str().c_str());
0750     pars->AddText((boost::format("Number of Bins used=%d") %m_nBinsFit[l][s]).str().c_str());
0751     pars->Draw("same");
0752     c1->SaveAs(m_QAFileName.c_str());
0753       }
0754     }
0755     c1->SaveAs((boost::format("%s]") %m_QAFileName).str().c_str());
0756   }
0757
0758   TFile *simDistortion = new TFile("/cvmfs/sphenix.sdcc.bnl.gov/gcc-12.1.0/release/release_new/new.10/share/calibrations/distortion_maps/average_minus_static_distortion_inverted_10-new.root", "READ");
0759   TH3 *hIntDistortionP_posz = (TH3 *) simDistortion->Get("hIntDistortionP_posz");
0760   hIntDistortionP_posz->GetZaxis()->SetRange(2, 2);
0761   TH2 *simPhiDistortion[2];
0762   simPhiDistortion[1] = (TH2 *) hIntDistortionP_posz->Project3D("yx");
0763   TH3 *hIntDistortionP_negz = (TH3 *) simDistortion->Get("hIntDistortionP_negz");
0764   hIntDistortionP_negz->GetZaxis()->SetRange(hIntDistortionP_negz->GetNbinsZ() - 1, hIntDistortionP_negz->GetNbinsZ() - 1);
0765   simPhiDistortion[0] = (TH2 *) hIntDistortionP_negz->Project3D("yx");
0766
0767   int interpolateSuccess = InterpolatePhiDistortions(simPhiDistortion);
0768   if (interpolateSuccess != Fun4AllReturnCodes::EVENT_OK)
0769   {
0770     std::cout << PHWHERE << " Return code for lamination interpolation was " << interpolateSuccess << " and not successful" << std::endl;
0771     return Fun4AllReturnCodes::ABORTRUN;
0772   }
0773
0774   for (int s = 0; s < 2; s++)
0775   {
0776     scaleFactorMap[s] = (TH2 *) m_dcc_out->m_hDPint[s]->Clone();
0777     scaleFactorMap[s]->SetName((boost::format("scaleFactorMap%d") %s).str().c_str());
0778     scaleFactorMap[s]->Divide(simPhiDistortion[s]);
0779   }
0780
0781   TH3 *hIntDistortionR_posz = (TH3 *) simDistortion->Get("hIntDistortionR_posz");
0782   hIntDistortionR_posz->GetZaxis()->SetRange(2, 2);
0783   TH2 *simRDistortion[2];
0784   simRDistortion[1] = (TH2 *) hIntDistortionR_posz->Project3D("yx");
0785   TH3 *hIntDistortionR_negz = (TH3 *) simDistortion->Get("hIntDistortionR_negz");
0786   hIntDistortionR_negz->GetZaxis()->SetRange(hIntDistortionR_negz->GetNbinsZ() - 1, hIntDistortionR_negz->GetNbinsZ() - 1);
0787   simRDistortion[0] = (TH2 *) hIntDistortionR_negz->Project3D("yx");
0788
0789   for (int s = 0; s < 2; s++)
0790   {
0791     /*
0792     for(int i=1; i<=m_dcc_out->m_hDRint[s]->GetNbinsX(); i++)
0793     {
0794             for(int j=1; j<=m_dcc_out->m_hDRint[s]->GetNbinsY(); j++)
0795             {
0796                     if(simRDistortion[s]->GetBinContent(i,j) != 0.0)
0797                     {
0798                             m_dcc_out->m_hDRint[s]->SetBinContent(i,j, simRDistortion[s]->GetBinContent(i,j));
0799                     }
0800             }
0801     }
0802     */
0803     m_dcc_out->m_hDRint[s] = (TH2 *) simRDistortion[s]->Clone();
0804     m_dcc_out->m_hDRint[s]->SetName((boost::format("hIntDistortionR%s") %(s == 0 ? "_negz" : "_posz")).str().c_str());
0805     m_dcc_out->m_hDRint[s]->Multiply(scaleFactorMap[s]);
0806   }
0807
0808
0809   fill_guarding_bins(m_dcc_out);
0810
0811
0812   for(int s=0; s<2; s++)
0813   {
0814     for(int l=0; l<18; l++)
0815     {
0816       m_side = s;
0817       m_lamIndex = s*18 + l;
0818       m_lamPhi = m_laminationCenter[l][s];
0819       m_goodFit = m_laminationGoodFit[l][s];
0820       m_A = m_fLamination[l][s]->GetParameter(0);
0821       m_B = m_fLamination[l][s]->GetParameter(1);
0822       m_C = m_fLamination[l][s]->GetParameter(2);
0823       m_A_err = m_fLamination[l][s]->GetParError(0);
0824       m_B_err = m_fLamination[l][s]->GetParError(1);
0825       m_C_err = m_fLamination[l][s]->GetParError(2);
0826       m_dist = m_distanceToFit[l][s];
0827       m_nBins = m_nBinsFit[l][s];
0828       m_laminationTree->Fill();
0829     }
0830   }
0831
0832   TFile *outputfile = new TFile(m_outputfile.c_str(), "RECREATE");
0833   outputfile->cd();
0834   for (int s = 0; s < 2; s++)
0835   {
0836     for (const auto &h : {m_dcc_out->m_hDRint[s], m_dcc_out->m_hDPint[s], m_dcc_out->m_hDZint[s], m_dcc_out->m_hentries[s]})
0837     {
0838       if (h)
0839       {
0840         h->Write();
0841       }
0842     }
0843     phiDistortionLamination[s]->Write();
0844     scaleFactorMap[s]->Write();
0845   }
0846   m_laminationTree->Write();
0847
0848   m_hLamination[13][0]->Write();
0849   m_hLamination[13][1]->Write();
0850   m_hLamination[14][1]->Write();
0851
0852
0853   outputfile->Close();
0854
0855   return Fun4AllReturnCodes::EVENT_OK;
0856 }
0857
0858 //________________________________
0859 void TpcLaminationFitting::fill_guarding_bins(TpcDistortionCorrectionContainer *dcc)
0860 {
0861   for (int s = 0; s < 2; s++)
0862   {
0863     for (const auto &h : {dcc->m_hDRint[s], dcc->m_hDPint[s]})
0864     {
0865       const auto phibins = h->GetNbinsX();
0866       const auto rbins = h->GetNbinsY();
0867
0868       int r30Bin = h->GetYaxis()->FindBin(30.0) + 1;
0869
0870       for (int ir = 0; ir < rbins; ir++)
0871       {
0872         h->SetBinContent(1, ir + 1, h->GetBinContent(phibins - 1, ir + 1));
0873         h->SetBinError(1, ir + 1, h->GetBinError(phibins - 1, ir + 1));
0874
0875         h->SetBinContent(phibins, ir + 1, h->GetBinContent(2, ir + 1));
0876         h->SetBinError(phibins, ir + 1, h->GetBinError(2, ir + 1));
0877       }
0878
0879       for (int iphi = 0; iphi < phibins; iphi++)
0880       {
0881         h->SetBinContent(iphi + 1, r30Bin-1, h->GetBinContent(iphi + 1, r30Bin));
0882         h->SetBinError(iphi + 1, r30Bin-1, h->GetBinError(iphi + 1, r30Bin));
0883
0884         h->SetBinContent(iphi + 1, rbins, h->GetBinContent(iphi + 1, rbins - 1));
0885         h->SetBinError(iphi + 1, rbins, h->GetBinError(iphi + 1, rbins - 1));
0886       }
0887     }
0888   }
0889 }
0890
0891
0892 //_____________________________
0893 void TpcLaminationFitting::setOutputfile(const std::string &outputfile)
0894 {
0895   m_outputfile = outputfile;
0896 }