PHGenFitPkg/Example/testPHGenFit.cc

0001 /*!
0002  *  \file       testPHGenFit.cc
0003  *  \brief      Program to demonstrate the usage of PHGenFit.
0004  *  \author     Haiwang Yu <yuhw@nmsu.edu>
0005  */
0006
0007 // STL
0008 #include <vector>
0009
0010 // BOOST
0011 #include <boost/make_shared.hpp>
0012
0013 #define SMART(expr) boost::shared_ptr<expr>
0014 #define NEW(expr) boost::make_shared<expr>
0015
0016 #include <phfield/PHFieldUtility.h>
0017
0018 // ROOT
0019 #include <TCanvas.h>
0020 #include <TF1.h>
0021 #include <TFile.h>
0022 #include <TH1D.h>
0023 #include <TH2D.h>
0024 #include <TMath.h>
0025 #include <TMatrixDSym.h>
0026 #include <TProfile.h>
0027 #include <TROOT.h>
0028 #include <TRandom.h>
0029 #include <TStyle.h>
0030 #include <TTree.h>
0031 #include <TVector3.h>
0032
0033 // GenFit
0034 #include <GenFit/AbsTrackRep.h>
0035 #include <GenFit/RKTrackRep.h>
0036 #include <GenFit/StateOnPlane.h>
0037
0038 // PHGenFit
0039 #include <phgenfit/Fitter.h>
0040 #include <phgenfit/Measurement.h>
0041 #include <phgenfit/PlanarMeasurement.h>
0042 #include <phgenfit/Track.h>
0043
0044 #define LogDEBUG std::cout << "DEBUG: " << __LINE__ << "\n"
0045
0046 // void pause() {
0047 //   std::cout << "Press ENTER to continue..." << std::flush;
0048 //   std::cin.clear();  // use only if a human is involved
0049 //   std::cin.flush();  // use only if a human is involved
0050 //   std::cin.ignore( std::numeric_limits<std::streamsize>::max(), '\n' );
0051 // }
0052
0053 int main(int /*argc*/, char ** /*argv*/)
0054 {
0055   TFile *fPHG4Hits = TFile::Open("AnaSvtxTracksForGenFit.root", "read");
0056   if (!fPHG4Hits)
0057   {
0058     std::cout << "No TFile Openned: " << __LINE__ << "\n";
0059     return -1;
0060   }
0061   TTree *T = (TTree *) fPHG4Hits->Get("tracks");
0062   if (!T)
0063   {
0064     std::cout << "No TTree Found: " << __LINE__ << "\n";
0065     return -1;
0066   }
0067
0068   //! Initiallize Geometry, Field, Fitter
0069   PHGenFit::Fitter *fitter = new PHGenFit::Fitter("sPHENIX_Geo.root", PHFieldUtility::BuildFieldMap(PHFieldUtility::DefaultFieldConfig(), 1));
0070
0071   double resolution_detector_xy = 0.005 / 3.;  // 50/3. micron
0072
0073   TH2D *hpT_residual_vs_pT = new TH2D("hpT_residual_vs_pT", "#Delta pT/pT; pT[GeV/c]; #Delta pT/pT", 40, 0.5, 40.5, 1000, -1, 1);
0074
0075   TH2D *hDCAr_vs_pT = new TH2D("hDCAr_vs_pT", "DCAr vs. p; p [GeV/c]; DCAr [cm]", 40, 0.5, 40.5, 1000, -0.1, 0.1);
0076
0077 #define NLAYERS 7
0078
0079   Float_t Cluster_x[NLAYERS];
0080   Float_t Cluster_y[NLAYERS];
0081   Float_t Cluster_z[NLAYERS];
0082   Float_t Cluster_size_dphi[NLAYERS];
0083   Float_t Cluster_size_dz[NLAYERS];
0084   Float_t True_px;
0085   Float_t True_py;
0086   Float_t True_pz;
0087   Float_t True_vx;
0088   Float_t True_vy;
0089   Float_t True_vz;
0090   Float_t AlanDion_px;
0091   Float_t AlanDion_py;
0092   Float_t AlanDion_pz;
0093   Float_t AlanDion_dca2d;
0094   Int_t nhits;
0095
0096   T->SetBranchAddress("nhits", &nhits);
0097   T->SetBranchAddress("gpx", &True_px);
0098   T->SetBranchAddress("gpy", &True_py);
0099   T->SetBranchAddress("gpz", &True_pz);
0100   T->SetBranchAddress("gvx", &True_vx);
0101   T->SetBranchAddress("gvy", &True_vy);
0102   T->SetBranchAddress("gvz", &True_vz);
0103   T->SetBranchAddress("px", &AlanDion_px);
0104   T->SetBranchAddress("py", &AlanDion_py);
0105   T->SetBranchAddress("pz", &AlanDion_pz);
0106   T->SetBranchAddress("dca2d", &AlanDion_dca2d);
0107   T->SetBranchAddress("x", Cluster_x);
0108   T->SetBranchAddress("y", Cluster_y);
0109   T->SetBranchAddress("z", Cluster_z);
0110   T->SetBranchAddress("size_dphi", Cluster_size_dphi);
0111   T->SetBranchAddress("size_dz", Cluster_size_dz);
0112
0113   double nentries = 10;
0114   // double nentries = T->GetEntries();
0115   for (unsigned int ientry = 0; ientry < nentries; ++ientry)
0116   {
0117     // T->GetEntry(atoi(argv[1]));
0118     if (ientry % 1000 == 0)
0119     {
0120       std::cout << "Processing: " << 100. * ientry / nentries << "%"
0121                 << "\n";
0122     }
0123
0124     T->GetEntry(ientry);
0125
0126     if (nhits < 0)
0127     {
0128       // LogDEBUG;
0129       continue;
0130     }
0131
0132     // true start values
0133     TVector3 init_pos(0, 0, 0);  // cm
0134     TVector3 True_mom(True_px, True_py, True_pz);
0135
0136     // Seed: use smeared values
0137     const bool smearPosMom = true;                        // init the Reps with smeared init_pos and True_mom
0138     const double posSmear = 10 * resolution_detector_xy;  // cm
0139     const double momSmear = 3. / 180. * TMath::Pi();      // rad
0140     const double momMagSmear = 0.1;                       // relative
0141
0142     TVector3 seed_pos(init_pos);
0143     TVector3 seed_mom(True_mom);
0144     if (smearPosMom)
0145     {
0146       seed_pos.SetX(gRandom->Gaus(seed_pos.X(), posSmear));
0147       seed_pos.SetY(gRandom->Gaus(seed_pos.Y(), posSmear));
0148       seed_pos.SetZ(gRandom->Gaus(seed_pos.Z(), posSmear));
0149
0150       seed_mom.SetPhi(gRandom->Gaus(True_mom.Phi(), momSmear));
0151       seed_mom.SetTheta(gRandom->Gaus(True_mom.Theta(), momSmear));
0152       seed_mom.SetMag(
0153           gRandom->Gaus(True_mom.Mag(),
0154                         momMagSmear * True_mom.Mag()));
0155     }
0156
0157     // approximate covariance
0158     TMatrixDSym seed_cov(6);
0159
0160     for (int idim = 0; idim < 3; ++idim)
0161     {
0162       seed_cov(idim, idim) = resolution_detector_xy * resolution_detector_xy;
0163     }
0164     for (int idim = 3; idim < 6; ++idim)
0165     {
0166       seed_cov(idim, idim) = pow(
0167           resolution_detector_xy / NLAYERS / sqrt(3), 2);
0168     }
0169
0170     //! Build TrackRep from particle assumption
0171     int pid = -13;  // mu+
0172     // SMART(genfit::AbsTrackRep) rep = NEW(genfit::RKTrackRep)(pid);
0173     genfit::AbsTrackRep *rep = new genfit::RKTrackRep(pid);
0174
0175     //! Initiallize track with seed from pattern recognition
0176     PHGenFit::Track *track = new PHGenFit::Track(rep, seed_pos,
0177                                                  seed_mom, seed_cov);
0178
0179     //! Create measurements
0180     std::vector<PHGenFit::Measurement *> measurements;
0181     for (int imeas = 0; imeas < NLAYERS; imeas++)
0182     {
0183       TVector3 pos(Cluster_x[imeas], Cluster_y[imeas], Cluster_z[imeas]);
0184       TVector3 n(Cluster_x[imeas], Cluster_y[imeas], 0);
0185       TVector3 v(0, 0, 1);
0186       TVector3 u = v.Cross(n);
0187
0188       PHGenFit::Measurement *meas = new PHGenFit::PlanarMeasurement(
0189           pos, u, v,
0190           Cluster_size_dphi[imeas], Cluster_size_dz[imeas]);
0191       measurements.push_back(meas);
0192     }
0193
0194     //! Add measurements to track
0195     track->addMeasurements(measurements);
0196
0197     //! Fit the track
0198     fitter->processTrack(track, false);
0199
0200     //!
0201     genfit::MeasuredStateOnPlane *state_at_beam_line = track->extrapolateToLine(
0202         TVector3(0, 0, 0), TVector3(0, 0, 1));
0203     // state_at_beam_line->Print();
0204
0205     //      genfit::MeasuredStateOnPlane* state_at_layer_6 = track->extrapolateToCylinder(80.,
0206     //                      TVector3(0, 0, 0), TVector3(0, 0, 1));
0207     //      //state_at_layer_6->Print();
0208     //      delete state_at_layer_6;
0209
0210     TVector3 GenFit_mom = state_at_beam_line->getMom();
0211
0212     TVector3 AlanDion_mom(AlanDion_px, AlanDion_py, AlanDion_pz);
0213
0214     hpT_residual_vs_pT->Fill(True_mom.Pt(), (GenFit_mom.Pt() - True_mom.Pt()) / True_mom.Pt());
0215
0216     hDCAr_vs_pT->Fill(True_mom.Mag(), state_at_beam_line->getState()[3]);
0217
0218     delete state_at_beam_line;
0219     delete track;
0220     measurements.clear();
0221   }
0222
0223   gStyle->SetOptFit();
0224   gStyle->SetOptStat(000000000);
0225
0226   TF1 *tf_pT_resolution = new TF1("tf_pT_resolution", "sqrt([0]*[0] + x*x*[1]*[1])", 0, 40);
0227   tf_pT_resolution->SetParameters(0, 0);
0228   TCanvas *c3 = new TCanvas("c3", "c3");
0229   c3->Divide(2, 1);
0230   c3->cd(1);
0231   hpT_residual_vs_pT->FitSlicesY();
0232   TH1D *hpT_resolution_vs_pT = (TH1D *) gDirectory->Get("hpT_residual_vs_pT_2");
0233   hpT_resolution_vs_pT->SetTitle("PHGenFit: #sigma_{p_{T}}/p_{T}; p_{T}[GeV/c]; #sigma_{p_{T}}/p_{T}");
0234   hpT_resolution_vs_pT->SetMarkerStyle(20);
0235   hpT_resolution_vs_pT->Draw("e");
0236   hpT_resolution_vs_pT->Fit(tf_pT_resolution);
0237   c3->cd(2);
0238   hDCAr_vs_pT->FitSlicesY();
0239   TH1D *hDCAr_resolution_vs_pT = (TH1D *) gDirectory->Get("hDCAr_vs_pT_2");
0240   hDCAr_resolution_vs_pT->SetTitle(
0241       "PHGenFit: #sigma_{DCAr} [cm]; p [GeV/c]; #sigma_{DCAr}");
0242   hDCAr_resolution_vs_pT->SetMarkerStyle(20);
0243   hDCAr_resolution_vs_pT->Draw("e");
0244   // hDCAr_resolution_vs_pT->Fit(tf_pT_resolution);
0245   c3->Print("pT_DCA_resolution.root");
0246
0247   //! Event display
0248   // fitter->displayEvent();
0249
0250   fPHG4Hits->Close();
0251
0252   delete fitter;
0253
0254   // pause();
0255
0256   std::cout << "SUCCESS! \n";
0257
0258   return 0;
0259 }