MvtxTelescope/AnaMvtxTelescopeHits/AnaMvtxTelescopeHits.C

0001 /*!
0002  *  \file     AnaMvtxTelescopeHits.C
0003  *  \brief    Analyze Mvtx 4 ladder cosmic telescope in simulations
0004  *  \details  Analyze simulations of hits in the 4 ladder
0005  *            Mvtx cosimic ray telescope
0006  *  \author   Darren McGlinchey
0007  */
0008
0009 #include "AnaMvtxTelescopeHits.h"
0010
0011 #include <phool/getClass.h>
0012 #include <phool/phool.h>
0013 #include <phool/PHCompositeNode.h>
0014 #include <phool/PHIODataNode.h>
0015 #include <phool/PHNodeIterator.h>
0016
0017 #include <phgeom/PHGeomUtility.h>
0018
0019 #include <fun4all/Fun4AllServer.h>
0020 #include <fun4all/Fun4AllReturnCodes.h>
0021 #include <fun4all/PHTFileServer.h>
0022
0023 #include <g4main/PHG4TruthInfoContainer.h>
0024 #include <g4main/PHG4Particle.h>
0025 #include <g4main/PHG4Hit.h>
0026 #include <g4main/PHG4HitContainer.h>
0027
0028 #include <g4hough/SvtxCluster.h>
0029 #include <g4hough/SvtxClusterMap.h>
0030 #include <g4hough/SvtxHitMap.h>
0031 #include <g4hough/SvtxTrack.h>
0032 #include <g4hough/SvtxVertex.h>
0033 #include <g4hough/SvtxTrackMap.h>
0034 #include <g4hough/SvtxVertexMap.h>
0035
0036 #include <g4jets/JetMap.h>
0037 #include <g4jets/Jet.h>
0038
0039 #include <g4detectors/PHG4CellContainer.h>
0040 #include <g4detectors/PHG4CylinderGeomContainer.h>
0041 #include <g4detectors/PHG4Cell.h>
0042 #include <g4detectors/PHG4CylinderGeom_MAPS.h>
0043 #include <g4detectors/PHG4CylinderGeom_Siladders.h>
0044
0045 #include <phgenfit/Fitter.h>
0046 #include <phgenfit/PlanarMeasurement.h>
0047 #include <phgenfit/Track.h>
0048 #include <phgenfit/SpacepointMeasurement.h>
0049
0050 //GenFit
0051 #include <GenFit/FieldManager.h>
0052 #include <GenFit/GFRaveConverters.h>
0053 #include <GenFit/GFRaveVertex.h>
0054 #include <GenFit/GFRaveVertexFactory.h>
0055 #include <GenFit/MeasuredStateOnPlane.h>
0056 #include <GenFit/RKTrackRep.h>
0057 #include <GenFit/StateOnPlane.h>
0058 #include <GenFit/Track.h>
0059
0060 //Rave
0061 #include <rave/Version.h>
0062 #include <rave/Track.h>
0063 #include <rave/VertexFactory.h>
0064 #include <rave/ConstantMagneticField.h>
0065
0066 #include <phhepmc/PHHepMCGenEvent.h>
0067 #include <HepMC/GenEvent.h>
0068 #include <HepMC/GenVertex.h>
0069
0070 #include <HFJetTruthGeneration/HFJetDefs.h>
0071
0072 #include "TTree.h"
0073 #include "TFile.h"
0074 #include "TLorentzVector.h"
0075 #include "TH1.h"
0076 #include "TH2.h"
0077 #include "TVectorD.h"
0078 #include "TMatrixD.h"
0079 #include "TDecompSVD.h"
0080
0081
0082 #include <iostream>
0083 #include <map>
0084 #include <utility>
0085 #include <vector>
0086 #include <memory>
0087
0088
0089
0090 #define LogDebug(exp)   std::cout<<"DEBUG: "  <<__FILE__<<": "<<__LINE__<<": "<< exp <<"\n"
0091 #define LogError(exp)   std::cout<<"ERROR: "  <<__FILE__<<": "<<__LINE__<<": "<< exp <<"\n"
0092 #define LogWarning(exp) std::cout<<"WARNING: "  <<__FILE__<<": "<<__LINE__<<": "<< exp <<"\n"
0093
0094
0095 AnaMvtxTelescopeHits::AnaMvtxTelescopeHits(const std::string &name,
0096     const std::string &ofName) :
0097   SubsysReco(name),
0098   _clustermap(NULL),
0099   _foutname(ofName),
0100   _f(NULL)
0101 {
0102
0103   hlayer = NULL;
0104   for (int il = 0; il < 4; il++)
0105   {
0106     hsize_phi[il] = NULL;
0107     hsize_z[il] = NULL;
0108     hphiz[il] = NULL;
0109     hdphi[il] = NULL;
0110     hdz[il] = NULL;
0111   }
0112 }
0113
0114 int AnaMvtxTelescopeHits::Init(PHCompositeNode *topNode)
0115 {
0116
0117   //-- Open file
0118   _f = new TFile(_foutname.c_str(), "RECREATE");
0119
0120   //-- Initialize histograms
0121   hlayer = new TH1D("hlayer", ";layer", 6, -0.5, 5.5);
0122
0123   for (int il = 0; il < 4; il++)
0124   {
0125     hsize_phi[il] = new TH1D(Form("hsize_phi_l%i", il),
0126                              ";cluster size #phi [cm]",
0127                              100, 0, 0.1);
0128
0129     hsize_z[il] = new TH1D(Form("hsize_z_l%i", il),
0130                            ";cluster size z [cm]",
0131                            100, 0, 0.1);
0132
0133     hphiz[il] = new TH2D(Form("hphiz_l%i", il),
0134                          ";z [cm]; #phi",
0135                          1000, -10, 10,
0136                          100, -0.1, 0.1);
0137
0138     hdphi[il] = new TH1D(Form("hdphi_l%i", il),
0139                          "; d#phi [#mum]",
0140                          2000, -100, 100);
0141
0142     hdz[il] = new TH1D(Form("hdz_l%i", il),
0143                          "; dz [#mum]",
0144                          2000, -100, 100);
0145   } // il
0146
0147
0148
0149   return 0;
0150
0151 }
0152
0153 int AnaMvtxTelescopeHits::InitRun(PHCompositeNode *topNode)
0154 {
0155
0156   // TGeoManager* tgeo_manager = PHGeomUtility::GetTGeoManager(topNode);
0157
0158   // _fitter = PHGenFit::Fitter::getInstance(tgeo_manager, _mag_field_file_name.data(),
0159   //                                         (_reverse_mag_field) ? -1. * _mag_field_re_scaling_factor : _mag_field_re_scaling_factor,
0160   //                                         _track_fitting_alg_name, "RKTrackRep", _do_evt_display);
0161
0162   // if (!_fitter) {
0163   //   std::cerr << PHWHERE << std::endl;
0164   //   return Fun4AllReturnCodes::ABORTRUN;
0165   // }
0166
0167   //-- set counters
0168   _ievent = -1; // since we count at the beginning of the event, start at -1
0169
0170
0171   return 0;
0172
0173 }
0174
0175 int AnaMvtxTelescopeHits::process_event(PHCompositeNode *topNode)
0176 {
0177
0178   if ( verbosity > VERBOSITY_MORE )
0179     std::cout << "AnaMvtxTelescopeHits::process_event()" << std::endl;
0180
0181   // count event
0182   _ievent++;
0183
0184   //------
0185   //-- Get the nodes
0186   //------
0187   int retval = GetNodes(topNode);
0188   if ( retval != Fun4AllReturnCodes::EVENT_OK )
0189     return retval;
0190
0191   if ( verbosity > VERBOSITY_MORE )
0192   {
0193     std::cout << "-- evnt:" << _ievent << std::endl;
0194     std::cout << "-- Nclus:" << _clustermap->size() << std::endl;
0195   }
0196
0197
0198   //------
0199   //-- Loop over clusters
0200   //------
0201
0202   if ( verbosity > VERBOSITY_MORE )
0203     std::cout << "-- Looping over clusters" << std::endl;
0204
0205   LyrClusMap _mmap_lyr_clus;
0206
0207   for (SvtxClusterMap::Iter iter = _clustermap->begin();
0208        iter != _clustermap->end();
0209        ++iter)
0210   {
0211     SvtxCluster *clus = iter->second;
0212
0213     hlayer->Fill(clus->get_layer());
0214
0215     unsigned int lyr = clus->get_layer();
0216     if (lyr < 0 || lyr > 3)
0217       continue;
0218     hsize_phi[lyr]->Fill(clus->get_phi_size());
0219     hsize_z[lyr]->Fill(clus->get_z_size());
0220
0221     double phi = TMath::ATan2(clus->get_y(), clus->get_x());
0222     hphiz[lyr]->Fill(clus->get_z(), phi);
0223
0224     // insert into map
0225     _mmap_lyr_clus.insert(std::make_pair(lyr, clus));
0226   }
0227
0228   if ( verbosity > VERBOSITY_MORE )
0229     std::cout << "-- _mmap_lyr_clus.size():" << _mmap_lyr_clus.size() << std::endl;
0230
0231   //------
0232   //-- Simple tracking
0233   //------
0234
0235   if ( verbosity > VERBOSITY_MORE )
0236     std::cout << "-- Performing simple tracking" << std::endl;
0237
0238   // loop over clusters in layer 0
0239   LyrClusMap::iterator lyr0_itr;
0240   LyrClusMap::iterator lyr_itr;
0241   LyrClusMap::iterator lyr3_itr;
0242   for ( lyr0_itr = _mmap_lyr_clus.lower_bound(0);
0243         lyr0_itr != _mmap_lyr_clus.upper_bound(0);
0244         ++lyr0_itr)
0245   {
0246     SvtxCluster* clus0 = lyr0_itr->second;
0247
0248     if ( !clus0 )
0249     {
0250       std::cout << "WARNING: clus0 not found" << std::endl;
0251       continue;
0252     }
0253
0254     // get (r, phi, z)
0255     double r0 = TMath::Sqrt(TMath::Power(clus0->get_x(), 2) + TMath::Power(clus0->get_y(), 2));
0256     double p0 = TMath::ATan2(clus0->get_y(), clus0->get_x());
0257     double z0 = clus0->get_z();
0258
0259     if ( verbosity > VERBOSITY_MORE )
0260     {
0261       std::cout << "-- clus0:(" << r0 << ", " << p0 << ", " << z0 << ")" << std::endl;
0262     }
0263
0264     // loop over all clusters in the outer layer
0265     for ( lyr3_itr = _mmap_lyr_clus.lower_bound(3);
0266           lyr3_itr != _mmap_lyr_clus.upper_bound(3);
0267           ++lyr3_itr)
0268     {
0269       SvtxCluster* clus3 = lyr3_itr->second;
0270
0271       if ( !clus3 )
0272       {
0273         std::cout << "WARNING: clus3 not found" << std::endl;
0274         continue;
0275       }
0276
0277       // get (r, phi, z)
0278       double r3 = TMath::Sqrt(TMath::Power(clus3->get_x(), 2) + TMath::Power(clus3->get_y(), 2));
0279       double p3 = TMath::ATan2(clus3->get_y(), clus3->get_x());
0280       double z3 = clus3->get_z();
0281
0282       // calc m, b in (phi, r) space
0283       double mpr = CalcSlope(r0, p0, r3, p3);
0284       double bpr = CalcIntecept(r0, p0, mpr);
0285
0286       // calc m, b in (z, r) space
0287       double mzr = CalcSlope(r0, z0, r3, z3);
0288       double bzr = CalcIntecept(r0, z0, mzr);
0289
0290
0291       if ( verbosity > VERBOSITY_MORE )
0292       {
0293         std::cout << "-- clus3:(" << r3 << ", " << p3 << ", " << z3 << ")" << std::endl;
0294         std::cout << "-- mpr:" << mpr << " bpr:" << bpr << std::endl;
0295         std::cout << "-- mzr:" << mzr << " bzr:" << bzr << std::endl;
0296       }
0297
0298
0299
0300       // loop over all clusters in layer 1 & 2 and calc displacement
0301       for (int ilyr = 1; ilyr <= 2; ilyr++)
0302       {
0303         for ( lyr_itr = _mmap_lyr_clus.lower_bound(ilyr);
0304               lyr_itr != _mmap_lyr_clus.upper_bound(ilyr);
0305               ++lyr_itr)
0306         {
0307           SvtxCluster* clus = lyr_itr->second;
0308
0309           if ( !clus )
0310             continue;
0311
0312           // get (r, phi, z)
0313           double r = TMath::Sqrt(TMath::Power(clus->get_x(), 2) + TMath::Power(clus->get_y(), 2));
0314           double phi = TMath::ATan2(clus->get_y(), clus->get_x());
0315           double z = clus->get_z();
0316
0317           // calc track projection in phi, z
0318           double phi_proj = CalcProjection(r, mpr, bpr);
0319           double z_proj = CalcProjection(r, mzr, bzr);
0320
0321           // calc proj - true
0322           double dphi = (phi_proj - phi) * 1e4;
0323           double dz = (z_proj - z) * 1e4;
0324
0325           if ( !hdphi[ilyr] || !hdz[ilyr] )
0326           {
0327             std::cout << "WARNING!! can't find hdphi or hdz for lyr:" << ilyr << std::endl;
0328           }
0329           else
0330           {
0331             hdphi[ilyr]->Fill(dphi);
0332             hdz[ilyr]->Fill(dz);
0333           }
0334
0335           if ( verbosity > VERBOSITY_MORE )
0336           {
0337             std::cout << "------------------" << std::endl;
0338             std::cout << "-- clus" << clus->get_layer() << ":(" << r3 << ", " << p3 << ", " << z3 << ")" << std::endl;
0339             std::cout << "-- proj:(" << phi_proj << ", " << z_proj << ")" << std::endl;
0340             std::cout << "-- dphi [micron]: " << dphi << std::endl;
0341             std::cout << "-- dz [micron]:   " << dz << std::endl;
0342             std::cout << "------------------" << std::endl;
0343           }
0344
0345         } // lyr1_itr
0346       } // ilyr;
0347
0348     } // lyr3_itr
0349
0350   } //lyr0_itr
0351
0352
0353   if ( verbosity > VERBOSITY_MORE )
0354     std::cout << "-- Done simple tracking" << std::endl;
0355
0356
0357
0358   // //------
0359   // //-- Track candidates using all layers
0360   // //------
0361
0362   // LyrClusMap::iterator lyr_itr;
0363   // TrkVec _vec_trk;
0364   // ClusVec _trk_candidate;
0365
0366   // for ( lyr_itr = _mmap_lyr_clus.lower_bound(0);
0367   //       lyr_itr != _mmap_lyr_clus.upper_bound(0);
0368   //       ++lyr_itr)
0369   // {
0370   //   SvtxCluster* clus0 = lyr_itr->second;
0371
0372   // }
0373
0374
0375   //-- Cleanup
0376
0377
0378   //-- End
0379   return 0;
0380
0381 }
0382
0383 int AnaMvtxTelescopeHits::End(PHCompositeNode * topNode)
0384 {
0385   _f->Write();
0386   _f->Close();
0387
0388   return 0;
0389 }
0390
0391 /*
0392  * Get all the necessary nodes from the node tree
0393  */
0394 int AnaMvtxTelescopeHits::GetNodes(PHCompositeNode * topNode)
0395 {
0396
0397   // Input Svtx Clusters
0398   _clustermap = findNode::getClass<SvtxClusterMap>(topNode, "SvtxClusterMap");
0399   if (!_clustermap && _ievent < 2) {
0400     std::cout << PHWHERE << " SvtxClusterMap node not found on node tree"
0401               << std::endl;
0402     return Fun4AllReturnCodes::ABORTEVENT;
0403   }
0404
0405   return Fun4AllReturnCodes::EVENT_OK;
0406
0407 }
0408
0409 /*
0410  * Simple helper function for calculating the slope between two points
0411  */
0412 double
0413 AnaMvtxTelescopeHits::CalcSlope(double x0, double y0, double x1, double y1)
0414 {
0415   return (y1 - y0) / (x1 - x0);
0416 }
0417
0418 /*
0419  * Simple helper function for calculating intercept
0420  */
0421 double
0422 AnaMvtxTelescopeHits::CalcIntecept(double x0, double y0, double m)
0423 {
0424   return y0 - x0 * m;
0425 }
0426
0427 /*
0428  * Simple helper function for calculating projection
0429  */
0430 double
0431 AnaMvtxTelescopeHits::CalcProjection(double x, double m, double b)
0432 {
0433   return m * x + b;
0434 }
0435
0436
0437 // /*
0438 //  *
0439 //  */
0440 // int AnaMvtxTelescopeHits::GetTrackCandidates(LyrClusMap* clusmap,
0441 //     TrkVec* trkvec)
0442 // {
0443
0444 //   //-- Make a vector to hold all track candidates for a given first cluster
0445 //   TrkVec trkcand;
0446
0447 //   //-- Loop over all clusters in the first layer of the cluster map
0448 //   unsigned int lyr = (clusmap->begin())->first;
0449 //   for (LyrClusMap::iterator itr = clusmap->lower_bound(lyr);
0450 //        itr != clusmap->upper_bound(lyr);
0451 //        ++itr)
0452 //   {
0453 //     trkcand.clear();
0454
0455 //     //-- Find all track candidates for this cluster
0456 //     ClusVec trk;
0457 //     trk.insert(itr->second);
0458 //     LinkClusters(&trk, lyr+1, clusmap, &trkcand);
0459
0460 //     //-- Choose the best track candidate using this cluster
0461 //     ClusVec best_trk = ChooseBestTrk(&trkcand);
0462 //     if ( best_trk.size() > 2 )
0463 //       trkvec->insert(best_trk);
0464 //   }
0465
0466 //   return trkvec->size();
0467 // }
0468
0469 // /*
0470 //  *
0471 //  */
0472 // void AnaMvtxTelescopeHits::LinkClusters(ClusVec* trk,
0473 //                                         unsigned int next_lyr,
0474 //                                         LyrClusMap* clusmap,
0475 //                                         TrkVec* trkvec)
0476 // {
0477
0478 //   //-- Check layer
0479 //   if ( next_lyr > 3 )
0480 //   {
0481 //     // add this trk to the trk vector, we're done with this trk
0482 //     // we only consider tracks with at least 3 clusters
0483 //     if ( trk->size() > 2 )
0484 //       trkvec->push_back(*trk);
0485 //     return;
0486 //   }
0487
0488 //   bool found = false;
0489
0490 //   //-- Search for clusters in the next layer
0491 //   for ( LyrClusMap::iterator itr = LyrClusMap.lower_bound(next_lyr);
0492 //         itr != LyrClusMap.upper_bound(next_lyr);
0493 //         ++itr)
0494 //   {
0495 //     // make a copy of this trk candidate and run the recursion
0496 //     ClusVec new_trk = *trk;
0497 //     new_trk.insert(itr->second);
0498 //     LinkClusters(&new_trk, next_lyr + 1, clusmap, trkvec);
0499 //     found = true;
0500 //   }
0501
0502 //   //-- if no cluster found in next layer, try next layer
0503 //   if ( !found )
0504 //   {
0505 //     LinkClusters(*trk, next_lyr + 2, clusmap, trkvec);
0506 //   }
0507
0508 // }
0509
0510 // /*
0511 //  * Choose the best track candidates which all start with the same cluster
0512 //  */
0513 // ClusVec AnaMvtxTelescopeHits::ChooseBestTrk(TrkVec* trkcnd)
0514 // {
0515
0516 //   int best_idx = -1;
0517 //   double best_chi2 = 100;
0518 //   for (TrkVec::iterator itr = trkcnd->begin();
0519 //        itr != trkcnd->end();
0520 //        ++itr)
0521 //   {
0522 //     double chi2 = FitTrk(&itr->second);
0523 //     if ( chi2 < best_chi2 )
0524 //     {
0525 //       best_chi2 = chi2;
0526 //       best_idx = itr->first;
0527 //     }
0528 //   }
0529
0530 //   if ( best_idx >= 0 && best_chi2 < 5 )
0531 //   {
0532 //     return trkcnd->at(best_idx);
0533 //   }
0534 //   else
0535 //   {
0536 //     return ClusVec;
0537 //   }
0538 // }
0539
0540 // /*
0541 //  *
0542 //  */
0543 // double AnaMvtxTelescopeHits::FitTrk(ClusVec* trk)
0544 // {
0545 //   return 0;
0546 // }
0547
0548 /*
0549  * Simple generalized linear least-squares fitter.
0550  * Solve y(X) = beta'*X + error(L) for beta (vector of regression coefs.)
0551  * L is the inverse covariance matrix for y.
0552  * Least-squares solution involves solving X' L X beta = X' L y
0553  *
0554  */
0555 /*
0556 TVectorD
0557 AnaMvtxTelescopeHits::SolveGLS(TMatrixD &X, TVectorD &y, TMatrixD &L, TMatrixD &cov)
0558 {
0559
0560   double tol = 1.e-15;
0561   TMatrixD XT(X); XT.T();
0562   TMatrixD A = XT * L * X;
0563   TVectorD b = XT * L * y;
0564
0565   // Now solve A*beta = b using SVD. Decompose A = U S V'.
0566   TDecompSVD svd(A);
0567   TMatrixD UT = svd.GetU(); UT.T();
0568   TMatrixD V  = svd.GetV();
0569
0570   // Construct Moore-Penrose pseudoinverse of S
0571   TVectorD s = svd.GetSig();
0572   int n = s.GetNrows();
0573   TMatrixD Sd(n, n);
0574   for (int i = 0; i < n; i++)
0575     Sd(i, i) = s(i) > tol ? 1. / s(i) : 0.;
0576
0577   // Result
0578   TVectorD beta = V * Sd * UT * b;
0579
0580   // and covariance matrix
0581   V *= Sd;
0582   TMatrixD C(V, TMatrixD::kMultTranspose, V);
0583   cov.ResizeTo(C);
0584   cov = C;
0585
0586   return beta;
0587 }
0588 */
0589