Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:13:28 

0001 /*!
0002  *  \file     AnaMvtxPrototype1.C
0003  *  \brief    Analyze Mvtx 4 chip telescope from 2018 fermilab beam test
0004  *  \author   Darren McGlinchey
0005  */
0006 
0007 #include "AnaMvtxPrototype1.h"
0008 
0009 #include <trackbase/TrkrDefUtil.h>
0010 #include <trackbase/TrkrCluster.h>
0011 #include <trackbase/TrkrClusterContainer.h>
0012 
0013 #include <mvtxprototype1/MvtxStandaloneTracking.h>
0014 
0015 #include <phool/getClass.h>
0016 #include <phool/phool.h>
0017 #include <phool/PHCompositeNode.h>
0018 #include <phool/PHIODataNode.h>
0019 #include <phool/PHNodeIterator.h>
0020 
0021 #include <fun4all/Fun4AllServer.h>
0022 #include <fun4all/Fun4AllReturnCodes.h>
0023 #include <fun4all/PHTFileServer.h>
0024 
0025 #include <TTree.h>
0026 #include <TFile.h>
0027 #include <TLorentzVector.h>
0028 #include <TH1.h>
0029 #include <TH2.h>
0030 #include <TVectorD.h>
0031 #include <TMatrixD.h>
0032 #include <TDecompSVD.h>
0033 
0034 #include <iostream>
0035 #include <map>
0036 #include <utility>
0037 #include <vector>
0038 #include <memory>
0039 
0040 
0041 
0042 #define LogDebug(exp)   std::cout<<"DEBUG: "  <<__FILE__<<": "<<__LINE__<<": "<< exp <<"\n"
0043 #define LogError(exp)   std::cout<<"ERROR: "  <<__FILE__<<": "<<__LINE__<<": "<< exp <<"\n"
0044 #define LogWarning(exp) std::cout<<"WARNING: "  <<__FILE__<<": "<<__LINE__<<": "<< exp <<"\n"
0045 
0046 
0047 AnaMvtxPrototype1::AnaMvtxPrototype1(const std::string &name,
0048                                      const std::string &ofName) :
0049   SubsysReco(name),
0050   clusters_(NULL),
0051   _foutname(ofName),
0052   _f(NULL),
0053   _ievent(0)
0054 {
0055 
0056   hlayer = NULL;
0057   for (int il = 0; il < 4; il++)
0058   {
0059     hsize[il] = NULL;
0060     hsize_phi[il] = NULL;
0061     hsize_z[il] = NULL;
0062     hxz[il] = NULL;
0063     hdx[il] = NULL;
0064     hdz[il] = NULL;
0065   }
0066 
0067   mvtxtracking_ = new MvtxStandaloneTracking();
0068   mvtxtracking_->SetWindowX(10);
0069   mvtxtracking_->SetWindowZ(10);
0070   mvtxtracking_->Verbosity(0);
0071 }
0072 
0073 int AnaMvtxPrototype1::Init(PHCompositeNode *topNode)
0074 {
0075 
0076   //-- Open file
0077   _f = new TFile(_foutname.c_str(), "RECREATE");
0078 
0079   //-- Initialize histograms
0080   hlayer = new TH1D("hlayer", ";layer", 6, -0.5, 5.5);
0081 
0082   for (int il = 0; il < 4; il++)
0083   {
0084     hsize[il] = new TH1D(Form("hsize_l%i", il),
0085                          ";cluster size total [pixels]",
0086                          100, -0.5, 99.5);
0087 
0088     hsize_phi[il] = new TH1D(Form("hsize_phi_l%i", il),
0089                              ";cluster size #phi [pixels]",
0090                              100, -0.5, 99.5);
0091 
0092     hsize_z[il] = new TH1D(Form("hsize_z_l%i", il),
0093                            ";cluster size z [pixels]",
0094                            100, -0.5, 100.5);
0095 
0096     hxz[il] = new TH2D(Form("hxz_l%i", il),
0097                        ";z;x",
0098                        1024, -0.5, 1023.5,
0099                        512, -0.5, 511.5);
0100 
0101     hdx[il] = new TH1D(Form("hdx_l%i", il),
0102                        "; dx [pixels]",
0103                        2000, -500, 500);
0104 
0105     hdz[il] = new TH1D(Form("hdz_l%i", il),
0106                        "; dz [pixels]",
0107                        2000, -500, 500);
0108   } // il
0109 
0110   //-- results from tracking
0111   htrk = new TH1D("htrk", ";trks / event", 16, -0.5, 15.5);
0112   htrk_cut = new TH1D("htrk_cut", ";trks / event", 16, -0.5, 15.5);
0113 
0114   for (int il = 0; il < 4; il++)
0115   {
0116     htrk_dx[il] = new TH1D(Form("htrk_dx_l%i", il),
0117                            ";track dx [pixels]",
0118                            500, -25, 25);
0119 
0120     htrk_dz[il] = new TH1D(Form("htrk_dz_l%i", il),
0121                            ";track dz [pixels]",
0122                            500, -25, 25);
0123 
0124     htrk_cut_dx[il] = new TH1D(Form("htrk_cut_dx_l%i", il),
0125                                ";track dx [pixels]",
0126                                500, -25, 25);
0127 
0128     htrk_cut_dz[il] = new TH1D(Form("htrk_cut_dz_l%i", il),
0129                                ";track dz [pixels]",
0130                                500, -25, 25);
0131 
0132   }
0133   htrk_chi2xy = new TH1D("htrk_chi2xy",
0134                          ";track chi2/ndf in x vs y",
0135                          500, 0, 100);
0136   htrk_chi2zy = new TH1D("htrk_chi2zy",
0137                          ";track chi2/ndf in z vs y",
0138                          500, 0, 100);
0139 
0140   htrk_cut_chi2xy = new TH1D("htrk_cut_chi2xy",
0141                              ";track chi2/ndf in x vs y",
0142                              500, 0, 100);
0143   htrk_cut_chi2zy = new TH1D("htrk_cut_chi2zy",
0144                              ";track chi2/ndf in z vs y",
0145                              500, 0, 100);
0146 
0147   return 0;
0148 
0149 }
0150 
0151 int AnaMvtxPrototype1::InitRun(PHCompositeNode *topNode)
0152 {
0153 
0154   //-- set counters
0155   _ievent = -1; // since we count at the beginning of the event, start at -1
0156 
0157 
0158   return 0;
0159 
0160 }
0161 
0162 int AnaMvtxPrototype1::process_event(PHCompositeNode *topNode)
0163 {
0164 
0165   if ( verbosity > VERBOSITY_MORE )
0166     std::cout << "AnaMvtxPrototype1::process_event()" << std::endl;
0167 
0168   // count event
0169   _ievent++;
0170 
0171   //------
0172   //-- Get the nodes
0173   //------
0174   int retval = GetNodes(topNode);
0175   if ( retval != Fun4AllReturnCodes::EVENT_OK )
0176     return retval;
0177 
0178   if ( verbosity > VERBOSITY_MORE )
0179   {
0180     std::cout << "-- evnt:" << _ievent << std::endl;
0181     std::cout << "-- Nclus:" << clusters_->size() << std::endl;
0182   }
0183 
0184   //------
0185   //-- Misalign clusters
0186   //------
0187   // Misalign();
0188 
0189   //------
0190   //-- Loop over clusters
0191   //------
0192 
0193   if ( verbosity > VERBOSITY_MORE )
0194     std::cout << "-- Looping over clusters" << std::endl;
0195 
0196   LyrClusMap _mmap_lyr_clus;
0197 
0198   TrkrDefUtil trkrutil;
0199 
0200   TrkrClusterContainer::ConstRange clusrange = clusters_->GetClusters();
0201   for (TrkrClusterContainer::ConstIterator iter = clusrange.first;
0202        iter != clusrange.second;
0203        ++iter)
0204   {
0205     TrkrCluster *clus = iter->second;
0206     // clus->identify();
0207 
0208     TrkrDefs::cluskey ckey = clus->GetClusKey();
0209 
0210     int lyr = trkrutil.GetLayer(ckey);
0211 
0212     hlayer->Fill(lyr);
0213 
0214     if (lyr < 0 || lyr > 3)
0215       continue;
0216     hsize[lyr]->Fill(clus->GetAdc());
0217     hsize_phi[lyr]->Fill(clus->GetPhiSize());
0218     hsize_z[lyr]->Fill(clus->GetZSize());
0219     hxz[lyr]->Fill(clus->GetZ(), clus->GetX());
0220 
0221     // insert into map
0222     _mmap_lyr_clus.insert(std::make_pair(lyr, clus));
0223   }
0224 
0225   if ( verbosity > VERBOSITY_MORE )
0226     std::cout << "-- _mmap_lyr_clus.size():" << _mmap_lyr_clus.size() << std::endl;
0227 
0228   //------
0229   //-- Simple tracking
0230   //------
0231 
0232   if ( verbosity > VERBOSITY_MORE )
0233     std::cout << "-- Performing simple tracking" << std::endl;
0234 
0235   static bool check_range = true;
0236   if ( check_range )
0237   {
0238     TrkrDefs::cluskey mvtxl = trkrutil.GetClusKeyLo(TrkrDefs::TRKRID::mvtx_id);
0239     TrkrDefs::cluskey mvtxh = trkrutil.GetClusKeyHi(TrkrDefs::TRKRID::mvtx_id);
0240 
0241     std::cout << PHWHERE << " Mvtx key range: 0x" << std::hex << mvtxl << " - 0x" << mvtxh << std::dec << std::endl;
0242 
0243     for (int i = 0; i < 4; i++)
0244     {
0245       TrkrDefs::cluskey lyrl = trkrutil.GetClusKeyLo(TrkrDefs::TRKRID::mvtx_id, i);
0246       TrkrDefs::cluskey lyrh = trkrutil.GetClusKeyHi(TrkrDefs::TRKRID::mvtx_id, i);
0247 
0248       std::cout << PHWHERE << " Mvtx key range lyr " << i << ": 0x" << std::hex << lyrl << " - 0x" << lyrh << std::dec << std::endl;
0249     }
0250 
0251     check_range = false;
0252   }
0253 
0254 
0255   // loop over clusters in layer 0
0256   TrkrClusterContainer::ConstRange lyr0range = clusters_->GetClusters(TrkrDefs::TRKRID::mvtx_id, 0);
0257   for ( TrkrClusterContainer::ConstIterator lyr0_itr = lyr0range.first;
0258         lyr0_itr != lyr0range.second;
0259         ++lyr0_itr)
0260   {
0261     TrkrCluster* clus0 = lyr0_itr->second;
0262 
0263     if ( !clus0 )
0264     {
0265       std::cout << "WARNING: clus0 not found" << std::endl;
0266       continue;
0267     }
0268 
0269     // check layer
0270     int lyr = trkrutil.GetLayer(clus0->GetClusKey());
0271     if ( lyr != 0 )
0272     {
0273       std::cout << PHWHERE << "ERRORL: Expected clusters from layer 0 only. Found lyr=" << lyr << std::endl;
0274       continue;
0275     }
0276 
0277 
0278 
0279     auto misiter = _misalign.find(0);
0280     // get (x, y, z)
0281     double x0 = misiter != _misalign.end() ? clus0->GetX() + (misiter->second).dx : clus0->GetX();
0282     double y0 = misiter != _misalign.end() ? clus0->GetY() + (misiter->second).dy : clus0->GetY();
0283     double z0 = misiter != _misalign.end() ? clus0->GetZ() + (misiter->second).dz : clus0->GetZ();
0284 
0285     if ( verbosity > VERBOSITY_MORE )
0286     {
0287       std::cout << "-- clus0:(" << x0 << ", " << y0 << ", " << z0 << ")" << std::endl;
0288     }
0289 
0290     // loop over all clusters in the outer layer
0291     TrkrClusterContainer::ConstRange lyr3range = clusters_->GetClusters(TrkrDefs::TRKRID::mvtx_id, 3);
0292     for ( TrkrClusterContainer::ConstIterator lyr3_itr = lyr3range.first;
0293           lyr3_itr != lyr3range.second;
0294           ++lyr3_itr)
0295     {
0296       TrkrCluster* clus3 = lyr3_itr->second;
0297 
0298       if ( !clus3 )
0299       {
0300         std::cout << "WARNING: clus3 not found" << std::endl;
0301         continue;
0302       }
0303 
0304       auto misiter = _misalign.find(3);
0305       // get (x, y, z)
0306       double x3 = misiter != _misalign.end() ? clus3->GetX() + (misiter->second).dx : clus3->GetX();
0307       double y3 = misiter != _misalign.end() ? clus3->GetY() + (misiter->second).dy : clus3->GetY();
0308       double z3 = misiter != _misalign.end() ? clus3->GetZ() + (misiter->second).dz : clus3->GetZ();
0309 
0310 
0311       // calc m, b in (phi, r) space
0312       double mpr = CalcSlope(y0, x0, y3, x3);
0313       double bpr = CalcIntecept(y0, x0, mpr);
0314 
0315       // calc m, b in (z, r) space
0316       double mzr = CalcSlope(y0, z0, y3, z3);
0317       double bzr = CalcIntecept(y0, z0, mzr);
0318 
0319 
0320       if ( verbosity > VERBOSITY_MORE )
0321       {
0322         std::cout << "-- clus3:(" << x3 << ", " << y3 << ", " << z3 << ")" << std::endl;
0323         std::cout << "-- mpr:" << mpr << " bpr:" << bpr << std::endl;
0324         std::cout << "-- mzr:" << mzr << " bzr:" << bzr << std::endl;
0325       }
0326 
0327 
0328 
0329       // loop over all clusters in layer 1 & 2 and calc displacement
0330       for (int ilyr = 1; ilyr <= 2; ilyr++)
0331       {
0332         TrkrClusterContainer::ConstRange lyrrange = clusters_->GetClusters(TrkrDefs::TRKRID::mvtx_id, ilyr);
0333         for ( TrkrClusterContainer::ConstIterator lyr_itr = lyrrange.first;
0334               lyr_itr != lyrrange.second;
0335               ++lyr_itr)
0336         {
0337           TrkrCluster* clus = lyr_itr->second;
0338 
0339           if ( !clus )
0340             continue;
0341 
0342           auto misiter = _misalign.find(ilyr);
0343           // get (x, y, z)
0344           double x = misiter != _misalign.end() ? clus->GetX() + (misiter->second).dx : clus->GetX();
0345           double y = misiter != _misalign.end() ? clus->GetY() + (misiter->second).dy : clus->GetY();
0346           double z = misiter != _misalign.end() ? clus->GetZ() + (misiter->second).dz : clus->GetZ();
0347 
0348 
0349           // calc track projection in phi, z
0350           double x_proj = CalcProjection(y, mpr, bpr);
0351           double z_proj = CalcProjection(y, mzr, bzr);
0352 
0353           // calc proj - true
0354           double dx = (x_proj - x);
0355           double dz = (z_proj - z);
0356 
0357           if ( !hdx[ilyr] || !hdz[ilyr] )
0358           {
0359             std::cout << "WARNING!! can't find hdphi or hdz for lyr:" << ilyr << std::endl;
0360           }
0361           else
0362           {
0363             hdx[ilyr]->Fill(dx);
0364             hdz[ilyr]->Fill(dz);
0365           }
0366 
0367           if ( verbosity > VERBOSITY_MORE )
0368           {
0369             std::cout << "------------------" << std::endl;
0370             std::cout << "-- clus" << trkrutil.GetLayer(clus->GetClusKey()) << ":(" << x3 << ", " << y3 << ", " << z3 << ")" << std::endl;
0371             std::cout << "-- proj:(" << x_proj << ", " << z_proj << ")" << std::endl;
0372             std::cout << "-- dx [pixels]: " << dx << std::endl;
0373             std::cout << "-- dz [pixels]: " << dz << std::endl;
0374             std::cout << "------------------" << std::endl;
0375           }
0376 
0377         } // lyr1_itr
0378       } // ilyr;
0379 
0380     } // lyr3_itr
0381 
0382   } //lyr0_itr
0383 
0384   // // loop over clusters in layer 0
0385   // LyrClusMap::iterator lyr0_itr;
0386   // LyrClusMap::iterator lyr_itr;
0387   // LyrClusMap::iterator lyr3_itr;
0388   // for ( lyr0_itr = _mmap_lyr_clus.lower_bound(0);
0389   //       lyr0_itr != _mmap_lyr_clus.upper_bound(0);
0390   //       ++lyr0_itr)
0391   // {
0392   //   TrkrCluster* clus0 = lyr0_itr->second;
0393 
0394   //   if ( !clus0 )
0395   //   {
0396   //     std::cout << "WARNING: clus0 not found" << std::endl;
0397   //     continue;
0398   //   }
0399   //   auto misiter = _misalign.find(0);
0400   //   // get (x, y, z)
0401   //   double x0 = misiter != _misalign.end() ? clus0->GetX() + (misiter->second).dx : clus0->GetX();
0402   //   double y0 = misiter != _misalign.end() ? clus0->GetY() + (misiter->second).dy : clus0->GetY();
0403   //   double z0 = misiter != _misalign.end() ? clus0->GetZ() + (misiter->second).dz : clus0->GetZ();
0404 
0405   //   if ( verbosity > VERBOSITY_MORE )
0406   //   {
0407   //     std::cout << "-- clus0:(" << x0 << ", " << y0 << ", " << z0 << ")" << std::endl;
0408   //   }
0409 
0410   //   // loop over all clusters in the outer layer
0411   //   for ( lyr3_itr = _mmap_lyr_clus.lower_bound(3);
0412   //         lyr3_itr != _mmap_lyr_clus.upper_bound(3);
0413   //         ++lyr3_itr)
0414   //   {
0415   //     TrkrCluster* clus3 = lyr3_itr->second;
0416 
0417   //     if ( !clus3 )
0418   //     {
0419   //       std::cout << "WARNING: clus3 not found" << std::endl;
0420   //       continue;
0421   //     }
0422 
0423   //     auto misiter = _misalign.find(3);
0424   //     // get (x, y, z)
0425   //     double x3 = misiter != _misalign.end() ? clus3->GetX() + (misiter->second).dx : clus3->GetX();
0426   //     double y3 = misiter != _misalign.end() ? clus3->GetY() + (misiter->second).dy : clus3->GetY();
0427   //     double z3 = misiter != _misalign.end() ? clus3->GetZ() + (misiter->second).dz : clus3->GetZ();
0428 
0429 
0430   //     // calc m, b in (phi, r) space
0431   //     double mpr = CalcSlope(y0, x0, y3, x3);
0432   //     double bpr = CalcIntecept(y0, x0, mpr);
0433 
0434   //     // calc m, b in (z, r) space
0435   //     double mzr = CalcSlope(y0, z0, y3, z3);
0436   //     double bzr = CalcIntecept(y0, z0, mzr);
0437 
0438 
0439   //     if ( verbosity > VERBOSITY_MORE )
0440   //     {
0441   //       std::cout << "-- clus3:(" << x3 << ", " << y3 << ", " << z3 << ")" << std::endl;
0442   //       std::cout << "-- mpr:" << mpr << " bpr:" << bpr << std::endl;
0443   //       std::cout << "-- mzr:" << mzr << " bzr:" << bzr << std::endl;
0444   //     }
0445 
0446 
0447 
0448   //     // loop over all clusters in layer 1 & 2 and calc displacement
0449   //     for (int ilyr = 1; ilyr <= 2; ilyr++)
0450   //     {
0451   //       for ( lyr_itr = _mmap_lyr_clus.lower_bound(ilyr);
0452   //             lyr_itr != _mmap_lyr_clus.upper_bound(ilyr);
0453   //             ++lyr_itr)
0454   //       {
0455   //         TrkrCluster* clus = lyr_itr->second;
0456 
0457   //         if ( !clus )
0458   //           continue;
0459 
0460   //         auto misiter = _misalign.find(ilyr);
0461   //         // get (x, y, z)
0462   //         double x = misiter != _misalign.end() ? clus->GetX() + (misiter->second).dx : clus->GetX();
0463   //         double y = misiter != _misalign.end() ? clus->GetY() + (misiter->second).dy : clus->GetY();
0464   //         double z = misiter != _misalign.end() ? clus->GetZ() + (misiter->second).dz : clus->GetZ();
0465 
0466 
0467   //         // calc track projection in phi, z
0468   //         double x_proj = CalcProjection(y, mpr, bpr);
0469   //         double z_proj = CalcProjection(y, mzr, bzr);
0470 
0471   //         // calc proj - true
0472   //         double dx = (x_proj - x);
0473   //         double dz = (z_proj - z);
0474 
0475   //         if ( !hdx[ilyr] || !hdz[ilyr] )
0476   //         {
0477   //           std::cout << "WARNING!! can't find hdphi or hdz for lyr:" << ilyr << std::endl;
0478   //         }
0479   //         else
0480   //         {
0481   //           hdx[ilyr]->Fill(dx);
0482   //           hdz[ilyr]->Fill(dz);
0483   //         }
0484 
0485   //         if ( verbosity > VERBOSITY_MORE )
0486   //         {
0487   //           std::cout << "------------------" << std::endl;
0488   //           std::cout << "-- clus" << trkrutil.GetLayer(clus->GetClusKey()) << ":(" << x3 << ", " << y3 << ", " << z3 << ")" << std::endl;
0489   //           std::cout << "-- proj:(" << x_proj << ", " << z_proj << ")" << std::endl;
0490   //           std::cout << "-- dx [pixels]: " << dx << std::endl;
0491   //           std::cout << "-- dz [pixels]: " << dz << std::endl;
0492   //           std::cout << "------------------" << std::endl;
0493   //         }
0494 
0495   //       } // lyr1_itr
0496   //     } // ilyr;
0497 
0498   //   } // lyr3_itr
0499 
0500   // } //lyr0_itr
0501 
0502 
0503   if ( verbosity > VERBOSITY_MORE )
0504     std::cout << "-- Done simple tracking" << std::endl;
0505 
0506 
0507 
0508   // //------
0509   // //-- Track candidates using all layers
0510   // //------
0511 
0512   // LyrClusMap::iterator lyr_itr;
0513   // TrkVec _vec_trk;
0514   // ClusVec _trk_candidate;
0515 
0516   // for ( lyr_itr = _mmap_lyr_clus.lower_bound(0);
0517   //       lyr_itr != _mmap_lyr_clus.upper_bound(0);
0518   //       ++lyr_itr)
0519   // {
0520   //   TrkrCluster* clus0 = lyr_itr->second;
0521 
0522   // }
0523 
0524 
0525   //------
0526   // Try full tracking
0527   //------
0528   MvtxStandaloneTracking::MvtxTrackList tracklist;
0529 
0530   std::vector<int> uselayers;
0531   uselayers.push_back(0);
0532   uselayers.push_back(1);
0533   uselayers.push_back(2);
0534   uselayers.push_back(3);
0535 
0536   mvtxtracking_->RunTracking(topNode, tracklist, uselayers);
0537 
0538   htrk->Fill(tracklist.size());
0539 
0540   int ncut = 0;
0541   for ( unsigned int i = 0; i < tracklist.size(); i++)
0542   {
0543 
0544     htrk_chi2xy->Fill(tracklist.at(i).chi2_xy);
0545     htrk_chi2zy->Fill(tracklist.at(i).chi2_zy);
0546 
0547     if ( tracklist.at(i).chi2_xy < 5 && tracklist.at(i).chi2_zy < 5)
0548     {
0549       ncut++;
0550       htrk_cut_chi2xy->Fill(tracklist.at(i).chi2_xy);
0551       htrk_cut_chi2zy->Fill(tracklist.at(i).chi2_zy);
0552     }
0553 
0554     for ( unsigned int j = 0; j < tracklist.at(i).ClusterList.size(); j++)
0555     {
0556       TrkrDefs::cluskey ckey = tracklist.at(i).ClusterList.at(j)->GetClusKey();
0557 
0558       int lyr = trkrutil.GetLayer(ckey);
0559 
0560       if ( lyr < 0 || lyr >= 4 )
0561       {
0562         std::cout << PHWHERE << " WARNING: bad layer from track cluster. lyr:" << lyr << std::endl;
0563         continue;
0564       }
0565 
0566       htrk_dx[lyr]->Fill(tracklist.at(i).dx.at(j));
0567       htrk_dz[lyr]->Fill(tracklist.at(i).dz.at(j));
0568 
0569       if ( tracklist.at(i).chi2_xy < 5 && tracklist.at(i).chi2_zy < 5)
0570       {
0571         htrk_cut_dx[lyr]->Fill(tracklist.at(i).dx.at(j));
0572         htrk_cut_dz[lyr]->Fill(tracklist.at(i).dz.at(j));
0573       }
0574     } // j
0575 
0576 
0577   } // i
0578 
0579   htrk_cut->Fill(ncut);
0580   
0581   //-- Cleanup
0582 
0583 
0584   //-- End
0585   return 0;
0586 
0587 }
0588 
0589 int AnaMvtxPrototype1::End(PHCompositeNode * topNode)
0590 {
0591 
0592   PrintMisalign();
0593 
0594   _f->Write();
0595   _f->Close();
0596 
0597   return 0;
0598 }
0599 
0600 /*
0601  * Get all the necessary nodes from the node tree
0602  */
0603 int AnaMvtxPrototype1::GetNodes(PHCompositeNode * topNode)
0604 {
0605 
0606   // Input Svtx Clusters
0607   clusters_ = findNode::getClass<TrkrClusterContainer>(topNode, "TrkrClusterContainer");
0608   if (!clusters_ && _ievent < 2)
0609   {
0610     std::cout << PHWHERE << " TrkrClusterContainer node not found on node tree"
0611               << std::endl;
0612     return Fun4AllReturnCodes::ABORTEVENT;
0613   }
0614 
0615   return Fun4AllReturnCodes::EVENT_OK;
0616 
0617 }
0618 
0619 /*
0620  * Simple helper function for calculating the slope between two points
0621  */
0622 double
0623 AnaMvtxPrototype1::CalcSlope(double x0, double y0, double x1, double y1)
0624 {
0625   return (y1 - y0) / (x1 - x0);
0626 }
0627 
0628 /*
0629  * Simple helper function for calculating intercept
0630  */
0631 double
0632 AnaMvtxPrototype1::CalcIntecept(double x0, double y0, double m)
0633 {
0634   return y0 - x0 * m;
0635 }
0636 
0637 /*
0638  * Simple helper function for calculating projection
0639  */
0640 double
0641 AnaMvtxPrototype1::CalcProjection(double x, double m, double b)
0642 {
0643   return m * x + b;
0644 }
0645 
0646 /*
0647  *
0648  */
0649 void
0650 AnaMvtxPrototype1::MisalignLayer(int lyr, double dx, double dy, double dz)
0651 {
0652   mis align;
0653   align.dx = dx;
0654   align.dy = dy;
0655   align.dz = dz;
0656 
0657   _misalign.insert( std::make_pair(lyr, align) );
0658 }
0659 
0660 /*
0661  * Misalign clusters for each layer
0662  */
0663 void
0664 AnaMvtxPrototype1::Misalign()
0665 {
0666   // only bother if we've set some misalignments to non-zero values
0667   if ( _misalign.size() < 1 )
0668     return;
0669 
0670   TrkrClusterContainer::ConstRange clusrange = clusters_->GetClusters();
0671   for (TrkrClusterContainer::ConstIterator iter = clusrange.first;
0672        iter != clusrange.second;
0673        ++iter)
0674   {
0675     // TrkrCluster *clus = iter->second;
0676     /*
0677     // get misalignments
0678     auto misitr = _misalign.find(clus->get_layer());
0679     if ( misitr != _misalign.end() )
0680     {
0681       mis align = misitr->second;
0682 
0683       float x = clus->GetX();
0684       float y = clus->GetY();
0685       float z = clus->GetZ();
0686       float r = TMath::Sqrt(x * x + y * y);
0687       float phi = TMath::ATan2(y, x);
0688       float s = r * phi;
0689 
0690       double zp = z;
0691       double rp = r;
0692       double sp = s;
0693 
0694       // --- pitch, yaw, roll
0695       // pitch is dr(z)
0696       rp += z * TMath::Sin(align.dpitch);
0697 
0698       // yaw is ds(z)
0699       sp += z * TMath::Sin(align.dyaw);
0700 
0701       // roll is dr(s)
0702       rp += s * TMath::Sin(align.droll);
0703 
0704       // --- translate (s, z, r)
0705       rp += align.dr;
0706       sp += align.ds;
0707       zp += align.dz;
0708 
0709       // --- recalculat x', y', z'
0710       float phip = sp / rp;
0711       float xp = rp * TMath::Cos(phip);
0712       float yp = rp * TMath::Sin(phip);
0713 
0714       // --- set new coordinates
0715       clus->set_x(xp);
0716       clus->set_y(yp);
0717       clus->set_z(zp);
0718 
0719       if ( false )
0720       {
0721         std::cout << " lyr:" << clus->get_layer()
0722                   << " ds:" << align.ds
0723                   << " dz:" << align.dz
0724                   << " dr:" << align.dr
0725                   << " dpitch:" << align.dpitch
0726                   << " dyaw:" << align.dyaw
0727                   << " droll:" << align.droll
0728                   << " phi:" << phi << " phip:" << phip
0729                   << " z:" << z << " zp:" << zp
0730                   << " x:" << x << " xp:" << xp
0731                   << " y:" << y << " yp:" << yp
0732                   << std::endl;
0733       }
0734 
0735     }
0736     */
0737   } // iter
0738 }
0739 
0740 /*
0741  *
0742  */
0743 void
0744 AnaMvtxPrototype1::PrintMisalign()
0745 {
0746   std::cout << "================================" << std::endl;
0747   std::cout << "== " << PHWHERE << std::endl;
0748 
0749   for (auto iter = _misalign.begin(); iter != _misalign.end(); ++iter)
0750   {
0751     std::cout << "== Lyr: " << iter->first
0752               << " dx[cm]:" << (iter->second).dx
0753               << " dy[cm]:" << (iter->second).dy
0754               << " dz[cm]:" << (iter->second).dz
0755               << std::endl;
0756   } // iter
0757 
0758   std::cout << "================================" << std::endl;
0759 }
0760
[ 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!