packages/MvtxPrototype1/MvtxQAHisto.C

0001 ////////////////////////////////////////////////////////////////////////////////
0002 //
0003 // This module is desgined to grab svtx tracks and put truth and cluster
0004 // information into a TTree for GenFit testing
0005 //
0006 ////////////////////////////////////////////////////////////////////////////////
0007 //
0008 // Darren McGlinchey
0009 // 1 Apr 2016
0010 //
0011 ////////////////////////////////////////////////////////////////////////////////
0012
0013
0014 #include "MvtxQAHisto.h"
0015
0016 #include <phool/phool.h>
0017 #include <phool/getClass.h>
0018 #include <g4main/PHG4HitContainer.h>
0019 #include <g4main/PHG4TruthInfoContainer.h>
0020 #include <g4main/PHG4Particle.h>
0021 #include <g4main/PHG4Hit.h>
0022 #include <g4main/PHG4VtxPoint.h>
0023 #include <fun4all/PHTFileServer.h>
0024 #include <fun4all/Fun4AllServer.h>
0025 #include <phool/recoConsts.h>
0026
0027 #include <trackbase/TrkrDefUtil.h>
0028 #include <trackbase/TrkrClusterContainer.h>
0029 #include <trackbase/TrkrCluster.h>
0030 #include <trackbase/TrkrClusterv1.h>
0031 #include <trackbase/TrkrHitSetContainer.h>
0032 #include <trackbase/TrkrHitSet.h>
0033 #include <trackbase/TrkrHitSetv1.h>
0034 #include <mvtx/MvtxDefUtil.h>
0035 #include <mvtx/MvtxHitSetv1.h>
0036
0037 #include <TTree.h>
0038 #include <TH1F.h>
0039 #include <TH2F.h>
0040 #include <TVector3.h>
0041
0042 #include <iostream>
0043
0044 using namespace std;
0045
0046 //----------------------------------------------------------------------------//
0047 //-- Constructor:
0048 //--  simple initialization
0049 //----------------------------------------------------------------------------//
0050 MvtxQAHisto::MvtxQAHisto(const string &name):
0051   SubsysReco( name )
0052 {
0053   //initialize
0054   _outfile = "MvtxQAHisto.root";
0055
0056     beam_x[0] = 826.693;
0057     beam_x[1] = 822.267;
0058     beam_x[2] = 818.413;
0059     beam_x[3] = 830.190;
0060
0061     beam_y[0] = 158.773;
0062     beam_y[1] = 167.537;
0063     beam_y[2] = 181.318;
0064     beam_y[3] = 190.988;
0065
0066 }
0067
0068 //----------------------------------------------------------------------------//
0069 //-- Init():
0070 //--   Intialize all histograms, trees, and ntuples
0071 //----------------------------------------------------------------------------//
0072 int MvtxQAHisto::Init(PHCompositeNode *topNode)
0073 {
0074
0075     //recoConsts *rc = recoConsts::instance();
0076   cout << PHWHERE << " Openning file " << _outfile << endl;
0077   PHTFileServer::get().open( _outfile, "RECREATE");
0078
0079     for (int ichip=0; ichip<4; ichip++){
0080         h2d_hit[ichip] = new TH2F(Form("h2d_hit_chip%d",ichip),"",1024,-0.5,1023.5,512,-0.5,511.5);
0081         h1d_hit_per_evt[ichip] = new TH1F(Form("h1d_hit_per_evt_chip%d",ichip),"",101,-0.5,100);
0082         h1d_clus_per_evt[ichip] = new TH1F(Form("h1d_clus_per_evt_chip%d",ichip),"",51,-0.5,50.5);
0083
0084         h2d_hit_beam[ichip] = new TH2F(Form("h2d_hit_beam_chip%d",ichip),"",1025,-512.5,512.5,513,-256.5,256.5);
0085         h2d_hit_trk[ichip] = new TH2F(Form("h2d_hit_trk_chip%d",ichip),"",1025,-512.5,512.5,513,-256.5,256.5);
0086         h2d_clus[ichip] = new TH2F(Form("h2d_clus_chip%d",ichip),"",1024,0,3.0,512,0,1.5);
0087         h2d_clus_beam[ichip] = new TH2F(Form("h2d_clus_beam_chip%d",ichip),"",1024,-3.0/2,3.0/2,512,-1.5/2,1.5/2);
0088
0089         h1d_clus_size_x[ichip] = new TH1F(Form("h1d_clus_size_x_chip%d",ichip),"",51,-0.5,50.5);
0090         h1d_clus_size_z[ichip] = new TH1F(Form("h1d_clus_size_z_chip%d",ichip),"",51,-0.5,50.5);
0091
0092     }
0093     h1d_trk_finder_x = new TH1F("h1d_trk_finder_x","",513,-256.5,256.5);
0094     h1d_trk_finder_z = new TH1F("h1d_trk_finder_z","",1025,-512.5,512.5);
0095     h2d_trk_finder = new TH2F("h2d_trk_finder","",1025,-512.5,512.5,513,-256.5,256.5);
0096
0097     h1d_clus_associated = new TH1F("h1d_clus_associated","",11,-0.5,10.5);
0098     h1d_clus_eff = new TH1F("h1d_clus_eff","",5,-0.5,4.5);
0099
0100   return 0;
0101 }
0102
0103 //----------------------------------------------------------------------------//
0104 //-- process_event():
0105 //--   Call user instructions for every event.
0106 //--   This function contains the analysis structure.
0107 //----------------------------------------------------------------------------//
0108 int MvtxQAHisto::process_event(PHCompositeNode *topNode)
0109 {
0110   _event++;
0111   //if (1)
0112   if (_event % 1000 == 0)
0113     cout << PHWHERE << "Events processed: " << _event << endl;
0114
0115   GetNodes(topNode);
0116
0117     TrkrDefUtil trkrdefutil;
0118     MvtxDefUtil mvtxdefutil;
0119
0120     int nhit_per_chip[4] = {0};
0121     int nclus_per_chip[4] = {0};
0122
0123     double avg_x = 0.0;
0124     double avg_y = 0.0;
0125     double norm = 0.0;
0126
0127     h1d_trk_finder_x->Reset();
0128     h1d_trk_finder_z->Reset();
0129     h2d_trk_finder->Reset();
0130
0131     TrkrHitSetContainer::ConstRange iter_range = hitsetcon->GetHitSets();
0132     for ( TrkrHitSetContainer::ConstIterator iter = iter_range.first; iter!=iter_range.second; ++iter){
0133
0134         int ichip = int(mvtxdefutil.GetChipId(iter->first));
0135
0136         MvtxHitSetv1 *hitset = static_cast<MvtxHitSetv1 *>(iter->second);
0137
0138         MvtxHitSetv1::ConstRange hit_iter_range = hitset->GetHits();
0139         for ( MvtxHitSetv1::ConstIterator hit_iter = hit_iter_range.first; hit_iter!=hit_iter_range.second; ++hit_iter){
0140
0141             int icol = int(hit_iter->first);
0142             int irow = int(hit_iter->second);
0143
0144             h2d_hit[ichip]->Fill(icol,irow);
0145
0146             h1d_trk_finder_z->Fill(icol-beam_x[ichip]);
0147             h1d_trk_finder_x->Fill(irow-beam_y[ichip]);
0148             h2d_trk_finder->Fill(icol-beam_x[ichip],irow-beam_y[ichip]);
0149
0150             nhit_per_chip[ichip]++;
0151
0152         }//hit_iter
0153     }//iter
0154
0155     TrkrClusterContainer::ConstRange clus_iter_range = cluscon->GetClusters();
0156     for (TrkrClusterContainer::ConstIterator clus_iter = clus_iter_range.first; clus_iter!=clus_iter_range.second; ++clus_iter){
0157
0158         int ichip  = trkrdefutil.GetLayer(clus_iter->first);
0159
0160         TrkrCluster *clus = clus_iter->second;
0161
0162         h2d_clus[ichip]->Fill(clus->GetZ(),clus->GetX());
0163         h2d_clus_beam[ichip]->Fill(clus->GetZ() - beam_x[ichip]*(28e-4),clus->GetX() - beam_y[ichip]*(28e-4));
0164
0165         h1d_clus_size_z[ichip]->Fill(clus->GetZSize());
0166         h1d_clus_size_x[ichip]->Fill(clus->GetPhiSize());
0167
0168         nclus_per_chip[ichip]++;
0169     }//
0170
0171     int nchip_w_clus = 0;
0172
0173     for (int ichip=0; ichip<4; ichip++){
0174         h1d_hit_per_evt[ichip]->Fill(nhit_per_chip[ichip]);
0175         h1d_clus_per_evt[ichip]->Fill(nclus_per_chip[ichip]);
0176
0177         if ( nclus_per_chip[ichip]>0 ) nchip_w_clus++;
0178     }
0179
0180     //if ( nchip_w_clus<3 ) return 0;
0181     //if ( 1 ) return 0;
0182
0183     //fast track reconstruction
0184     /*
0185     double trk_z = h1d_trk_finder_z->GetBinCenter(h1d_trk_finder_z->GetMaximumBin());
0186     double trk_x = h1d_trk_finder_x->GetBinCenter(h1d_trk_finder_x->GetMaximumBin());
0187     double trk_z_w = h1d_trk_finder_z->GetBinContent(h1d_trk_finder_z->GetMaximumBin());
0188     double trk_x_w = h1d_trk_finder_x->GetBinContent(h1d_trk_finder_x->GetMaximumBin());
0189     */
0190
0191     /*
0192     double trk_z = 0;
0193     double trk_x = 0;
0194     double trk_w = 0.1;
0195     for (int ix=0; ix<h2d_trk_finder->GetNbinsX(); ix++){
0196         for (int iy=0; iy<h2d_trk_finder->GetNbinsY(); iy++){
0197             if ( h2d_trk_finder->GetBinContent(ix+1,iy+1)<trk_w ) continue;
0198
0199             trk_z = h1d_trk_finder_z->GetBinCenter(ix+1);
0200             trk_x = h1d_trk_finder_x->GetBinCenter(iy+1);
0201             trk_w = h2d_trk_finder->GetBinContent(ix+1,iy+1);
0202
0203         }
0204     }
0205
0206     if ( trk_w<3 ) return 0;
0207
0208     int nclus_associated = 0;
0209     int nclus_associated_per_chip[4] = {0};
0210
0211     //TrkrClusterContainer::ConstRange clus_iter_range = cluscon->GetClusters();
0212     for (TrkrClusterContainer::ConstIterator clus_iter = clus_iter_range.first; clus_iter!=clus_iter_range.second; ++clus_iter){
0213
0214         int ichip  = trkrdefutil.GetLayer(clus_iter->first);
0215
0216         TrkrCluster *clus = clus_iter->second;
0217
0218         double clusz = clus->GetZ()/(28e-4) - beam_x[ichip];
0219         double clusx = clus->GetX()/(28e-4) - beam_y[ichip];
0220
0221         if ( fabs(clusz-trk_z)<5 && fabs(clusx-trk_x)<5 ){
0222             nclus_associated_per_chip[ichip]++;
0223             nclus_associated++;
0224         }
0225
0226     }//
0227
0228     h1d_clus_associated->Fill(nclus_associated);
0229
0230     if ( nclus_associated<3 ) return 0;
0231
0232     for (int ichip=0; ichip<4; ichip++){
0233         if ( nclus_associated_per_chip[ichip]>0 ) h1d_clus_eff->Fill(ichip);
0234     }
0235
0236     //cout << "N Associated Clus: " << nclus_associated << endl;
0237
0238
0239     //cout << "Z: " << trk_z << ", X: " << trk_x << ", S: " << trk_w << endl;
0240
0241     //TrkrHitSetContainer::ConstRange iter_range = hitsetcon->GetHitSets();
0242     for ( TrkrHitSetContainer::ConstIterator iter = iter_range.first; iter!=iter_range.second; ++iter){
0243
0244         int ichip = int(mvtxdefutil.GetChipId(iter->first));
0245
0246         MvtxHitSetv1 *hitset = static_cast<MvtxHitSetv1 *>(iter->second);
0247
0248         MvtxHitSetv1::ConstRange hit_iter_range = hitset->GetHits();
0249         for ( MvtxHitSetv1::ConstIterator hit_iter = hit_iter_range.first; hit_iter!=hit_iter_range.second; ++hit_iter){
0250
0251             int icol = int(hit_iter->first);
0252             int irow = int(hit_iter->second);
0253
0254             h2d_hit_trk[ichip]->Fill(icol-beam_x[ichip]-trk_z,irow-beam_y[ichip]-trk_x);
0255
0256         }//hit_iter
0257     }//iter
0258     */
0259
0260
0261   return 0;
0262 }
0263
0264 //----------------------------------------------------------------------------//
0265 //-- End():
0266 //--   End method, wrap everything up
0267 //----------------------------------------------------------------------------//
0268 int MvtxQAHisto::End(PHCompositeNode *topNode)
0269 {
0270
0271     cout << "-----MvtxQAHisto::End------" << endl;
0272
0273   PHTFileServer::get().cd( _outfile );
0274     PHTFileServer::get().write( _outfile );
0275   //PHTFileServer::get().close();
0276
0277   return 0;
0278 }
0279
0280
0281 //----------------------------------------------------------------------------//
0282 //-- GetNodes():
0283 //--   Get all the all the required nodes off the node tree
0284 //----------------------------------------------------------------------------//
0285 void MvtxQAHisto::GetNodes(PHCompositeNode * topNode)
0286 {
0287   //DST objects
0288     //
0289     //
0290     cluscon = findNode::getClass<TrkrClusterContainer>(topNode, "TrkrClusterContainer");
0291     if (cluscon==0)
0292     {
0293         cout << "MvtxQAHisto::Process_Event - TrkrClusterContainer not found" << endl;
0294         exit(-1);
0295     }
0296
0297     hitsetcon = findNode::getClass<TrkrHitSetContainer>(topNode, "TrkrHitSetContainer");
0298     if (hitsetcon==0)
0299     {
0300         cout << "MvtxQAHisto::Process_Event - TrkrHitSetContainer not found" << endl;
0301         exit(-1);
0302     }
0303
0304 }
0305
0306
0307