10 #include <g4hough/SvtxVertexMap.h>
11 #include <g4hough/SvtxVertex.h>
12 #include <g4hough/SvtxTrackMap.h>
13 #include <g4hough/SvtxTrack.h>
14 #include <g4hough/SvtxClusterMap.h>
15 #include <g4hough/SvtxCluster.h>
36 fReconstructable_TPCHits(30),
37 fFair_ClustersContribution(20),
39 fGTrack_TPCClusters(20)
47 const int nptbinsPONE = 52;
48 Int_t nptbins = nptbinsPONE - 1;
49 Double_t ptbins[nptbinsPONE];
51 for(
int i=1;
i!=9; ++
i)
52 ptbins[
i] = 0.2 +
i*0.2;
53 for(
int i=0;
i!=11; ++
i)
54 ptbins[9+
i] = 2.0 +
i*0.5;
55 for(
int i=0;
i!=10; ++
i)
56 ptbins[20+
i] = 8.0 +
i*1.0;
57 for(
int i=0;
i!=12; ++
i)
58 ptbins[30+
i] = 18.0 +
i*2.0;
59 for(
int i=0;
i!=10; ++
i)
60 ptbins[42+
i] = 45.0 +
i*5.0;
62 for(
int i=0;
i!=nptbinsPONE; ++
i) std::cout <<
i <<
" " << ptbins[
i] << std::endl;
64 fHNEvents =
new TH1F(
"Events",
"Events",1,-0.5,0.5);
67 TString sTname[5] = {
"AllG4Particles",
"Primaries",
"Embedded",
"Reconstructables",
"RecoMatched"};
68 for(
int i=0;
i!=5; ++
i) {
69 fHTN[
i] =
new TH1F(Form(
"Truths%d_N",
i),Form(
"Number of %s",sTname[
i].
Data()),50,-0.5,49.5);
70 fHTPt[
i] =
new TH1F(Form(
"Truths%d_Pt",
i),Form(
"%s Pt",sTname[
i].
Data()),nptbins,ptbins);
72 fHTEta[
i] =
new TH1F(Form(
"Truths%d_Eta",
i),Form(
"%s Eta",sTname[
i].
Data()),100,-10,+10);
78 fHTNHits[
i] =
new TH1F(Form(
"Truths%d_Hits",i),Form(
"Number of TPCHits in %s",sTname[i].
Data()),61,-0.5,60.5);
81 fHTChi2[
i] =
new TH2F(Form(
"Truths%d_Chi2NDF",i),Form(
"%s Chi2NDF;PT;CHI2NDF",sTname[i].
Data()),nptbins,ptbins,100,0,5);
82 fHTDca2D[
i] =
new TH2F( Form(
"Truths%d_Dca2D",i),Form(
"%s Dca2D;PT;DCA2D",sTname[i].
Data()),nptbins,ptbins,200,-0.003,+0.003);
83 fHTNClustersContribution[
i] =
new TH2F(Form(
"Truths%d_NClusters",i), Form(
"%s Number of Clusters Matched",sTname[i].
Data()), nptbins,ptbins,70,-0.5,69.5);
84 fHTPtResolution[
i] =
new TH2F(Form(
"Truths%d_ResPt",i), Form(
"%s PtResolution;PT;REL DIFF",sTname[i].
Data()), nptbins,ptbins,400,-1.1,+1.1);
85 fHTPtResolution2[
i] =
new TH2F(Form(
"Truths%d_Res2Pt",i), Form(
"%s PtResolution2;PT;DIFF",sTname[i].
Data()), nptbins,ptbins,400,-0.2,+0.2);
86 fHTPhiResolution[
i] =
new TH2F(Form(
"Truths%d_ResPhi",i), Form(
"%s PhiResolution;PT;REL DIFF",sTname[i].
Data()), nptbins,ptbins,200,-0.1,+0.1);
87 fHTEtaResolution[
i] =
new TH2F(Form(
"Truths%d_ResEta",i), Form(
"%s EtaResolution;PT;REL DIFF",sTname[i].
Data()), nptbins,ptbins,200,-0.1,+0.1);
97 TString sRname[4] = {
"AllTracks",
"GoodTracks",
"GTracksM2R",
"GTracksExcess"};
98 for(
int i=0;
i!=4; ++
i) {
99 fHRN[
i] =
new TH1F(Form(
"Tracks%d_N",
i),Form(
"Number of %s",sRname[
i].
Data()),100,-0.5,1999.5);
100 fHRPt[
i] =
new TH1F(Form(
"Tracks%d_Pt",
i),Form(
"%s Pt",sRname[
i].
Data()),nptbins,ptbins);
102 fHREta[
i] =
new TH1F(Form(
"Tracks%d_Eta",
i),Form(
"%s Eta",sRname[
i].
Data()),100,-1.5,+1.5);
103 fHRTPCClus[
i] =
new TH2F(Form(
"Tracks%d_TPCClus",
i),Form(
"%s TPCCLus;PT;TPCClusters",sRname[
i].
Data()),nptbins,ptbins,70,-0.5,69.5);
104 fHRChi2[
i] =
new TH2F(Form(
"Tracks%d_Chi2NDF",
i),Form(
"%s Chi2NDF;PT;CHI2NDF",sRname[
i].
Data()),nptbins,ptbins,100,0,5);
105 fHRDca2D[
i] =
new TH2F( Form(
"Tracks%d_Dca2D",
i),Form(
"%s Dca2D;PT;DCA2D",sRname[
i].
Data()),nptbins,ptbins,200,-0.003,+0.003);
114 fHRPtResolution[
i] =
new TH2F(Form(
"Tracks%d_ResPt",i), Form(
"%s PtResolution;PT;REL DIFF",sRname[i].
Data()), nptbins,ptbins,400,-1.5,+1.5);
115 fHRPtResolution2[
i] =
new TH2F(Form(
"Tracks%d_Res2Pt",i), Form(
"%s PtResolution2;PT;DIFF",sRname[i].
Data()), nptbins,ptbins,400,-0.2,+0.2);
116 fHRPhiResolution[
i] =
new TH2F(Form(
"Tracks%d_ResPhi",i), Form(
"%s PhiResolution;PT;REL DIFF",sRname[i].
Data()), nptbins,ptbins,200,-0.1,+0.1);
117 fHREtaResolution[
i] =
new TH2F(Form(
"Tracks%d_ResEta",i), Form(
"%s EtaResolution;PT;REL DIFF",sRname[i].
Data()), nptbins,ptbins,200,-0.1,+0.1);
118 fHRNClustersContribution[
i] =
new TH2F(Form(
"Tracks%d_NClustersContri",i), Form(
"%s Number of Clusters Contribution",sRname[i].
Data()), nptbins,ptbins,70,-0.5,69.5);
126 fHNVertexes =
new TH1F(
"Vertexes_N",
"Number of Vertexes",10,-0.5,9.5);
136 SvtxTrackMap *trackmap = findNode::getClass<SvtxTrackMap>(topNode,
"SvtxTrackMap");
137 SvtxClusterMap *clustermap = findNode::getClass<SvtxClusterMap>(topNode,
"SvtxClusterMap");
138 SvtxVertexMap* vertexmap = findNode::getClass<SvtxVertexMap>(topNode,
"SvtxVertexMap");
140 cerr <<
PHWHERE <<
" ERROR: Can't find G4TruthInfo" << endl;
144 cerr <<
PHWHERE <<
" ERROR: Can't find SvtxTrackMap" << endl;
148 cerr <<
PHWHERE <<
" ERROR: Can't find SvtxVertexMap" << endl;
153 std::pair< std::map<int,int>::const_iterator, std::map<int,int>::const_iterator > rangeE = truthinfo->
GetEmbeddedTrkIds();
154 for(std::map<int,int>::const_iterator iterE=rangeE.first; iterE!=rangeE.second; ++iterE) {
155 int tid = (*iterE).first;
156 int tem = (*iterE).second;
158 fEmbedded.insert( std::make_pair(tid,tem) );
172 iter != range.second;
175 if(!g4particle)
continue;
177 float px = g4particle->
get_px();
178 float py = g4particle->
get_py();
179 float pz = g4particle->
get_pz();
180 float pt = TMath::Sqrt( px*px + py*py );
181 float p = TMath::Sqrt( px*px + py*py + pz*pz );
182 float phi = TMath::Pi()+TMath::ATan2(-py,-px);
184 if(p != TMath::Abs(pz)) eta = 0.5*
TMath::Log((p+pz)/(p-pz));
188 if( iter->first < 0 )
continue;
194 std::map<int,int>::iterator embeddedflag =
fEmbedded.find(
id);
195 if(embeddedflag==
fEmbedded.end())
continue;
196 if((*embeddedflag).second==0)
continue;
200 std::set<PHG4Hit*> g4hits = trutheval->
all_truth_hits(g4particle);
202 for(
int i=0;
i!=70; ++
i) hit[
i] =
false;
203 for(std::set<PHG4Hit*>::iterator iter = g4hits.begin();
204 iter != g4hits.end();
207 unsigned int lyr = candidate->
get_layer();
218 (*embeddedflag).second = 2;
230 float rpt = TMath::Sqrt( rpx*rpx + rpy*rpy );
231 float rp = TMath::Sqrt( rpx*rpx + rpy*rpy + rpz*rpz );
232 float rphi = TMath::Pi()+TMath::ATan2(-rpy,-rpx);
234 if(rp != TMath::Abs(rpz)) reta = 0.5*
TMath::Log((rp+rpz)/(rp-rpz));
238 int ndf = int(best->
get_ndf());
240 if(ndf!=0) chi2ndf = best->
get_chisq() / ndf;
241 float ptreldiff = (rpt-
pt)/pt;
242 float ptdiff = rpt-
pt;
243 float phireldiff = (rphi-
phi)/phi;
244 float etareldiff = (reta-
eta)/eta;
249 fHTChi2[4]->Fill(chi2ndf,ptreldiff);
269 iter != trackmap->end(); ++iter) {
273 if (!track)
continue;
278 unsigned int cluster_id = *iter;
279 SvtxCluster* cluster = clustermap->get(cluster_id);
280 if(!cluster)
continue;
281 int lyr = cluster->get_layer();
287 float rpt = TMath::Sqrt( rpx*rpx + rpy*rpy );
288 float rp = TMath::Sqrt( rpx*rpx + rpy*rpy + rpz*rpz );
289 float rphi = TMath::Pi()+TMath::ATan2(-rpy,-rpx);
291 if(rp != TMath::Abs(rpz)) reta = 0.5*
TMath::Log((rp+rpz)/(rp-rpz));
293 int ndf = int(track->
get_ndf());
295 if(ndf!=0) chi2ndf = track->
get_chisq() / ndf;
314 if(!g4particle)
continue;
321 std::map<int,int>::iterator embeddedflag =
fEmbedded.find(
id);
322 if(embeddedflag==
fEmbedded.end())
continue;
324 if((*embeddedflag).second==0)
continue;
326 if((*embeddedflag).second==1)
continue;
327 if((*embeddedflag).second>1) (*embeddedflag).second++;
331 float px = g4particle->
get_px();
332 float py = g4particle->
get_py();
333 float pz = g4particle->
get_pz();
334 float pt = TMath::Sqrt( px*px + py*py );
335 float p = TMath::Sqrt( px*px + py*py + pz*pz );
336 float phi = TMath::Pi()+TMath::ATan2(-py,-px);
338 if(p != TMath::Abs(pz)) eta = 0.5*
TMath::Log((p+pz)/(p-pz));
340 float ptreldiff = (rpt-
pt)/pt;
341 float ptdiff = rpt-
pt;
342 float phireldiff = (rphi-
phi)/phi;
343 float etareldiff = (reta-
eta)/eta;
355 if((*embeddedflag).second<4)
continue;
386 iter != vertexmap->end();