hpstr
The Heavy Photon Search Toolkit for Reconstruction (hpstr) provides an interface to physics data from the HPS experiment saved in the LCIO format and converts it into an ROOT based format.
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
ThreeProngHistos.cxx
Go to the documentation of this file.
1#include "ThreeProngHistos.h"
2#include "TLorentzVector.h"
3#include "TVector3.h"
4#include <iostream>
5
6
7//fill plots for 3-prong tridents....
8//the "ele" always has higher energy cluster than "recoil"
10 Particle* pos,
11 Particle* recoil,
12 float weight) {
13 CalCluster eleClu=ele->getCluster();
14 CalCluster posClu=pos->getCluster();
15 CalCluster recClu=recoil->getCluster();
16 Track eleTrk=ele->getTrack();
17 Track posTrk=pos->getTrack();
18 Track recTrk=recoil->getTrack();
19 bool hasElectronTrack=false;
20 bool hasPositronTrack=false;
21 bool hasRecoilTrack=false;
22 if(eleTrk.getTrackerHitCount()>0)
23 hasElectronTrack=true;
24 if(posTrk.getTrackerHitCount()>0)
25 hasPositronTrack=true;
26 if(recTrk.getTrackerHitCount()>0)
27 hasRecoilTrack=true;
28 if(hasElectronTrack)
29 if(eleTrk.getCharge()==1){
30 std::cout<<"FillThreeProngPlots::Oops...positron pointing to electron-side cluster...lets skip it...return without filling"<<std::endl;
31 return;
32 }
33
34 if(hasPositronTrack)
35 if(posTrk.getCharge()==-1){
36 std::cout<<"FillThreeProngPlots::Oops...electron pointing to positron-side cluster...lets skip it...return without filling"<<std::endl;
37 return;
38 }
39
40 if(hasRecoilTrack)
41 if(recTrk.getCharge()==1){
42 std::cout<<"FillThreeProngPlots::Oops...positron pointing to recoil cluster...lets skip it for now...return without filling"<<std::endl;
43 return;
44 }
45
46
47 int evtCharge=0;
48 if(hasElectronTrack)
49 evtCharge+=eleTrk.getCharge();
50 if(hasPositronTrack)
51 evtCharge+=posTrk.getCharge();
52 if(hasRecoilTrack)
53 evtCharge+=recTrk.getCharge();
54
55
56 int eleTrkNHits=eleTrk.getTrackerHitCount();
57 int posTrkNHits=posTrk.getTrackerHitCount();
58 int recTrkNHits=recTrk.getTrackerHitCount();
59 double eleE=eleClu.getEnergy();
60 double eleX=eleClu.getPosition().at(0);
61 double eleY=eleClu.getPosition().at(1);
62 double eleZ=eleClu.getPosition().at(2);
63 double posE=posClu.getEnergy();
64 double posX=posClu.getPosition().at(0);
65 double posY=posClu.getPosition().at(1);
66 double posZ=posClu.getPosition().at(2);
67 double recE=recClu.getEnergy();
68 double recX=recClu.getPosition().at(0);
69 double recY=recClu.getPosition().at(1);
70 double recZ=recClu.getPosition().at(2);
71 double esum=eleE+posE+recE;
72 //track momentum
73 TVector3 eleP(eleTrk.getMomentum()[0],eleTrk.getMomentum()[1],eleTrk.getMomentum()[2]);
74 TVector3 posP(posTrk.getMomentum()[0],posTrk.getMomentum()[1],posTrk.getMomentum()[2]);
75 TVector3 recP(recTrk.getMomentum()[0],recTrk.getMomentum()[1],recTrk.getMomentum()[2]);
76
77 double eleTrkMom=eleP.Mag();
78 double posTrkMom=posP.Mag();
79 double recTrkMom=recP.Mag();
80 if(eleTrkMom>666)
81 eleP.SetXYZ(0,0,0);
82 if(posTrkMom>666)
83 posP.SetXYZ(0,0,0);
84 if(recTrkMom>666)
85 recP.SetXYZ(0,0,0);
86
87 TVector3 sumP=eleP+posP+recP;
88
89 double pSum=sumP.Mag();
90 double pYSum=sumP.Py();
91 double pXSum=sumP.Px();
92 //std::cout<<"p Sum = "<<pSum<<"; E beam = "<<eBeam<<std::endl;
93 //roughly the amount of "energy in E"...true for photons, approximation for charged tracks
94 double eleClYEne=sin(atan2(eleY,eleZ))*eleE;
95 double posClYEne=sin(atan2(posY,posZ))*posE;
96 double recClYEne=sin(atan2(recY,recZ))*recE;
97 double netClYEne=posClYEne+eleClYEne+recClYEne;
98 bool electronsSameHalf=false;
99 if(eleY*recY>0)
100 electronsSameHalf=true;
101
102
103 Fill1DHisto("clE_sum_allpos_allele_allrec_h",esum);
104 Fill1DHisto("clNet_EY_allpos_allele_allrec_h",netClYEne);
105 Fill1DHisto("clE_ele_allpos_allele_allrec_h",eleE);
106 Fill1DHisto("clX_ele_allpos_allele_allrec_h",eleX);
107 Fill1DHisto("clY_ele_allpos_allele_allrec_h",eleY);
108 Fill1DHisto("clE_pos_allpos_allele_allrec_h",posE);
109 Fill1DHisto("clX_pos_allpos_allele_allrec_h",posX);
110 Fill1DHisto("clY_pos_allpos_allele_allrec_h",posY);
111 Fill1DHisto("clE_rec_allpos_allele_allrec_h",recE);
112 Fill1DHisto("clX_rec_allpos_allele_allrec_h",recX);
113 Fill1DHisto("clY_rec_allpos_allele_allrec_h",recY);
114
115 if(hasRecoilTrack&&hasPositronTrack&&hasElectronTrack){
116 Fill1DHisto("clE_sum_foundpos_foundele_foundrec_h",esum);
117 Fill1DHisto("clNet_EY_foundpos_foundele_foundrec_h",netClYEne);
118 Fill1DHisto("clE_ele_foundpos_foundele_foundrec_h",eleE);
119 Fill1DHisto("clX_ele_foundpos_foundele_foundrec_h",eleX);
120 Fill1DHisto("clY_ele_foundpos_foundele_foundrec_h",eleY);
121 Fill1DHisto("clE_pos_foundpos_foundele_foundrec_h",posE);
122 Fill1DHisto("clX_pos_foundpos_foundele_foundrec_h",posX);
123 Fill1DHisto("clY_pos_foundpos_foundele_foundrec_h",posY);
124 Fill1DHisto("clE_rec_foundpos_foundele_foundrec_h",recE);
125 Fill1DHisto("clX_rec_foundpos_foundele_foundrec_h",recX);
126 Fill1DHisto("clY_rec_foundpos_foundele_foundrec_h",recY);
127
128 Fill2DHisto("clY_vs_clE_ele_foundpos_foundele_foundrec_hh",eleE,eleY);
129 Fill2DHisto("clY_vs_clX_ele_foundpos_foundele_foundrec_hh",eleX,eleY);
130 Fill2DHisto("clE_vs_clX_ele_foundpos_foundele_foundrec_hh",eleX,eleE);
131 Fill2DHisto("clY_vs_clE_pos_foundpos_foundele_foundrec_hh",posE,posY);
132 Fill2DHisto("clY_vs_clX_pos_foundpos_foundele_foundrec_hh",posX,posY);
133 Fill2DHisto("clE_vs_clX_pos_foundpos_foundele_foundrec_hh",posX,posE);
134 Fill2DHisto("clY_vs_clE_rec_foundpos_foundele_foundrec_hh",recE,recY);
135 Fill2DHisto("clY_vs_clX_rec_foundpos_foundele_foundrec_hh",recX,recY);
136 Fill2DHisto("clE_vs_clX_rec_foundpos_foundele_foundrec_hh",recX,recE);
137
138 Fill1DHisto("nHits_ele_foundpos_foundele_foundrec_h",eleTrkNHits,weight);
139 Fill1DHisto("trkMom_ele_foundpos_foundele_foundrec_h",eleTrkMom,weight);
140 Fill1DHisto("trkEoverP_ele_foundpos_foundele_foundrec_h",eleE/eleTrkMom,weight);
141 Fill1DHisto("nHits_pos_foundpos_foundele_foundrec_h",posTrkNHits,weight);
142 Fill1DHisto("trkMom_pos_foundpos_foundele_foundrec_h",posTrkMom,weight);
143 Fill1DHisto("trkEoverP_pos_foundpos_foundele_foundrec_h",posE/posTrkMom,weight);
144 Fill1DHisto("nHits_rec_foundpos_foundele_foundrec_h",recTrkNHits,weight);
145 Fill1DHisto("trkMom_rec_foundpos_foundele_foundrec_h",recTrkMom,weight);
146 Fill1DHisto("trkEoverP_rec_foundpos_foundele_foundrec_h",recE/recTrkMom,weight);
147
148 Fill1DHisto("trkP_sum_foundpos_foundele_foundrec_h",pSum);
149 Fill1DHisto("trkPX_sum_foundpos_foundele_foundrec_h",pXSum);
150 Fill1DHisto("trkPY_sum_foundpos_foundele_foundrec_h",pYSum);
151
152 Fill2DHisto("clE_sum_vs_trkP_sum_foundpos_foundele_foundrec_hh",pSum,esum);
153
154 }
155 if(!hasRecoilTrack&&hasPositronTrack&&hasElectronTrack){
156 Fill1DHisto("clE_sum_foundpos_foundele_missrec_h",esum);
157 Fill1DHisto("clNet_EY_foundpos_foundele_missrec_h",netClYEne);
158 Fill1DHisto("clE_ele_foundpos_foundele_missrec_h",eleE);
159 Fill1DHisto("clX_ele_foundpos_foundele_missrec_h",eleX);
160 Fill1DHisto("clY_ele_foundpos_foundele_missrec_h",eleY);
161 Fill1DHisto("clE_pos_foundpos_foundele_missrec_h",posE);
162 Fill1DHisto("clX_pos_foundpos_foundele_missrec_h",posX);
163 Fill1DHisto("clY_pos_foundpos_foundele_missrec_h",posY);
164 Fill1DHisto("clE_rec_foundpos_foundele_missrec_h",recE);
165 Fill1DHisto("clX_rec_foundpos_foundele_missrec_h",recX);
166 Fill1DHisto("clY_rec_foundpos_foundele_missrec_h",recY);
167
168
169 Fill2DHisto("clY_vs_clE_ele_foundpos_foundele_missrec_hh",eleE,eleY);
170 Fill2DHisto("clY_vs_clX_ele_foundpos_foundele_missrec_hh",eleX,eleY);
171 Fill2DHisto("clE_vs_clX_ele_foundpos_foundele_missrec_hh",eleX,eleE);
172 Fill2DHisto("clY_vs_clE_pos_foundpos_foundele_missrec_hh",posE,posY);
173 Fill2DHisto("clY_vs_clX_pos_foundpos_foundele_missrec_hh",posX,posY);
174 Fill2DHisto("clE_vs_clX_pos_foundpos_foundele_missrec_hh",posX,posE);
175 Fill2DHisto("clY_vs_clE_rec_foundpos_foundele_missrec_hh",recE,recY);
176 Fill2DHisto("clY_vs_clX_rec_foundpos_foundele_missrec_hh",recX,recY);
177 Fill2DHisto("clE_vs_clX_rec_foundpos_foundele_missrec_hh",recX,recE);
178
179 Fill1DHisto("nHits_ele_foundpos_foundele_missrec_h",eleTrkNHits,weight);
180 Fill1DHisto("trkMom_ele_foundpos_foundele_missrec_h",eleTrkMom,weight);
181 Fill1DHisto("trkEoverP_ele_foundpos_foundele_missrec_h",eleE/eleTrkMom,weight);
182 Fill1DHisto("nHits_pos_foundpos_foundele_missrec_h",posTrkNHits,weight);
183 Fill1DHisto("trkMom_pos_foundpos_foundele_missrec_h",posTrkMom,weight);
184 Fill1DHisto("trkEoverP_pos_foundpos_foundele_missrec_h",posE/posTrkMom,weight);
185
186 Fill1DHisto("trkP_miss_foundpos_foundele_missrec_h",eBeam-pSum,weight); //this should be recoil p
187 Fill1DHisto("trkPX_miss_foundpos_foundele_missrec_h",-pXSum,weight); //I think I need to rotate
188 Fill1DHisto("trkPY_miss_foundpos_foundele_missrec_h",-pYSum,weight); //this should be -ive recoil pY
189 }
190 if(hasRecoilTrack&&!hasPositronTrack&&hasElectronTrack){
191 Fill1DHisto("clE_sum_misspos_foundele_foundrec_h",esum);
192 Fill1DHisto("clNet_EY_misspos_foundele_foundrec_h",netClYEne);
193 Fill1DHisto("clE_ele_misspos_foundele_foundrec_h",eleE);
194 Fill1DHisto("clX_ele_misspos_foundele_foundrec_h",eleX);
195 Fill1DHisto("clY_ele_misspos_foundele_foundrec_h",eleY);
196 Fill1DHisto("clE_pos_misspos_foundele_foundrec_h",posE);
197 Fill1DHisto("clX_pos_misspos_foundele_foundrec_h",posX);
198 Fill1DHisto("clY_pos_misspos_foundele_foundrec_h",posY);
199 Fill1DHisto("clE_rec_misspos_foundele_foundrec_h",recE);
200 Fill1DHisto("clX_rec_misspos_foundele_foundrec_h",recX);
201 Fill1DHisto("clY_rec_misspos_foundele_foundrec_h",recY);
202
203
204 Fill2DHisto("clY_vs_clE_ele_misspos_foundele_foundrec_hh",eleE,eleY);
205 Fill2DHisto("clY_vs_clX_ele_misspos_foundele_foundrec_hh",eleX,eleY);
206 Fill2DHisto("clE_vs_clX_ele_misspos_foundele_foundrec_hh",eleX,eleE);
207 Fill2DHisto("clY_vs_clE_pos_misspos_foundele_foundrec_hh",posE,posY);
208 Fill2DHisto("clY_vs_clX_pos_misspos_foundele_foundrec_hh",posX,posY);
209 Fill2DHisto("clE_vs_clX_pos_misspos_foundele_foundrec_hh",posX,posE);
210 Fill2DHisto("clY_vs_clE_rec_misspos_foundele_foundrec_hh",recE,recY);
211 Fill2DHisto("clY_vs_clX_rec_misspos_foundele_foundrec_hh",recX,recY);
212 Fill2DHisto("clE_vs_clX_rec_misspos_foundele_foundrec_hh",recX,recE);
213
214 Fill1DHisto("nHits_ele_misspos_foundele_foundrec_h",eleTrkNHits,weight);
215 Fill1DHisto("trkMom_ele_misspos_foundele_foundrec_h",eleTrkMom,weight);
216 Fill1DHisto("trkEoverP_ele_misspos_foundele_foundrec_h",eleE/eleTrkMom,weight);
217 Fill1DHisto("nHits_rec_misspos_foundele_foundrec_h",recTrkNHits,weight);
218 Fill1DHisto("trkMom_rec_misspos_foundele_foundrec_h",recTrkMom,weight);
219 Fill1DHisto("trkEoverP_rec_misspos_foundele_foundrec_h",recE/recTrkMom,weight);
220
221 Fill1DHisto("trkP_miss_misspos_foundele_foundrec_h",eBeam-pSum,weight); //this should be recoil p
222 Fill1DHisto("trkPX_miss_misspos_foundele_foundrec_h",-pXSum,weight); //I think I need to rotate
223 Fill1DHisto("trkPY_miss_misspos_foundele_foundrec_h",-pYSum,weight); //this should be -ive recoil pY
224 }
225 if(hasRecoilTrack&&hasPositronTrack&&!hasElectronTrack){
226 Fill1DHisto("clE_sum_foundpos_missele_foundrec_h",esum);
227 Fill1DHisto("clNet_EY_foundpos_missele_foundrec_h",netClYEne);
228 Fill1DHisto("clE_ele_foundpos_missele_foundrec_h",eleE);
229 Fill1DHisto("clX_ele_foundpos_missele_foundrec_h",eleX);
230 Fill1DHisto("clY_ele_foundpos_missele_foundrec_h",eleY);
231 Fill1DHisto("clE_pos_foundpos_missele_foundrec_h",posE);
232 Fill1DHisto("clX_pos_foundpos_missele_foundrec_h",posX);
233 Fill1DHisto("clY_pos_foundpos_missele_foundrec_h",posY);
234 Fill1DHisto("clE_rec_foundpos_missele_foundrec_h",recE);
235 Fill1DHisto("clX_rec_foundpos_missele_foundrec_h",recX);
236 Fill1DHisto("clY_rec_foundpos_missele_foundrec_h",recY);
237
238
239 Fill2DHisto("clY_vs_clE_ele_foundpos_missele_foundrec_hh",eleE,eleY);
240 Fill2DHisto("clY_vs_clX_ele_foundpos_missele_foundrec_hh",eleX,eleY);
241 Fill2DHisto("clE_vs_clX_ele_foundpos_missele_foundrec_hh",eleX,eleE);
242 Fill2DHisto("clY_vs_clE_pos_foundpos_missele_foundrec_hh",posE,posY);
243 Fill2DHisto("clY_vs_clX_pos_foundpos_missele_foundrec_hh",posX,posY);
244 Fill2DHisto("clE_vs_clX_pos_foundpos_missele_foundrec_hh",posX,posE);
245 Fill2DHisto("clY_vs_clE_rec_foundpos_missele_foundrec_hh",recE,recY);
246 Fill2DHisto("clY_vs_clX_rec_foundpos_missele_foundrec_hh",recX,recY);
247 Fill2DHisto("clE_vs_clX_rec_foundpos_missele_foundrec_hh",recX,recE);
248
249 Fill1DHisto("nHits_pos_foundpos_missele_foundrec_h",posTrkNHits,weight);
250 Fill1DHisto("trkMom_pos_foundpos_missele_foundrec_h",posTrkMom,weight);
251 Fill1DHisto("trkEoverP_pos_foundpos_missele_foundrec_h",posE/posTrkMom,weight);
252 Fill1DHisto("nHits_rec_foundpos_missele_foundrec_h",recTrkNHits,weight);
253 Fill1DHisto("trkMom_rec_foundpos_missele_foundrec_h",recTrkMom,weight);
254 Fill1DHisto("trkEoverP_rec_foundpos_missele_foundrec_h",recE/recTrkMom,weight);
255
256 Fill1DHisto("trkP_miss_foundpos_missele_foundrec_h",eBeam-pSum,weight); //this should be recoil p
257 Fill1DHisto("trkPX_miss_foundpos_missele_foundrec_h",-pXSum,weight); //I think I need to rotate
258 Fill1DHisto("trkPY_miss_foundpos_missele_foundrec_h",-pYSum,weight); //this should be -ive recoil pY
259
260
261 }
262 if(!hasRecoilTrack&&hasPositronTrack&&!hasElectronTrack){
263 Fill1DHisto("clE_sum_foundpos_missele_missrec_h",esum);
264 Fill1DHisto("clNet_EY_foundpos_missele_missrec_h",netClYEne);
265 Fill1DHisto("clE_ele_foundpos_missele_missrec_h",eleE);
266 Fill1DHisto("clX_ele_foundpos_missele_missrec_h",eleX);
267 Fill1DHisto("clY_ele_foundpos_missele_missrec_h",eleY);
268 Fill1DHisto("clE_pos_foundpos_missele_missrec_h",posE);
269 Fill1DHisto("clX_pos_foundpos_missele_missrec_h",posX);
270 Fill1DHisto("clY_pos_foundpos_missele_missrec_h",posY);
271 Fill1DHisto("clE_rec_foundpos_missele_missrec_h",recE);
272 Fill1DHisto("clX_rec_foundpos_missele_missrec_h",recX);
273 Fill1DHisto("clY_rec_foundpos_missele_missrec_h",recY);
274
275 Fill1DHisto("nHits_pos_foundpos_missele_missrec_h",posTrkNHits,weight);
276 Fill1DHisto("trkMom_pos_foundpos_missele_missrec_h",posTrkMom,weight);
277 Fill1DHisto("trkEoverP_pos_foundpos_missele_missrec_h",posE/posTrkMom,weight);
278 }
279 if(!hasRecoilTrack&&!hasPositronTrack&&!hasElectronTrack){
280 Fill1DHisto("clE_sum_misspos_missele_missrec_h",esum);
281 Fill1DHisto("clNet_EY_misspos_missele_missrec_h",netClYEne);
282 Fill1DHisto("clE_ele_misspos_missele_missrec_h",eleE);
283 Fill1DHisto("clX_ele_misspos_missele_missrec_h",eleX);
284 Fill1DHisto("clY_ele_misspos_missele_missrec_h",eleY);
285 Fill1DHisto("clE_pos_misspos_missele_missrec_h",posE);
286 Fill1DHisto("clX_pos_misspos_missele_missrec_h",posX);
287 Fill1DHisto("clY_pos_misspos_missele_missrec_h",posY);
288 Fill1DHisto("clE_rec_misspos_missele_missrec_h",recE);
289 Fill1DHisto("clX_rec_misspos_missele_missrec_h",recX);
290 Fill1DHisto("clY_rec_misspos_missele_missrec_h",recY);
291
292
293 Fill2DHisto("clY_vs_clE_ele_misspos_missele_missrec_hh",eleE,eleY);
294 Fill2DHisto("clY_vs_clX_ele_misspos_missele_missrec_hh",eleX,eleY);
295 Fill2DHisto("clE_vs_clX_ele_misspos_missele_missrec_hh",eleX,eleE);
296 Fill2DHisto("clY_vs_clE_pos_misspos_missele_missrec_hh",posE,posY);
297 Fill2DHisto("clY_vs_clX_pos_misspos_missele_missrec_hh",posX,posY);
298 Fill2DHisto("clE_vs_clX_pos_misspos_missele_missrec_hh",posX,posE);
299 Fill2DHisto("clY_vs_clE_rec_misspos_missele_missrec_hh",recE,recY);
300 Fill2DHisto("clY_vs_clX_rec_misspos_missele_missrec_hh",recX,recY);
301 Fill2DHisto("clE_vs_clX_rec_misspos_missele_missrec_hh",recX,recE);
302 }
303}
double getEnergy() const
Definition CalCluster.h:73
std::vector< double > getPosition() const
Definition CalCluster.h:63
void Fill2DHisto(const std::string &histoName, float valuex, float valuey, float weight=1.)
description
void Fill1DHisto(const std::string &histoName, float value, float weight=1.)
description
Track getTrack() const
Definition Particle.h:48
CalCluster getCluster() const
Definition Particle.h:62
void FillThreeProngPlots(Particle *ele, Particle *pos, Particle *rec, float weight=1.0)
Definition Track.h:32
int getTrackerHitCount() const
Definition Track.h:345
int getCharge() const
Definition Track.h:247
std::vector< double > getMomentum()
Definition Track.h:285