Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
testPHGenFit.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file testPHGenFit.cc
1 
7 //STL
8 #include <vector>
9 
10 //BOOST
11 #include<boost/make_shared.hpp>
12 
13 #define SMART(expr) boost::shared_ptr<expr>
14 #define NEW(expr) boost::make_shared<expr>
15 
16 #include <phfield/PHFieldUtility.h>
17 
18 //ROOT
19 #include <TVector3.h>
20 #include <TMatrixDSym.h>
21 #include <TF1.h>
22 #include <TH1D.h>
23 #include <TH2D.h>
24 #include <TProfile.h>
25 #include <TFile.h>
26 #include <TTree.h>
27 #include <TCanvas.h>
28 #include <TROOT.h>
29 #include <TStyle.h>
30 #include <TMath.h>
31 #include <TRandom.h>
32 
33 //GenFit
34 #include <GenFit/AbsTrackRep.h>
35 #include <GenFit/RKTrackRep.h>
36 #include <GenFit/StateOnPlane.h>
37 
38 //PHGenFit
39 #include <phgenfit/Fitter.h>
40 #include <phgenfit/Track.h>
41 #include <phgenfit/Measurement.h>
42 #include <phgenfit/PlanarMeasurement.h>
43 
44 #define LogDEBUG std::cout<<"DEBUG: "<<__LINE__<<"\n"
45 
46 //void pause() {
47 // std::cout << "Press ENTER to continue..." << std::flush;
48 // std::cin.clear(); // use only if a human is involved
49 // std::cin.flush(); // use only if a human is involved
50 // std::cin.ignore( std::numeric_limits<std::streamsize>::max(), '\n' );
51 //}
52 
53 int main(int argc, char**argv) {
54 
55  TFile *fPHG4Hits = TFile::Open("AnaSvtxTracksForGenFit.root", "read");
56  if (!fPHG4Hits) {
57  std::cout << "No TFile Openned: " << __LINE__ << "\n";
58  return -1;
59  }
60  TTree *T = (TTree*) fPHG4Hits->Get("tracks");
61  if (!T) {
62  std::cout << "No TTree Found: " << __LINE__ << "\n";
63  return -1;
64  }
65 
67  PHGenFit::Fitter* fitter = new PHGenFit::Fitter("sPHENIX_Geo.root",PHFieldUtility::BuildFieldMap(PHFieldUtility::DefaultFieldConfig(), 1));
68 
69 
70  double resolution_detector_xy = 0.005/3.; //50/3. micron
71 
72 
73  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);
74 
75  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);
76 
77 
78 #define NLAYERS 7
79 
80 
81  Float_t Cluster_x[NLAYERS];
82  Float_t Cluster_y[NLAYERS];
83  Float_t Cluster_z[NLAYERS];
84  Float_t Cluster_size_dphi[NLAYERS];
85  Float_t Cluster_size_dz[NLAYERS];
86  Float_t True_px;
87  Float_t True_py;
88  Float_t True_pz;
89  Float_t True_vx;
90  Float_t True_vy;
91  Float_t True_vz;
92  Float_t AlanDion_px;
93  Float_t AlanDion_py;
94  Float_t AlanDion_pz;
95  Float_t AlanDion_dca2d;
96  Int_t nhits;
97 
98  T->SetBranchAddress("nhits", &nhits);
99  T->SetBranchAddress("gpx", &True_px);
100  T->SetBranchAddress("gpy", &True_py);
101  T->SetBranchAddress("gpz", &True_pz);
102  T->SetBranchAddress("gvx", &True_vx);
103  T->SetBranchAddress("gvy", &True_vy);
104  T->SetBranchAddress("gvz", &True_vz);
105  T->SetBranchAddress("px", &AlanDion_px);
106  T->SetBranchAddress("py", &AlanDion_py);
107  T->SetBranchAddress("pz", &AlanDion_pz);
108  T->SetBranchAddress("dca2d", &AlanDion_dca2d);
109  T->SetBranchAddress("x", Cluster_x);
110  T->SetBranchAddress("y", Cluster_y);
111  T->SetBranchAddress("z", Cluster_z);
112  T->SetBranchAddress("size_dphi", Cluster_size_dphi);
113  T->SetBranchAddress("size_dz", Cluster_size_dz);
114 
115  double nentries = 10;
116  //double nentries = T->GetEntries();
117  for (unsigned int ientry = 0; ientry < nentries; ++ientry) {
118  //T->GetEntry(atoi(argv[1]));
119  if(ientry%1000==0) std::cout<<"Processing: "<<100.*ientry/nentries <<"%"<<"\n";
120 
121  T->GetEntry(ientry);
122 
123  if (nhits < 0) {
124  //LogDEBUG;
125  continue;
126  }
127 
128  // true start values
129  TVector3 init_pos(0, 0, 0); //cm
130  TVector3 True_mom(True_px, True_py, True_pz);
131 
132  // Seed: use smeared values
133  const bool smearPosMom = true; // init the Reps with smeared init_pos and True_mom
134  const double posSmear = 10 * resolution_detector_xy; // cm
135  const double momSmear = 3. / 180. * TMath::Pi(); // rad
136  const double momMagSmear = 0.1; // relative
137 
138  TVector3 seed_pos(init_pos);
139  TVector3 seed_mom(True_mom);
140  if (smearPosMom) {
141  seed_pos.SetX(gRandom->Gaus(seed_pos.X(), posSmear));
142  seed_pos.SetY(gRandom->Gaus(seed_pos.Y(), posSmear));
143  seed_pos.SetZ(gRandom->Gaus(seed_pos.Z(), posSmear));
144 
145  seed_mom.SetPhi(gRandom->Gaus(True_mom.Phi(), momSmear));
146  seed_mom.SetTheta(gRandom->Gaus(True_mom.Theta(), momSmear));
147  seed_mom.SetMag(
148  gRandom->Gaus(True_mom.Mag(),
149  momMagSmear * True_mom.Mag()));
150  }
151 
152  // approximate covariance
153  TMatrixDSym seed_cov(6);
154 
155  for (int idim = 0; idim < 3; ++idim)
156  seed_cov(idim, idim) = resolution_detector_xy
157  * resolution_detector_xy;
158  for (int idim = 3; idim < 6; ++idim)
159  seed_cov(idim, idim) = pow(
160  resolution_detector_xy / NLAYERS / sqrt(3), 2);
161 
162 
164  int pid = -13; //mu+
165  //SMART(genfit::AbsTrackRep) rep = NEW(genfit::RKTrackRep)(pid);
166  genfit::AbsTrackRep* rep = new genfit::RKTrackRep(pid);
167 
169  PHGenFit::Track* track = new PHGenFit::Track(rep, seed_pos,
170  seed_mom, seed_cov);
171 
173  std::vector<PHGenFit::Measurement*> measurements;
174  for (int imeas = 0; imeas < NLAYERS; imeas++) {
175  TVector3 pos(Cluster_x[imeas],Cluster_y[imeas],Cluster_z[imeas]);
176  TVector3 n(Cluster_x[imeas],Cluster_y[imeas],0);
177  TVector3 v(0,0,1);
178  TVector3 u = v.Cross(n);
179 
180  PHGenFit::Measurement* meas = new PHGenFit::PlanarMeasurement(
181  pos,u,v,
182  Cluster_size_dphi[imeas],Cluster_size_dz[imeas]);
183  measurements.push_back(meas);
184  }
185 
187  track->addMeasurements(measurements);
188 
190  fitter->processTrack(track, false);
191 
193  genfit::MeasuredStateOnPlane* state_at_beam_line = track->extrapolateToLine(
194  TVector3(0, 0, 0), TVector3(0, 0, 1));
195  //state_at_beam_line->Print();
196 
197 // genfit::MeasuredStateOnPlane* state_at_layer_6 = track->extrapolateToCylinder(80.,
198 // TVector3(0, 0, 0), TVector3(0, 0, 1));
199 // //state_at_layer_6->Print();
200 // delete state_at_layer_6;
201 
202 
203  TVector3 GenFit_mom = state_at_beam_line->getMom();
204 
205  TVector3 AlanDion_mom(AlanDion_px,AlanDion_py,AlanDion_pz);
206 
207  hpT_residual_vs_pT->Fill(True_mom.Pt(),(GenFit_mom.Pt() - True_mom.Pt())/True_mom.Pt());
208 
209  hDCAr_vs_pT->Fill(True_mom.Mag(),state_at_beam_line->getState()[3]);
210 
211  delete state_at_beam_line;
212  delete track;
213  measurements.clear();
214  }
215 
216  gStyle->SetOptFit();
217  gStyle->SetOptStat(000000000);
218 
219  TF1 *tf_pT_resolution = new TF1("tf_pT_resolution","sqrt([0]*[0] + x*x*[1]*[1])", 0, 40);
220  tf_pT_resolution->SetParameters(0,0);
221  TCanvas *c3 = new TCanvas("c3","c3");
222  c3->Divide(2,1);
223  c3->cd(1);
224  hpT_residual_vs_pT->FitSlicesY();
225  TH1D *hpT_resolution_vs_pT = (TH1D*)gDirectory->Get("hpT_residual_vs_pT_2");
226  hpT_resolution_vs_pT->SetTitle("PHGenFit: #sigma_{p_{T}}/p_{T}; p_{T}[GeV/c]; #sigma_{p_{T}}/p_{T}");
227  hpT_resolution_vs_pT->SetMarkerStyle(20);
228  hpT_resolution_vs_pT->Draw("e");
229  hpT_resolution_vs_pT->Fit(tf_pT_resolution);
230  c3->cd(2);
231  hDCAr_vs_pT->FitSlicesY();
232  TH1D *hDCAr_resolution_vs_pT = (TH1D*) gDirectory->Get("hDCAr_vs_pT_2");
233  hDCAr_resolution_vs_pT->SetTitle(
234  "PHGenFit: #sigma_{DCAr} [cm]; p [GeV/c]; #sigma_{DCAr}");
235  hDCAr_resolution_vs_pT->SetMarkerStyle(20);
236  hDCAr_resolution_vs_pT->Draw("e");
237  //hDCAr_resolution_vs_pT->Fit(tf_pT_resolution);
238  c3->Print("pT_DCA_resolution.root");
239 
241  //fitter->displayEvent();
242 
243  fPHG4Hits->Close();
244 
245  delete fitter;
246 
247  //pause();
248 
249  std::cout<<"SUCCESS! \n";
250 
251  return 0;
252 }