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
TridentHistos.cxx
Go to the documentation of this file.
1#include "TridentHistos.h"
2#include "TLorentzVector.h"
3#include "TVector3.h"
4#include <iostream>
5#include <regex>
6#include <cstddef>
7void TridentHistos::BuildAxes(){}
8
9void TridentHistos::DefineHistos(){
10 if(debug_)
11 std::cout<<" Making layer-hit histo"<<std::endl;
12 std::vector<std::string> layerTags;
13 unsigned char lbits=15;
14 int counter=15;
15 while(counter>=0){
16 std::string lstring="";
17 for (int i = 4; i < 8; i++) { //just doing sensor layers 1-4
18 // char* tmpSt;
19 // printf("%d", !!((lbits << i) & 0x80));
20 int tmpCh=!!((lbits << i) & 0x80);
21 std::string tmpSt="0";
22 if(tmpCh==1)
23 tmpSt="1";
24 lstring+=tmpSt;
25 }
26 // printf("\n");
27 layerTags.push_back(lstring);
28 counter--;
29 lbits--;
30 }
31
32 std::vector<std::string> halfTags{"","top","bot"};
33 std::vector<std::string> trkTags{"pos","ele"};
34 std::vector<std::string> trkCombos;
35 std::vector<std::string>::iterator halfIter = halfTags.begin();
36 for(halfIter; halfIter < halfTags.end(); halfIter++){
37 std::string half=*halfIter;
38 std::vector<std::string>::iterator trkIter = trkTags.begin();
39 for(trkIter; trkIter < trkTags.end(); trkIter++){
40 std::string trk=*trkIter;
41 std::string comboStr=half+"_"+trk;
42 if(half=="")
43 comboStr=trk;
44 if(debug_)
45 std::cout<<" splitting by "<<comboStr<<std::endl;
46 trkCombos.push_back(comboStr);
47 }
48 }
49
50 /*
51 std::vector<std::string> trkLayerCombos;
52 std::vector<std::string>::iterator trkIter = trkCombos.begin();
53 for(trkIter;trkIter<trkCombos.end();trkIter++){
54 std::string trkSt=*trkIter;
55 std::vector<std::string>::iterator lIter = layerTags.begin();
56 for(lIter;lIter<layerTags.end();lIter++){
57 std::string lSt=*lIter;
58 std::string comboStr="L"+lSt+"_"+trkSt;
59 if(trkSt=="")
60 comboStr="L"+lSt;
61 trkLayerCombos.push_back(comboStr);
62 // std::cout<<comboStr<<std::endl;
63
64 }
65 }
66 */
67 // std::vector<std::string> layerTagsForV0{"1111","0011","0111"};
68 std::vector<std::string> layerTagsForV0{"1111","0011","0111","1011"};
69 std::vector<std::string> trkLayerCombos;
70 std::vector<std::string>::iterator trkIter = trkCombos.begin();
71 for(trkIter;trkIter<trkCombos.end();trkIter++){
72 std::string trkSt=*trkIter;
73 std::vector<std::string>::iterator lIter = layerTagsForV0.begin();
74 for(lIter;lIter<layerTagsForV0.end();lIter++){
75 std::string l1St=*lIter;
76 std::vector<std::string>::iterator l2Iter = layerTagsForV0.begin();
77 for(l2Iter;l2Iter<layerTagsForV0.end();l2Iter++){
78 std::string l2St=*l2Iter;
79 std::string comboStr="pos"+l1St+"_ele"+l2St+"_"+trkSt;
80 trkLayerCombos.push_back(comboStr);
81 if(debug_)
82 std::cout<<" splitting by "<<comboStr<<std::endl;
83 }
84
85 }
86 }
87
88
89 std::vector<std::string> V0Combos{"","top_pos","bot_pos"};
90 std::vector<std::string> v0LayerCombos;
91 std::vector<std::string>::iterator v0Iter = V0Combos.begin();
92 for(v0Iter;v0Iter<V0Combos.end();v0Iter++){
93 std::string v0St=*v0Iter;
94 std::vector<std::string>::iterator lIter = layerTagsForV0.begin();
95 for(lIter;lIter<layerTagsForV0.end();lIter++){
96 std::string l1St=*lIter;
97 std::vector<std::string>::iterator l2Iter = layerTagsForV0.begin();
98 for(l2Iter;l2Iter<layerTagsForV0.end();l2Iter++){
99 std::string l2St=*l2Iter;
100 std::string comboStr="pos"+l1St+"_ele"+l2St+"_"+v0St;
101 if(v0St=="")
102 comboStr="pos"+l1St+"_ele"+l2St;
103 if(debug_)
104 std::cout<<" splitting by "<<comboStr<<std::endl;
105 v0LayerCombos.push_back(comboStr);
106 }
107 }
108 }
109
110
111
112 DefineHistosFromTemplateOnly(trkCombos,"_tc_h"); //should match both 1d and 2d histos
113 DefineHistosFromTemplateOnly(halfTags, "_hc_h"); //should match both 1d and 2d histos (though, I think there are o
114 // DefineHistosFromTemplateOnly(trkLayerCombos,"_tc_h"); //should match both 1d and 2d histos
115 DefineHistosFromTemplateOnly(V0Combos,"_vc_h"); //should match both 1d and 2d histos
116 // DefineHistosFromTemplateOnly(V0Combos,"_tc_h"); //should match both 1d and 2d histos
117 DefineHistosFromTemplateOnly(v0LayerCombos,"_vc_h"); //should match both 1d and 2d histos
118 DefineHistosFromTemplateOnly(trkLayerCombos,"_tc_h"); //should match both 1d and 2d histos
119
120 DefineOneTimeHistos();
121}
122
123
124void TridentHistos::Define2DHistos() {
125
126}//define 2dhistos
127
128void TridentHistos::AssignLayerCode(Track* ele_trk,
129 Track* pos_trk){
130 std::string eleLayerCode=getLayerCodeFromTrack(ele_trk);
131 std::string posLayerCode=getLayerCodeFromTrack(pos_trk);
132 this->layerCode="pos"+posLayerCode+"_ele"+eleLayerCode+"_";
133}
134
135void TridentHistos::Fill1DVertex(Vertex* vtx,
136 Particle* ele,
137 Particle* pos,
138 Track* ele_trk,
139 Track* pos_trk,
140 double trkTimeOffset,
141 float weight) {
142
143
144 if (ele_trk)
145 Fill1DTrack(ele_trk,trkTimeOffset,weight,"ele_");
146 if (pos_trk)
147 Fill1DTrack(pos_trk,trkTimeOffset,weight,"pos_");
148 // std::cout<<"TridentHistos::Fill1DVertex 1d track histos filled"<<std::endl;
149 std::string half="top";
150 if(pos_trk->getTanLambda()<0)
151 half="bot";
152 std::string posTag=half+"_pos_";
153 half="top";
154 if(ele_trk->getTanLambda()<0)
155 half="bot";
156 std::string eleTag=half+"_ele_";
157
158 // std::string eleLayerCode=getLayerCodeFromTrack(ele_trk);
159 // std::string posLayerCode=getLayerCodeFromTrack(pos_trk);
160 // std::string layerCode="pos"+eleLayerCode+"_ele"+posLayerCode+"_";
161
162 Fill1DHisto("vtx_chi2_vc_h" ,vtx->getChi2(),weight);
163 Fill1DHisto("vtx_X_vc_h" ,vtx->getX(),weight);
164 Fill1DHisto("vtx_Y_vc_h" ,vtx->getY(),weight);
165 Fill1DHisto("vtx_Z_vc_h" ,vtx->getZ(),weight);
166
167 Fill1DHisto(posTag+"vtx_chi2_vc_h" ,vtx->getChi2(),weight);
168 Fill1DHisto(posTag+"vtx_X_vc_h" ,vtx->getX(),weight);
169 Fill1DHisto(posTag+"vtx_Y_vc_h" ,vtx->getY(),weight);
170 Fill1DHisto(posTag+"vtx_Z_vc_h" ,vtx->getZ(),weight);
171
172 Fill1DHisto(layerCode+"vtx_chi2_vc_h" ,vtx->getChi2(),weight);
173 Fill1DHisto(layerCode+"vtx_X_vc_h" ,vtx->getX(),weight);
174 Fill1DHisto(layerCode+"vtx_Y_vc_h" ,vtx->getY(),weight);
175 Fill1DHisto(layerCode+"vtx_Z_vc_h" ,vtx->getZ(),weight);
176
177 Fill1DHisto(layerCode+posTag+"vtx_chi2_vc_h" ,vtx->getChi2(),weight);
178 Fill1DHisto(layerCode+posTag+"vtx_X_vc_h" ,vtx->getX(),weight);
179 Fill1DHisto(layerCode+posTag+"vtx_Y_vc_h" ,vtx->getY(),weight);
180 Fill1DHisto(layerCode+posTag+"vtx_Z_vc_h" ,vtx->getZ(),weight);
181
182 TVector3 vtxPosSvt;
183 vtxPosSvt.SetX(vtx->getX());
184 vtxPosSvt.SetY(vtx->getY());
185 vtxPosSvt.SetZ(vtx->getZ());
186
187 vtxPosSvt.RotateY(-0.0305);
188
189 //Fill1DHisto("vtx_X_svt_vc_h",vtxPosSvt.X(),weight);
190 //Fill1DHisto("vtx_Y_svt_vc_h",vtxPosSvt.Y(),weight);
191 //Fill1DHisto("vtx_Z_svt_vc_h",vtxPosSvt.Z(),weight);
192
193 // 0 xx 1 xy 2 xz 3 yy 4 yz 5 zz
194 Fill1DHisto("vtx_sigma_X_vc_h",sqrt(vtx->getCovariance()[0]),weight);
195 Fill1DHisto("vtx_sigma_Y_vc_h",sqrt(vtx->getCovariance()[3]),weight);
196 Fill1DHisto("vtx_sigma_Z_vc_h",sqrt(vtx->getCovariance()[5]),weight);
197 Fill1DHisto("vtx_InvM_vc_h" ,vtx->getInvMass(),weight);
198 Fill1DHisto("vtx_InvMErr_Z_vc_h",vtx->getInvMassErr(),weight);
199 Fill1DHisto("vtx_px_vc_h",vtx->getP().X());
200 Fill1DHisto("vtx_py_vc_h",vtx->getP().Y());
201 Fill1DHisto("vtx_pz_vc_h",vtx->getP().Z());
202 Fill1DHisto("vtx_p_vc_h" ,vtx->getP().Mag());
203 Fill1DHisto("vtx_p_std_vc_h" ,vtx->getP().Mag()*(stdBeamEnergy_/eBeam_));
204
205 Fill1DHisto(posTag+"vtx_sigma_X_vc_h",sqrt(vtx->getCovariance()[0]),weight);
206 Fill1DHisto(posTag+"vtx_sigma_Y_vc_h",sqrt(vtx->getCovariance()[3]),weight);
207 Fill1DHisto(posTag+"vtx_sigma_Z_vc_h",sqrt(vtx->getCovariance()[5]),weight);
208 Fill1DHisto(posTag+"vtx_InvM_vc_h" ,vtx->getInvMass(),weight);
209 Fill1DHisto(posTag+"vtx_InvMErr_Z_vc_h",vtx->getInvMassErr(),weight);
210 Fill1DHisto(posTag+"vtx_px_vc_h",vtx->getP().X());
211 Fill1DHisto(posTag+"vtx_py_vc_h",vtx->getP().Y());
212 Fill1DHisto(posTag+"vtx_pz_vc_h",vtx->getP().Z());
213 Fill1DHisto(posTag+"vtx_p_vc_h" ,vtx->getP().Mag());
214 Fill1DHisto(posTag+"vtx_p_std_vc_h" ,vtx->getP().Mag()*(stdBeamEnergy_/eBeam_));
215
216 Fill1DHisto(layerCode+"vtx_sigma_X_vc_h",sqrt(vtx->getCovariance()[0]),weight);
217 Fill1DHisto(layerCode+"vtx_sigma_Y_vc_h",sqrt(vtx->getCovariance()[3]),weight);
218 Fill1DHisto(layerCode+"vtx_sigma_Z_vc_h",sqrt(vtx->getCovariance()[5]),weight);
219 Fill1DHisto(layerCode+"vtx_InvM_vc_h" ,vtx->getInvMass(),weight);
220 Fill1DHisto(layerCode+"vtx_InvMErr_Z_vc_h",vtx->getInvMassErr(),weight);
221 Fill1DHisto(layerCode+"vtx_px_vc_h",vtx->getP().X());
222 Fill1DHisto(layerCode+"vtx_py_vc_h",vtx->getP().Y());
223 Fill1DHisto(layerCode+"vtx_pz_vc_h",vtx->getP().Z());
224 Fill1DHisto(layerCode+"vtx_p_vc_h" ,vtx->getP().Mag());
225 Fill1DHisto(layerCode+"vtx_p_std_vc_h" ,vtx->getP().Mag()*(stdBeamEnergy_/eBeam_));
226
227 Fill1DHisto(layerCode+posTag+"vtx_sigma_X_vc_h",sqrt(vtx->getCovariance()[0]),weight);
228 Fill1DHisto(layerCode+posTag+"vtx_sigma_Y_vc_h",sqrt(vtx->getCovariance()[3]),weight);
229 Fill1DHisto(layerCode+posTag+"vtx_sigma_Z_vc_h",sqrt(vtx->getCovariance()[5]),weight);
230 Fill1DHisto(layerCode+posTag+"vtx_InvM_vc_h" ,vtx->getInvMass(),weight);
231 Fill1DHisto(layerCode+posTag+"vtx_InvMErr_Z_vc_h",vtx->getInvMassErr(),weight);
232 Fill1DHisto(layerCode+posTag+"vtx_px_vc_h",vtx->getP().X());
233 Fill1DHisto(layerCode+posTag+"vtx_py_vc_h",vtx->getP().Y());
234 Fill1DHisto(layerCode+posTag+"vtx_pz_vc_h",vtx->getP().Z());
235 Fill1DHisto(layerCode+posTag+"vtx_p_vc_h" ,vtx->getP().Mag());
236 Fill1DHisto(layerCode+posTag+"vtx_p_std_vc_h" ,vtx->getP().Mag()*(stdBeamEnergy_/eBeam_));
237 TLorentzVector p_ele;
238 //p_ele.SetPxPyPzE(ele->getMomentum()[0], ele->getMomentum()[1],ele->getMomentum()[2],ele->getEnergy());
239 p_ele.SetPxPyPzE(ele_trk->getMomentum()[0],ele_trk->getMomentum()[1],ele_trk->getMomentum()[2],ele->getEnergy());
240
241 TLorentzVector p_pos;
242 //p_pos.SetPxPyPzE(pos->getMomentum()[0], pos->getMomentum()[1],pos->getMomentum()[2],pos->getEnergy());
243 p_pos.SetPxPyPzE(pos_trk->getMomentum()[0],pos_trk->getMomentum()[1],pos_trk->getMomentum()[2],pos->getEnergy());
244
245 //Fill ele and pos information
246 // Fill1DHisto("ele_p_h",p_ele.P(),weight);
247 // Fill1DHisto("pos_p_h",p_pos.P(),weight);
248
249 //Compute some extra variables
250
251 //TODO::Rotate them
252 p_ele.RotateY(-0.0305);
253 p_pos.RotateY(-0.0305);
254
255 //Massless electrons. TODO fix initialization
256 TLorentzVector p_beam(0.,0.,2.3,2.3);
257 TLorentzVector p_v0 = p_ele+p_pos;
258 TLorentzVector p_miss = p_beam - p_v0;
259
260 double thetax_v0_val = TMath::ATan2(p_v0.X(),p_v0.Z());
261 double thetay_v0_val = TMath::ATan2(p_v0.Y(),p_v0.Z());
262
263 double thetax_miss_val = TMath::ATan2(p_miss.X(),p_miss.Z());
264 double thetay_miss_val = TMath::ATan2(p_miss.Y(),p_miss.Z());
265
266 double pt_ele = p_ele.Pt();
267 double pt_pos = p_pos.Pt();
268
269 double pt_asym_val = (pt_ele-pt_pos) / (pt_ele+pt_pos);
270
271
272 //Fill event information
273
274 //1D histos
275 Fill1DHisto("Esum_vc_h",ele->getEnergy() + pos->getEnergy(),weight);
276 Fill1DHisto("Psum_vc_h",p_ele.P() + p_pos.P());
277 Fill1DHisto("Psum_std_vc_h",(p_ele.P() + p_pos.P())*(stdBeamEnergy_/eBeam_));
278 Fill1DHisto("PtAsym_vc_h",pt_asym_val,weight);
279 Fill1DHisto("Pmiss_vc_h",p_miss.P());
280 Fill1DHisto("thetax_v0_vc_h",thetax_v0_val,weight);
281 Fill1DHisto("thetay_v0_vc_h",thetay_v0_val,weight);
282 Fill1DHisto("thetax_miss_vc_h",thetax_miss_val,weight);
283 Fill1DHisto("thetay_miss_vc_h",thetay_miss_val,weight);
284
285
286 //1D histos
287 Fill1DHisto(posTag+"Esum_vc_h",ele->getEnergy() + pos->getEnergy(),weight);
288 Fill1DHisto(posTag+"Psum_vc_h",p_ele.P() + p_pos.P());
289 Fill1DHisto(posTag+"Psum_std_vc_h",(p_ele.P() + p_pos.P())*(stdBeamEnergy_/eBeam_));
290 Fill1DHisto(posTag+"PtAsym_vc_h",pt_asym_val,weight);
291 Fill1DHisto(posTag+"Pmiss_vc_h",p_miss.P());
292 Fill1DHisto(posTag+"thetax_v0_vc_h",thetax_v0_val,weight);
293 Fill1DHisto(posTag+"thetay_v0_vc_h",thetay_v0_val,weight);
294 Fill1DHisto(posTag+"thetax_miss_vc_h",thetax_miss_val,weight);
295 Fill1DHisto(posTag+"thetay_miss_vc_h",thetay_miss_val,weight);
296
297 //1D histos
298 Fill1DHisto(layerCode+"Esum_vc_h",ele->getEnergy() + pos->getEnergy(),weight);
299 Fill1DHisto(layerCode+"Psum_vc_h",p_ele.P() + p_pos.P());
300 Fill1DHisto(layerCode+"Psum_std_vc_h",(p_ele.P() + p_pos.P())*(stdBeamEnergy_/eBeam_));
301 Fill1DHisto(layerCode+"PtAsym_vc_h",pt_asym_val,weight);
302 Fill1DHisto(layerCode+"Pmiss_vc_h",p_miss.P());
303 Fill1DHisto(layerCode+"thetax_v0_vc_h",thetax_v0_val,weight);
304 Fill1DHisto(layerCode+"thetay_v0_vc_h",thetay_v0_val,weight);
305 Fill1DHisto(layerCode+"thetax_miss_vc_h",thetax_miss_val,weight);
306 Fill1DHisto(layerCode+"thetay_miss_vc_h",thetay_miss_val,weight);
307
308
309 //1D histos
310 Fill1DHisto(layerCode+posTag+"Esum_vc_h",ele->getEnergy() + pos->getEnergy(),weight);
311 Fill1DHisto(layerCode+posTag+"Psum_vc_h",p_ele.P() + p_pos.P());
312 Fill1DHisto(layerCode+posTag+"Psum_std_vc_h",(p_ele.P() + p_pos.P())*(stdBeamEnergy_/eBeam_));
313 Fill1DHisto(layerCode+posTag+"PtAsym_vc_h",pt_asym_val,weight);
314 Fill1DHisto(layerCode+posTag+"Pmiss_vc_h",p_miss.P());
315 Fill1DHisto(layerCode+posTag+"thetax_v0_vc_h",thetax_v0_val,weight);
316 Fill1DHisto(layerCode+posTag+"thetay_v0_vc_h",thetay_v0_val,weight);
317 Fill1DHisto(layerCode+posTag+"thetax_miss_vc_h",thetax_miss_val,weight);
318 Fill1DHisto(layerCode+posTag+"thetay_miss_vc_h",thetay_miss_val,weight);
319
320 Fill2DHisto("vtx_InvM_vtx_z_vc_hh",vtx->getInvMass(),vtx->getZ(),weight);
321 Fill2DHisto(posTag+"vtx_InvM_vtx_z_vc_hh",vtx->getInvMass(),vtx->getZ(),weight);
322 Fill2DHisto(layerCode+"vtx_InvM_vtx_z_vc_hh",vtx->getInvMass(),vtx->getZ(),weight);
323 Fill2DHisto(layerCode+posTag+"vtx_InvM_vtx_z_vc_hh",vtx->getInvMass(),vtx->getZ(),weight);
324
325 //
326 Fill2DHisto("ele_vtx_z_vs_z0_over_tanLambda_tc_hh",ele_trk->getZ0()/ele_trk->getTanLambda(), vtx->getZ());
327 Fill2DHisto("pos_vtx_z_vs_z0_over_tanLambda_tc_hh",pos_trk->getZ0()/pos_trk->getTanLambda(), vtx->getZ());
328 Fill2DHisto(eleTag+"vtx_z_vs_z0_over_tanLambda_tc_hh",ele_trk->getZ0()/ele_trk->getTanLambda(), vtx->getZ());
329 Fill2DHisto(posTag+"vtx_z_vs_z0_over_tanLambda_tc_hh",pos_trk->getZ0()/pos_trk->getTanLambda(), vtx->getZ());
330 if(fabs(vtx->getInvMass()-0.0925)<0.010){
331 Fill2DHisto("ele_vtx_z_vs_z0_over_tanLambda_m_eq_92_tc_hh",ele_trk->getZ0()/ele_trk->getTanLambda(), vtx->getZ());
332 Fill2DHisto("pos_vtx_z_vs_z0_over_tanLambda_m_eq_92_tc_hh",pos_trk->getZ0()/pos_trk->getTanLambda(), vtx->getZ());
333 Fill2DHisto(eleTag+"vtx_z_vs_z0_over_tanLambda_m_eq_92_tc_hh",ele_trk->getZ0()/ele_trk->getTanLambda(), vtx->getZ());
334 Fill2DHisto(posTag+"vtx_z_vs_z0_over_tanLambda_m_eq_92_tc_hh",pos_trk->getZ0()/pos_trk->getTanLambda(), vtx->getZ());
335 }
336
337}
338
339
340void TridentHistos::Fill2DTrack(Track* track, float weight, const std::string& trkname) {
341
342
343 if (track) {
344 std::string half="top";
345
346 if(track->getTanLambda()<0)
347 half="bot";
348 std::string tag=half+"_"+trkname;
349 // std::string layerCode="L"+getLayerCodeFromTrack(track)+"_";
350 double d0 = track->getD0();
351 double z0 = track->getZ0();
352 Fill2DHisto(trkname+"tanlambda_vs_phi0_tc_hh",track->getPhi(),track->getTanLambda(), weight);
353 Fill2DHisto(trkname+"d0_vs_p_tc_hh",track->getP(),d0,weight);
354 Fill2DHisto(trkname+"d0_vs_phi0_tc_hh",track->getPhi(),d0,weight);
355 Fill2DHisto(trkname+"d0_vs_tanlambda_tc_hh",track->getTanLambda(),d0,weight);
356
357 Fill2DHisto(trkname+"z0_vs_p_tc_hh",track->getP(),z0,weight);
358 Fill2DHisto(trkname+"z0_vs_phi0_tc_hh",track->getPhi(),z0,weight);
359 Fill2DHisto(trkname+"z0_vs_tanlambda_tc_hh",track->getTanLambda(),z0,weight);
360
361 Fill2DHisto(tag+"tanlambda_vs_phi0_tc_hh",track->getPhi(),track->getTanLambda(), weight);
362 Fill2DHisto(tag+"d0_vs_p_tc_hh",track->getP(),d0,weight);
363 Fill2DHisto(tag+"d0_vs_phi0_tc_hh",track->getPhi(),d0,weight);
364 Fill2DHisto(tag+"d0_vs_tanlambda_tc_hh",track->getTanLambda(),d0,weight);
365
366 Fill2DHisto(tag+"z0_vs_p_tc_hh",track->getP(),z0,weight);
367 Fill2DHisto(tag+"z0_vs_phi0_tc_hh",track->getPhi(),z0,weight);
368 Fill2DHisto(tag+"z0_vs_tanlambda_tc_hh",track->getTanLambda(),z0,weight);
369
370
371 Fill2DHisto(layerCode+trkname+"tanlambda_vs_phi0_tc_hh",track->getPhi(),track->getTanLambda(), weight);
372 Fill2DHisto(layerCode+trkname+"d0_vs_p_tc_hh",track->getP(),d0,weight);
373 Fill2DHisto(layerCode+trkname+"d0_vs_phi0_tc_hh",track->getPhi(),d0,weight);
374 Fill2DHisto(layerCode+trkname+"d0_vs_tanlambda_tc_hh",track->getTanLambda(),d0,weight);
375
376 Fill2DHisto(layerCode+trkname+"z0_vs_p_tc_hh",track->getP(),z0,weight);
377 Fill2DHisto(layerCode+trkname+"z0_vs_phi0_tc_hh",track->getPhi(),z0,weight);
378 Fill2DHisto(layerCode+trkname+"z0_vs_tanlambda_tc_hh",track->getTanLambda(),z0,weight);
379
380 Fill2DHisto(layerCode+tag+"tanlambda_vs_phi0_tc_hh",track->getPhi(),track->getTanLambda(), weight);
381 Fill2DHisto(layerCode+tag+"d0_vs_p_tc_hh",track->getP(),d0,weight);
382 Fill2DHisto(layerCode+tag+"d0_vs_phi0_tc_hh",track->getPhi(),d0,weight);
383 Fill2DHisto(layerCode+tag+"d0_vs_tanlambda_tc_hh",track->getTanLambda(),d0,weight);
384
385 Fill2DHisto(layerCode+tag+"z0_vs_p_tc_hh",track->getP(),z0,weight);
386 Fill2DHisto(layerCode+tag+"z0_vs_phi0_tc_hh",track->getPhi(),z0,weight);
387 Fill2DHisto(layerCode+tag+"z0_vs_tanlambda_tc_hh",track->getTanLambda(),z0,weight);
388
389
390 }
391}
392
393void TridentHistos::Fill1DTrack(Track* track, double trkTimeOffset,float weight, const std::string& trkname) {
394
395 double charge = (double) track->getCharge();
396
397 //2D hits
398 int n_hits_2d = track->getTrackerHitCount();
399
400 TVector3 p_trk;
401 p_trk.SetXYZ(track->getMomentum()[0],track->getMomentum()[1],track->getMomentum()[2]);
402 std::string half="top";
403
404 if(track->getTanLambda()<0)
405 half="bot";
406
407 std::string tag=half+"_"+trkname;
408
409 // std::string layerCode="L"+getLayerCodeFromTrack(track)+"_"+tag;
410 // std::string layerCode="L"+getLayerCodeFromTrack(track)+"_";
411
412 Fill1DHisto(trkname+"p_tc_h",p_trk.Mag(),weight);
413 Fill1DHisto(trkname+"d0_tc_h" ,track->getD0() ,weight);
414 Fill1DHisto(trkname+"Phi_tc_h" ,track->getPhi() ,weight);
415 Fill1DHisto(trkname+"Omega_tc_h" ,track->getOmega() ,weight);
416 Fill1DHisto(trkname+"pT_tc_h" ,charge*track->getPt() ,weight);
417 Fill1DHisto(trkname+"invpT_tc_h" ,charge/track->getPt() ,weight);
418 Fill1DHisto(trkname+"TanLambda_tc_h",track->getTanLambda() ,weight);
419 Fill1DHisto(trkname+"Z0_tc_h" ,track->getZ0() ,weight);
420 Fill1DHisto(trkname+"Z0_over_TanLambda_tc_h",track->getZ0()/track->getTanLambda() ,weight);
421 Fill1DHisto(trkname+"time_tc_h" ,track->getTrackTime()-trkTimeOffset ,weight);
422 Fill1DHisto(trkname+"chi2_tc_h" ,track->getChi2() ,weight);
423 Fill1DHisto(trkname+"chi2ndf_tc_h" ,track->getChi2Ndf() ,weight);
424 Fill1DHisto(trkname+"nShared_tc_h" ,track->getNShared() ,weight);
425 Fill1DHisto(trkname+"nHits_2d_tc_h" ,n_hits_2d ,weight);
426
427
428 Fill1DHisto(tag+"p_tc_h",p_trk.Mag(),weight);
429 Fill1DHisto(tag+"d0_tc_h" ,track->getD0() ,weight);
430 Fill1DHisto(tag+"Phi_tc_h" ,track->getPhi() ,weight);
431 Fill1DHisto(tag+"Omega_tc_h" ,track->getOmega() ,weight);
432 Fill1DHisto(tag+"pT_tc_h" ,-1*charge*track->getPt() ,weight);
433 Fill1DHisto(tag+"invpT_tc_h" ,-1*charge/track->getPt() ,weight);
434 Fill1DHisto(tag+"TanLambda_tc_h",track->getTanLambda() ,weight);
435 Fill1DHisto(tag+"Z0_tc_h" ,track->getZ0() ,weight);
436 Fill1DHisto(tag+"Z0_over_TanLambda_tc_h",track->getZ0()/track->getTanLambda() ,weight);
437 Fill1DHisto(tag+"time_tc_h" ,track->getTrackTime()-trkTimeOffset ,weight);
438 Fill1DHisto(tag+"chi2_tc_h" ,track->getChi2() ,weight);
439 Fill1DHisto(tag+"chi2ndf_tc_h" ,track->getChi2Ndf() ,weight);
440 Fill1DHisto(tag+"nShared_tc_h" ,track->getNShared() ,weight);
441 Fill1DHisto(tag+"nHits_2d_tc_h" ,n_hits_2d ,weight);
442
443 Fill1DHisto(layerCode+trkname+"p_tc_h",p_trk.Mag(),weight);
444 Fill1DHisto(layerCode+trkname+"d0_tc_h" ,track->getD0() ,weight);
445 Fill1DHisto(layerCode+trkname+"Phi_tc_h" ,track->getPhi() ,weight);
446 Fill1DHisto(layerCode+trkname+"Omega_tc_h" ,track->getOmega() ,weight);
447 Fill1DHisto(layerCode+trkname+"pT_tc_h" ,charge*track->getPt() ,weight);
448 Fill1DHisto(layerCode+trkname+"invpT_tc_h" ,charge/track->getPt() ,weight);
449 Fill1DHisto(layerCode+trkname+"TanLambda_tc_h",track->getTanLambda() ,weight);
450 Fill1DHisto(layerCode+trkname+"Z0_tc_h" ,track->getZ0() ,weight);
451 Fill1DHisto(layerCode+trkname+"Z0_over_TanLambda_tc_h",track->getZ0()/track->getTanLambda() ,weight);
452 Fill1DHisto(layerCode+trkname+"time_tc_h" ,track->getTrackTime()-trkTimeOffset ,weight);
453 Fill1DHisto(layerCode+trkname+"chi2_tc_h" ,track->getChi2() ,weight);
454 Fill1DHisto(layerCode+trkname+"chi2ndf_tc_h" ,track->getChi2Ndf() ,weight);
455 Fill1DHisto(layerCode+trkname+"nShared_tc_h" ,track->getNShared() ,weight);
456 Fill1DHisto(layerCode+trkname+"nHits_2d_tc_h" ,n_hits_2d ,weight);
457
458
459 Fill1DHisto(layerCode+tag+"p_tc_h",p_trk.Mag(),weight);
460 Fill1DHisto(layerCode+tag+"d0_tc_h" ,track->getD0() ,weight);
461 Fill1DHisto(layerCode+tag+"Phi_tc_h" ,track->getPhi() ,weight);
462 Fill1DHisto(layerCode+tag+"Omega_tc_h" ,track->getOmega() ,weight);
463 Fill1DHisto(layerCode+tag+"pT_tc_h" ,-1*charge*track->getPt() ,weight);
464 Fill1DHisto(layerCode+tag+"invpT_tc_h" ,-1*charge/track->getPt() ,weight);
465 Fill1DHisto(layerCode+tag+"TanLambda_tc_h",track->getTanLambda() ,weight);
466 Fill1DHisto(layerCode+tag+"Z0_tc_h" ,track->getZ0() ,weight);
467 Fill1DHisto(layerCode+tag+"Z0_over_TanLambda_tc_h",track->getZ0()/track->getTanLambda() ,weight);
468 Fill1DHisto(layerCode+tag+"time_tc_h" ,track->getTrackTime()-trkTimeOffset ,weight);
469 Fill1DHisto(layerCode+tag+"chi2_tc_h" ,track->getChi2() ,weight);
470 Fill1DHisto(layerCode+tag+"chi2ndf_tc_h" ,track->getChi2Ndf() ,weight);
471 Fill1DHisto(layerCode+tag+"nShared_tc_h" ,track->getNShared() ,weight);
472 Fill1DHisto(layerCode+tag+"nHits_2d_tc_h" ,n_hits_2d ,weight);
473
474
475
476 /*
477 for (int ihit=0; ihit<track->getTrackerHitCount();++ihit) {
478 TrackerHit* hit = (TrackerHit*) track->getSvtHits().At(ihit);
479 RawSvtHit* rhit=(RawSvtHit*)(hit->getRawHits()).At(0);
480 int layer=rhit->getLayer();
481 if(layer==4){
482 if(hasLayer4)
483 std::cout<<"What...I already counted layer 4!"<<std::endl;
484 else
485 hasLayer4=true;
486 }
487 */
488
489 bool hasLayer4=false;
490 for (auto & layer : track->getHitLayers()) {
491 if(layer==4){
492 if(hasLayer4)
493 std::cout<<"What...I already counted layer 4!"<<std::endl;
494 else
495 hasLayer4=true;
496 }
497 Fill1DHisto(trkname+"layersHit_tc_h",layer,weight);
498 Fill1DHisto(tag+"layersHit_tc_h",layer,weight);
499 Fill1DHisto(layerCode+trkname+"layersHit_tc_h",layer,weight);
500 Fill1DHisto(layerCode+tag+"layersHit_tc_h",layer,weight);
501 }
502
503}
504
505
506void TridentHistos::Fill1DHistograms(Track *track, Vertex* vtx, float weight ) {
507
508 if (track) {
509 Fill1DTrack(track,-666);
510 }
511
512 //Vertices
513
514 if (vtx) {
515 Fill1DVertex(vtx);
516 }
517}
518
519
520void TridentHistos::Fill1DTrackTruth(Track *track, Track* truth_track, float weight, const std::string& trkname) {
521
522 if (!track || !truth_track)
523 return;
524
525 //Momentum
526 std::vector<double> trk_mom = track->getMomentum();
527 std::vector<double> trk_truth_mom = truth_track->getMomentum();
528
529 double d0 = track->getD0();
530 double d0err = track->getD0Err();
531 double d0_truth = truth_track->getD0();
532 double phi = track->getPhi();
533 double phi_truth = truth_track->getPhi();
534 double phierr = track->getPhiErr();
535 double omega = track->getOmega();
536 double omega_truth = truth_track->getOmega();
537 double omegaerr = track->getOmegaErr();
538 double tanLambda = track->getTanLambda();
539 double tanLambda_truth = truth_track->getTanLambda();
540 double tanLambdaerr = track->getTanLambdaErr();
541 double z0 = track->getZ0();
542 double z0_truth = truth_track->getZ0();
543 double z0err = track->getZ0Err();
544 double p = track->getP();
545 //Charge different wrt Robert's plots.
546 double invPt = -1.*(double) track->getCharge()/track->getPt();
547 double p_truth = truth_track->getP();
548 double invPt_truth = -1*(double) track->getCharge()/truth_track->getPt();
549
550 double diff_percent_invpT = ((invPt - invPt_truth) / invPt_truth) * 100.;
551
552 // truth residuals
553 Fill1DHisto(trkname+"d0_truth_res_h", d0 - d0_truth , weight);
554 Fill1DHisto(trkname+"Phi_truth_res_h", phi - phi_truth , weight);
555 Fill1DHisto(trkname+"Omega_truth_res_h", omega - omega_truth , weight);
556 Fill1DHisto(trkname+"TanLambda_truth_res_h",tanLambda - tanLambda_truth , weight);
557 Fill1DHisto(trkname+"Z0_truth_res_h", z0 - z0_truth , weight);
558 Fill1DHisto(trkname+"p_truth_res_h", p - p_truth , weight);
559 Fill1DHisto(trkname+"invpT_truth_res_h", invPt - invPt_truth , weight);
560 Fill1DHisto(trkname+"invpT_truth_res_percent_h", diff_percent_invpT , weight);
561 Fill1DHisto(trkname+"px_truth_res_h", trk_mom[0] - trk_truth_mom[0] , weight);
562 Fill1DHisto(trkname+"py_truth_res_h", trk_mom[1] - trk_truth_mom[1] , weight);
563 Fill1DHisto(trkname+"pz_truth_res_h", trk_mom[2] - trk_truth_mom[2] , weight);
564
565 // truth pulls
566 Fill1DHisto(trkname+"d0_truth_pull_h", (d0 - d0_truth) / d0err, weight);
567 Fill1DHisto(trkname+"Phi_truth_pull_h", (phi - phi_truth) / phierr , weight);
568 Fill1DHisto(trkname+"Omega_truth_pull_h", (omega - omega_truth) / omegaerr, weight);
569 Fill1DHisto(trkname+"TanLambda_truth_pull_h",(tanLambda - tanLambda_truth) / tanLambdaerr, weight);
570 Fill1DHisto(trkname+"Z0_truth_pull_h", (z0 - z0_truth) / z0err, weight);
571
572}
573
574
575
576void TridentHistos::Fill2DHistograms(Vertex* vtx, float weight) {
577 /*
578 if (vtx) {
579
580 //TODO Improve this.
581 TVector3 vtxPosSvt;
582 vtxPosSvt.SetX(vtx->getX());
583 vtxPosSvt.SetY(vtx->getY());
584 vtxPosSvt.SetZ(vtx->getZ());
585
586 vtxPosSvt.RotateY(-0.0305);
587
588
589 double vtxP = vtx->getP().Mag();
590
591 Fill2DHisto("vtx_InvM_vtx_z_hh",vtx->getInvMass(),vtx->getZ(),weight);
592 Fill2DHisto("vtx_InvM_vtx_svt_z_hh",vtx->getInvMass(),vtxPosSvt.Z(),weight);
593
594 //Fill2DHisto("vtx_p_svt_z_hh",vtxP,vtxPosSvt.Z(),weight);
595 // Fill2DHisto("vtx_p_svt_x_hh",vtxP,vtxPosSvt.X(),weight);
596 //Fill2DHisto("vtx_p_svt_y_hh",vtxP,vtxPosSvt.Y(),weight);
597
598 //Fill2DHisto("vtx_p_sigmaZ_hh",vtxP,vtx->getCovariance()[5],weight);
599 // Fill2DHisto("vtx_p_sigmaX_hh",vtxP,vtx->getCovariance()[3],weight);
600 //Fill2DHisto("vtx_p_sigmaY_hh",vtxP,vtx->getCovariance()[0],weight);
601 }
602 */
603}
604
605void TridentHistos::FillTrackComparisonHistograms(Track* track_x, Track* track_y, float weight) {
606
607 if (doTrkCompPlots) {
608 /*
609 histos2d[m_name+"_d0_vs_d0" ]->Fill(track_x->getD0(),track_y->getD0(),weight);
610 histos2d[m_name+"_Phi_vs_Phi" ]->Fill(track_x->getPhi(),track_y->getPhi(),weight);
611 histos2d[m_name+"_Omega_vs_Omega" ]->Fill(track_x->getOmega(),track_y->getOmega(),weight);
612 histos2d[m_name+"_TanLambda_vs_TanLambda"]->Fill(track_x->getTanLambda(),track_y->getTanLambda(),weight);
613 histos2d[m_name+"_Z0_vs_Z0" ]->Fill(track_x->getZ0(),track_y->getZ0(),weight);
614 histos2d[m_name+"_time_vs_time" ]->Fill(track_x->getTrackTime(),track_y->getTrackTime(),weight);
615 histos2d[m_name+"_chi2_vs_chi2" ]->Fill(track_x->getChi2Ndf(),
616 track_y->getChi2Ndf(),
617 weight);
618 */
619 }
620
621}
622
623
624//Residual Plots ============ They should probably go somewhere else ====================
625
626
627void TridentHistos::FillResidualHistograms(Track* track, int ly, double res, double sigma) {
628
629 double trk_mom = track->getP();
630 std::string lyr = std::to_string(ly);
631
632 TrackerHit* hit = nullptr;
633 //Get the hits on track
634 for (int ihit = 0; ihit<track->getSvtHits().GetEntries();++ihit) {
635 TrackerHit* tmphit = (TrackerHit*) track->getSvtHits().At(ihit);
636 if (tmphit->getLayer() == ly) {
637 hit = tmphit;
638 break;
639 }
640 }
641
642 if (!hit) {
643 std::cout<<"Hit-on-track residual infos not found on hit on track list for ly="<<ly<<std::endl;
644 }
645
646 double hit_y = -9999.;
647 if (hit) {
648 hit_y = hit->getPosition()[1];
649 }
650
651 //General Plots
652 Fill1DHisto("u_res_ly_"+lyr+"_h",res);
653 Fill2DHisto("u_res_ly_"+lyr+"_vsp_hh",trk_mom,res);
654 Fill2DHisto("u_res_ly_"+lyr+"_vsy_hh",hit_y,res);
655
656 //Top = 0 bottom=1 - Per Volume
657 std::string vol = track->getTanLambda()>0 ? "top" : "bot";
658 Fill1DHisto("u_res_ly_"+lyr+"_"+vol+"_h",res);
659 Fill2DHisto("u_res_ly_"+lyr+"_"+vol+"_vsp_hh",trk_mom,res);
660
661}
662
663
664std::pair<CalCluster*, Track*> TridentHistos::getTrackClusterPair(Track* trk,std::vector<CalCluster*>& clusters, float weight){
665
666 double minDelta=666666;
667 double minDeltaX=666666;
668 double minDeltaY=666666;
669 CalCluster* bestCluster=NULL;
670 double minDeltaCut=20.0;
671 double clXMin=666666;
672 double clYMin=666666;
673 for(auto clu:clusters){
674 double clX=(clu->getPosition()).at(0);
675 double clY=(clu->getPosition()).at(1);
676 if(trk->getTanLambda()*clY<0)//trk and cluster in same half
677 continue;
678 int trkCh=trk->getCharge();
679 if(trkCh*clX>0) //positrons on positron side and vice/versa...sign is flipped in track collection?
680 continue;
681 std::vector<double> posAtECal=trk->getPositionAtEcal();
682 double delY=clY-posAtECal.at(1);
683 double delX=clX-posAtECal.at(0);
684 if(abs(delY)<minDelta){//right now, just take delta Y
685 bestCluster=clu;
686 minDelta=abs(delY);
687 minDeltaX=delX;
688 minDeltaY=delY;
689 clXMin=clX;
690 clYMin=clY;
691 }
692 }
693
694 if(minDelta<minDeltaCut){
695 Fill1DHisto("trkClMinDelta_h",minDelta,weight); //should equal deltaY for now
696 Fill1DHisto("trkClDeltaX_h",minDeltaX,weight);
697 Fill1DHisto("trkClDeltaY_h",minDeltaY,weight);
698 Fill2DHisto("trkClDeltaXY_hh",minDeltaX,minDeltaY,weight);
699 Fill2DHisto("trkClDeltaXvsX_hh",clXMin,minDeltaX,weight);
700 Fill2DHisto("trkClDeltaYvsY_hh",clYMin,minDeltaY,weight);
701 Fill2DHisto("trkClDeltaYvsX_hh",clXMin,minDeltaY,weight);
702 Fill2DHisto("trkClDeltaXvsY_hh",clYMin,minDeltaX,weight);
703 }
704
705 return std::pair<CalCluster*,Track*>(bestCluster,trk);
706}
707
708/*
709 * fill cluster/track times and other ecal stuff for both WAB and trident events
710 * mg...5/9/20 currently just do time
711 */
712//void TridentHistos::FillTrackClusterHistos(std::pair<CalCluster, Track*> ele, std::pair<CalCluster, Track*> posOrGamma, double calTimeOffset, double trkTimeOffset,std::vector<CalCluster*> * clusterList, double weight){
713void TridentHistos::FillTrackClusterHistos(std::pair<CalCluster, Track*> ele, std::pair<CalCluster, Track*> posOrGamma, double calTimeOffset, double trkTimeOffset, double weight){
714 CalCluster eleClu=ele.first;
715 Track* eleTrk=ele.second;
716 CalCluster posClu=posOrGamma.first;
717 Track* posTrk=posOrGamma.second; //these "positrons" may be gammas
718
719 std::string half="top";
720 if(eleTrk->getTanLambda()<0)
721 half="bot";
722 std::string eleTag=half+"_ele_";
723
724 half="top";
725 if(posTrk->getTanLambda()<0)
726 half="bot";
727 std::string posTag=half+"_pos_";
728
729 // this does not work, but maybe in new tuple???
730 //std::string eleLayerCode="L"+getLayerCodeFromTrack(&eleTrk)+"_";
731 //std::string posLayerCode="L"+getLayerCodeFromTrack(&posTrk)+"_";
732
733
734 double ele_cluTime=eleClu.getTime()-calTimeOffset;
735 double pos_cluTime=posClu.getTime()-calTimeOffset;
736 double ele_trkTime=eleTrk->getTrackTime()-trkTimeOffset;
737 double pos_trkTime=posTrk->getTrackTime()-trkTimeOffset;
738
739
740 double ele_cluX=eleClu.getPosition().at(0);
741 double ele_cluY=eleClu.getPosition().at(1);
742 double pos_cluX=posClu.getPosition().at(0);
743 double pos_cluY=posClu.getPosition().at(1);
744
746 double ele_trkX=eleTrk->getPositionAtEcal().at(0);
747 double ele_trkY=eleTrk->getPositionAtEcal().at(1);
748 double pos_trkX=posTrk->getPositionAtEcal().at(0);
749 double pos_trkY=posTrk->getPositionAtEcal().at(1);
751
752 double ele_clu_trk_deltaX=ele_cluX-ele_trkX;
753 double pos_clu_trk_deltaX=pos_cluX-pos_trkX;
754
755 double ele_clu_trk_deltaY=ele_cluY-ele_trkY;
756 double pos_clu_trk_deltaY=pos_cluY-pos_trkY;
757
758 // if(eleClu.getTime()>-300){
759 Fill1DHisto("ele_cl_time_tc_h", eleClu.getTime()-calTimeOffset,weight);
760 Fill1DHisto("ele_cl_ene_tc_h",eleClu.getEnergy(),weight);
761 Fill1DHisto("ele_clu_trk_deltaX_tc_h",ele_clu_trk_deltaX,weight);
762 Fill1DHisto("ele_clu_trk_deltaY_tc_h",ele_clu_trk_deltaY,weight);
763
764
765 Fill1DHisto(layerCode+"ele_cl_time_tc_h", eleClu.getTime()-calTimeOffset,weight);
766 Fill1DHisto(layerCode+"ele_cl_ene_tc_h",eleClu.getEnergy(),weight);
767 Fill1DHisto(layerCode+"ele_clu_trk_deltaX_tc_h",ele_clu_trk_deltaX,weight);
768 Fill1DHisto(layerCode+"ele_clu_trk_deltaY_tc_h",ele_clu_trk_deltaY,weight);
769
770 Fill2DHisto("ele_clu_trk_deltaX_vs_cluX_hc_hh",ele_cluX,ele_clu_trk_deltaX,weight);
771 Fill2DHisto("ele_clu_trk_deltaY_vs_cluX_hc_hh",ele_cluX,ele_clu_trk_deltaY,weight);
772 Fill2DHisto("ele_cluY_vs_cluX_hc_hh",ele_cluX,ele_cluY,weight);
773
774 Fill1DHisto(eleTag+"cl_time_tc_h", eleClu.getTime()-calTimeOffset,weight);
775 Fill1DHisto(eleTag+"cl_ene_tc_h",eleClu.getEnergy(),weight);
776 Fill1DHisto(eleTag+"clu_trk_deltaX_tc_h",ele_clu_trk_deltaX,weight);
777 Fill1DHisto(eleTag+"clu_trk_deltaY_tc_h",ele_clu_trk_deltaY,weight);
778
779 Fill1DHisto(layerCode+eleTag+"cl_time_tc_h", eleClu.getTime()-calTimeOffset,weight);
780 Fill1DHisto(layerCode+eleTag+"cl_ene_tc_h",eleClu.getEnergy(),weight);
781 Fill1DHisto(layerCode+eleTag+"clu_trk_deltaX_tc_h",ele_clu_trk_deltaX,weight);
782 Fill1DHisto(layerCode+eleTag+"clu_trk_deltaY_tc_h",ele_clu_trk_deltaY,weight);
783
784 Fill2DHisto(eleTag+"clu_trk_deltaX_vs_cluX_hc_hh",ele_cluX,ele_clu_trk_deltaX,weight);
785 Fill2DHisto(eleTag+"clu_trk_deltaY_vs_cluX_hc_hh",ele_cluX,ele_clu_trk_deltaY,weight);
786 Fill2DHisto(eleTag+"cluY_vs_cluX_hc_hh",ele_cluX,ele_cluY,weight);
787
788 double pele=sqrt(eleTrk->getMomentum()[0]*eleTrk->getMomentum()[0]+
789 eleTrk->getMomentum()[1]*eleTrk->getMomentum()[1]+
790 eleTrk->getMomentum()[2]*eleTrk->getMomentum()[2]);
791 Fill1DHisto("ele_cltrk_time_diff_tc_h", eleClu.getTime()-calTimeOffset-eleTrk->getTrackTime()+trkTimeOffset,weight);
792 Fill1DHisto("ele_EOverp_tc_h",eleClu.getEnergy()/pele,weight);
793 Fill2DHisto("ele_EOverp_vs_cluX_hc_hh",ele_cluX,eleClu.getEnergy()/pele,weight);
794 // if(eleClu.getEnergy()/pele<=0.02){
795 // std::cout<<"ele E/p<=0.02: "<<eleClu.getEnergy()<<" "<<pele<<std::endl;
796 // std::cout<<"electron momentum:: x = "<<eleTrk->getMomentum()[0]<<" y = "<<eleTrk->getMomentum()[1]<<" z = "<<eleTrk->getMomentum()[2]<<std::endl;
797 // std::cout<<"cluster time = "<<ele_cluTime<<" track time = "<<ele_trkTime<<std::endl;
798 //}
799
800 Fill1DHisto(layerCode+"ele_cltrk_time_diff_tc_h", eleClu.getTime()-calTimeOffset-eleTrk->getTrackTime()+trkTimeOffset,weight);
801 Fill1DHisto(layerCode+"ele_EOverp_tc_h",eleClu.getEnergy()/pele,weight);
802
803 Fill1DHisto(eleTag+"cltrk_time_diff_tc_h", eleClu.getTime()-calTimeOffset-eleTrk->getTrackTime()+trkTimeOffset,weight);
804 Fill1DHisto(eleTag+"EOverp_tc_h",eleClu.getEnergy()/pele,weight);
805 Fill2DHisto(eleTag+"EOverp_vs_cluX_hc_hh",ele_cluX,eleClu.getEnergy()/pele,weight);
806
807 Fill1DHisto(layerCode+eleTag+"cltrk_time_diff_tc_h", eleClu.getTime()-calTimeOffset-eleTrk->getTrackTime()+trkTimeOffset,weight);
808 Fill1DHisto(layerCode+eleTag+"EOverp_tc_h",eleClu.getEnergy()/pele,weight);
809
810 Fill1DHisto("pos_cl_time_tc_h", posClu.getTime()-calTimeOffset,weight);
811 Fill1DHisto("pos_cl_ene_tc_h",posClu.getEnergy(),weight);
812 Fill1DHisto("pos_clu_trk_deltaX_tc_h",pos_clu_trk_deltaX,weight);
813 Fill1DHisto("pos_clu_trk_deltaY_tc_h",pos_clu_trk_deltaY,weight);
814 Fill1DHisto("pos_cltrk_time_diff_tc_h", posClu.getTime()-calTimeOffset-posTrk->getTrackTime()+trkTimeOffset,weight);
815
816 Fill1DHisto(posTag+"cl_time_tc_h", posClu.getTime()-calTimeOffset,weight);
817 Fill1DHisto(posTag+"cl_ene_tc_h",posClu.getEnergy(),weight);
818 Fill1DHisto(posTag+"clu_trk_deltaX_tc_h",pos_clu_trk_deltaX,weight);
819 Fill1DHisto(posTag+"clu_trk_deltaY_tc_h",pos_clu_trk_deltaY,weight);
820 Fill1DHisto(posTag+"cltrk_time_diff_tc_h", posClu.getTime()-calTimeOffset-posTrk->getTrackTime()+trkTimeOffset,weight);
821
822 Fill1DHisto(layerCode+"pos_cl_time_tc_h", posClu.getTime()-calTimeOffset,weight);
823 Fill1DHisto(layerCode+"pos_cl_ene_tc_h",posClu.getEnergy(),weight);
824 Fill1DHisto(layerCode+"pos_clu_trk_deltaX_tc_h",pos_clu_trk_deltaX,weight);
825 Fill1DHisto(layerCode+"pos_clu_trk_deltaY_tc_h",pos_clu_trk_deltaY,weight);
826 Fill1DHisto(layerCode+"pos_cltrk_time_diff_tc_h", posClu.getTime()-calTimeOffset-posTrk->getTrackTime()+trkTimeOffset,weight);
827
828 Fill1DHisto(layerCode+posTag+"cl_time_tc_h", posClu.getTime()-calTimeOffset,weight);
829 Fill1DHisto(layerCode+posTag+"cl_ene_tc_h",posClu.getEnergy(),weight);
830 Fill1DHisto(layerCode+posTag+"clu_trk_deltaX_tc_h",pos_clu_trk_deltaX,weight);
831 Fill1DHisto(layerCode+posTag+"clu_trk_deltaY_tc_h",pos_clu_trk_deltaY,weight);
832 Fill1DHisto(layerCode+posTag+"cltrk_time_diff_tc_h", posClu.getTime()-calTimeOffset-posTrk->getTrackTime()+trkTimeOffset,weight);
833
834 double ppos=sqrt(posTrk->getMomentum()[0]*posTrk->getMomentum()[0]+
835 posTrk->getMomentum()[1]*posTrk->getMomentum()[1]+
836 posTrk->getMomentum()[2]*posTrk->getMomentum()[2]);
837 Fill1DHisto("pos_EOverp_tc_h",posClu.getEnergy()/ppos,weight);
838 Fill1DHisto(layerCode+"pos_EOverp_tc_h",posClu.getEnergy()/ppos,weight);
839 Fill2DHisto("pos_EOverp_vs_cluX_hc_hh",pos_cluX,posClu.getEnergy()/ppos,weight);
840 Fill2DHisto("pos_cluY_vs_cluX_hc_hh",pos_cluX,pos_cluY,weight);
841 Fill2DHisto("pos_clu_trk_deltaX_vs_cluX_hc_hh",pos_cluX,pos_clu_trk_deltaX,weight);
842 Fill2DHisto("pos_clu_trk_deltaY_vs_cluX_hc_hh",pos_cluX,pos_clu_trk_deltaY,weight);
843
844 Fill1DHisto(posTag+"EOverp_tc_h",posClu.getEnergy()/ppos,weight);
845 Fill1DHisto(layerCode+posTag+"EOverp_tc_h",posClu.getEnergy()/ppos,weight);
846 Fill2DHisto(posTag+"EOverp_vs_cluX_hc_hh",pos_cluX,posClu.getEnergy()/ppos,weight);
847 Fill2DHisto(posTag+"cluY_vs_cluX_hc_hh",pos_cluX,pos_cluY,weight);
848 Fill2DHisto(posTag+"clu_trk_deltaX_vs_cluX_hc_hh",pos_cluX,pos_clu_trk_deltaX,weight);
849 Fill2DHisto(posTag+"clu_trk_deltaY_vs_cluX_hc_hh",pos_cluX,pos_clu_trk_deltaY,weight);
850 //find the closest cluster that makes sense
851 /*
852 double minDistX=99999;
853 double minDistY=99999;
854 double minDist=99999;
855 int ele_minDistCluIndex=-666;
856 for(int i_clu=0; i_clu<clusterList->size(); i_clu++){
857 CalCluster* clu=clusterList->at(i_clu);
858 double clX=clu->getPosition().at(0);
859 double clY=clu->getPosition().at(1);
860 double clTime=clu->getTime()-calTimeOffset;
861 if(ele_trkY*clY<0)// make sure track and cluster in same half of ecal
862 continue;
863 //add in a cluster-track time cut to get rid of garbage
864 double absClTrkTimeDiff=abs(clTime-ele_trkTime);
865 if(absClTrkTimeDiff>10)
866 continue;
867
868 //find the min distance in
869 double dist=sqrt(pow(ele_trkX-clX,2)+pow(ele_trkY-clY,2));
870 if(dist<minDist){
871 ele_minDistCluIndex=i_clu;
872 }
873 }
874
875 minDistX=99999;
876 minDistY=99999;
877 minDist=99999;
878 int pos_minDistCluIndex=-666;
879 for(int i_clu=0; i_clu<clusterList->size(); i_clu++){
880 CalCluster* clu=clusterList->at(i_clu);
881 double clX=clu->getPosition().at(0);
882 double clY=clu->getPosition().at(1);
883 double clTime=clu->getTime()-calTimeOffset;
884 if(pos_trkY*clY<0)// make sure track and cluster in same half of ecal
885 continue;
886 //add in a cluster-track time cut to get rid of garbage
887 double absClTrkTimeDiff=abs(clTime-pos_trkTime);
888 if(absClTrkTimeDiff>10)
889 continue;
890
891 //find the min distance in
892 double dist=sqrt(pow(pos_trkX-clX,2)+pow(pos_trkY-clY,2));
893 if(dist<minDist){
894 pos_minDistCluIndex=i_clu;
895 }
896 }
897 */
899 /*
900 bool foundOgCluster=false;
901 bool foundNewCluster=false;
902 bool foundSameCluster=false;
903 if(eleClu.getTime()>-300)
904 foundOgCluster=true;
905 if(ele_minDistCluIndex>-1){
906 foundNewCluster=true;
907 CalCluster* ele_myBestClu=clusterList->at(ele_minDistCluIndex);
908 //see if this cluster is the same as cluster found by hps-java
909 double clX=ele_myBestClu->getPosition().at(0);
910 double clY=ele_myBestClu->getPosition().at(1);
911 double clE=ele_myBestClu->getEnergy();
912 double clTime=ele_myBestClu->getTime()-calTimeOffset;;
913 // if(clTime!= ele_cluTime){
914 if(clTime!= ele_cluTime && !foundOgCluster){ //only plot if there was no match
915 //found a different cluster...plot new parameters
916 double new_deltaX=clX-ele_trkX;
917 double new_deltaY=clY-ele_trkY;
918 Fill1DHisto("ele_new_cl_time_tc_h",clTime,weight);
919 Fill1DHisto("ele_new_cl_ene_tc_h",clE,weight);
920 Fill1DHisto("ele_new_clu_trk_deltaX_tc_h",new_deltaX,weight);
921 Fill1DHisto("ele_new_clu_trk_deltaY_tc_h",new_deltaY,weight);
922 Fill1DHisto("ele_new_EOverp_tc_h",clE/pele,weight);
923 */
924 /*
925 Fill1DHisto(layerCode+"ele_new_cl_time_tc_h",clTime,weight);
926 Fill1DHisto(layerCode+"ele_new_cl_ene_tc_h",clE,weight);
927 Fill1DHisto(layerCode+"ele_new_clu_trk_deltaX_tc_h",new_deltaX,weight);
928 Fill1DHisto(layerCode+"ele_new_clu_trk_deltaY_tc_h",new_deltaY,weight);
929 Fill1DHisto(layerCode+"ele_new_EOverp_tc_h",clE/pele,weight);
930 */
931 /*
932 Fill2DHisto("ele_new_clu_trk_deltaX_vs_cluX_hc_hh",clX,new_deltaX,weight);
933 Fill2DHisto("ele_new_clu_trk_deltaY_vs_cluX_hc_hh",clX,new_deltaY,weight);
934 Fill2DHisto("ele_new_EOverp_vs_cluX_hc_hh",clX,clE/pele,weight);
935 Fill2DHisto("ele_new_cluY_vs_cluX_hc_hh",clX,clY,weight);
936
937 Fill1DHisto(eleTag+"new_cl_time_tc_h",clTime,weight);
938 Fill1DHisto(eleTag+"new_cl_ene_tc_h",clE,weight);
939 Fill1DHisto(eleTag+"new_clu_trk_deltaX_tc_h",new_deltaX,weight);
940 Fill1DHisto(eleTag+"new_clu_trk_deltaY_tc_h",new_deltaY,weight);
941 Fill1DHisto(eleTag+"new_EOverp_tc_h",clE/pele,weight);
942 */
943 /*
944 Fill1DHisto(layerCode+eleTag+"new_cl_time_tc_h",clTime,weight);
945 Fill1DHisto(layerCode+eleTag+"new_cl_ene_tc_h",clE,weight);
946 Fill1DHisto(layerCode+eleTag+"new_clu_trk_deltaX_tc_h",new_deltaX,weight);
947 Fill1DHisto(layerCode+eleTag+"new_clu_trk_deltaY_tc_h",new_deltaY,weight);
948 Fill1DHisto(layerCode+eleTag+"new_EOverp_tc_h",clE/pele,weight);
949 */
950 /*
951 Fill2DHisto(eleTag+"new_clu_trk_deltaX_vs_cluX_hc_hh",clX,new_deltaX,weight);
952 Fill2DHisto(eleTag+"new_clu_trk_deltaY_vs_cluX_hc_hh",clX,new_deltaY,weight);
953 Fill2DHisto(eleTag+"new_EOverp_vs_cluX_hc_hh",clX,clE/pele,weight);
954 Fill2DHisto(eleTag+"new_cluY_vs_cluX_hc_hh",clX,clY,weight);
955
956
957 }else{
958 foundSameCluster=true;
959 }
960 }
961 */
962 /*
963 int clMatchCode=-666;
964 if(foundSameCluster)
965 clMatchCode=0;
966 else if(foundOgCluster and !foundNewCluster)
967 clMatchCode=-1;
968 else if(!foundOgCluster and foundNewCluster)
969 clMatchCode=1;
970 else if(foundOgCluster and foundNewCluster)
971 clMatchCode=2;
972 else
973 clMatchCode=-2; // didn't find cluster match with either old or new algorithm
974 Fill1DHisto("ele_new_cluster_match_code_tc_h",clMatchCode,weight);
975 Fill1DHisto(eleTag+"new_cluster_match_code_tc_h",clMatchCode,weight);
976
977 Fill1DHisto(layerCode+"ele_new_cluster_match_code_tc_h",clMatchCode,weight);
978 Fill1DHisto(layerCode+eleTag+"new_cluster_match_code_tc_h",clMatchCode,weight);
979
980
982 foundOgCluster=false;
983 foundNewCluster=false;
984 foundSameCluster=false;
985 if(posClu.getTime()>-300)
986 foundOgCluster=true;
987 if(pos_minDistCluIndex>-1){
988 foundNewCluster=true;
989 CalCluster* pos_myBestClu=clusterList->at(pos_minDistCluIndex);
990 //see if this cluster is the same as cluster found by hps-java
991 double clX=pos_myBestClu->getPosition().at(0);
992 double clY=pos_myBestClu->getPosition().at(1);
993 double clE=pos_myBestClu->getEnergy();
994 double clTime=pos_myBestClu->getTime()-calTimeOffset;;
995 // if(clTime!= pos_cluTime){
996 if(clTime!= pos_cluTime&& !foundOgCluster){ //only plot if there was no match
997 //found a different cluster...plot new parameters
998 double new_deltaX=clX-pos_trkX;
999 double new_deltaY=clY-pos_trkY;
1000 Fill1DHisto("pos_new_cl_time_tc_h",clTime,weight);
1001 Fill1DHisto("pos_new_cl_ene_tc_h",clE,weight);
1002 Fill1DHisto("pos_new_clu_trk_deltaX_tc_h",new_deltaX,weight);
1003 Fill1DHisto("pos_new_clu_trk_deltaY_tc_h",new_deltaY,weight);
1004 Fill1DHisto("pos_new_EOverp_tc_h",clE/ppos,weight);
1005
1006 Fill1DHisto(layerCode+"pos_new_cl_time_tc_h",clTime,weight);
1007 Fill1DHisto(layerCode+"pos_new_cl_ene_tc_h",clE,weight);
1008 Fill1DHisto(layerCode+"pos_new_clu_trk_deltaX_tc_h",new_deltaX,weight);
1009 Fill1DHisto(layerCode+"pos_new_clu_trk_deltaY_tc_h",new_deltaY,weight);
1010 Fill1DHisto(layerCode+"pos_new_EOverp_tc_h",clE/ppos,weight);
1011
1012 Fill2DHisto("pos_new_clu_trk_deltaX_vs_cluX_hc_hh",clX,new_deltaX,weight);
1013 Fill2DHisto("pos_new_clu_trk_deltaY_vs_cluX_hc_hh",clX,new_deltaY,weight);
1014 Fill2DHisto("pos_new_EOverp_vs_cluX_hc_hh",clX,clE/ppos,weight);
1015 Fill2DHisto("pos_new_cluY_vs_cluX_hc_hh",clX,clY,weight);
1016
1017 Fill1DHisto(posTag+"new_cl_time_tc_h",clTime,weight);
1018 Fill1DHisto(posTag+"new_cl_ene_tc_h",clE,weight);
1019 Fill1DHisto(posTag+"new_clu_trk_deltaX_tc_h",new_deltaX,weight);
1020 Fill1DHisto(posTag+"new_clu_trk_deltaY_tc_h",new_deltaY,weight);
1021 Fill1DHisto(posTag+"new_EOverp_tc_h",clE/ppos,weight);
1022
1023 Fill1DHisto(layerCode+posTag+"new_cl_time_tc_h",clTime,weight);
1024 Fill1DHisto(layerCode+posTag+"new_cl_ene_tc_h",clE,weight);
1025 Fill1DHisto(layerCode+posTag+"new_clu_trk_deltaX_tc_h",new_deltaX,weight);
1026 Fill1DHisto(layerCode+posTag+"new_clu_trk_deltaY_tc_h",new_deltaY,weight);
1027 Fill1DHisto(layerCode+posTag+"new_EOverp_tc_h",clE/ppos,weight);
1028
1029 Fill2DHisto(posTag+"new_clu_trk_deltaX_vs_cluX_hc_hh",clX,new_deltaX,weight);
1030 Fill2DHisto(posTag+"new_clu_trk_deltaY_vs_cluX_hc_hh",clX,new_deltaY,weight);
1031 Fill2DHisto(posTag+"new_EOverp_vs_cluX_hc_hh",clX,clE/ppos,weight);
1032 Fill2DHisto(posTag+"new_cluY_vs_cluX_hc_hh",clX,clY,weight);
1033
1034 }else{
1035 foundSameCluster=true;
1036 }
1037 }
1038 */
1039 /*
1040 clMatchCode=-666;
1041 if(foundSameCluster)
1042 clMatchCode=0;
1043 else if(foundOgCluster and !foundNewCluster)
1044 clMatchCode=-1;
1045 else if(!foundOgCluster and foundNewCluster)
1046 clMatchCode=1;
1047 else if(foundOgCluster and foundNewCluster)
1048 clMatchCode=2;
1049 else
1050 clMatchCode=-2; // didn't find cluster match with either old or new algorithm
1051 Fill1DHisto("pos_new_cluster_match_code_tc_h",clMatchCode,weight);
1052 Fill1DHisto(posTag+"new_cluster_match_code_tc_h",clMatchCode,weight);
1053 // Fill1DHisto(layerCode+"pos_new_cluster_match_code_tc_h",clMatchCode,weight);
1054 // Fill1DHisto(layerCode+posTag+"new_cluster_match_code_tc_h",clMatchCode,weight);
1055 */
1056 Fill1DHisto("cl_time_diff_vc_h", posClu.getTime()-eleClu.getTime(),weight);
1057 Fill1DHisto("trk_time_diff_vc_h", posTrk->getTrackTime()-eleTrk->getTrackTime(),weight);
1058
1059 Fill1DHisto(posTag+"cl_time_diff_vc_h", posClu.getTime()-eleClu.getTime(),weight);
1060 Fill1DHisto(posTag+"trk_time_diff_vc_h", posTrk->getTrackTime()-eleTrk->getTrackTime(),weight);
1061
1062
1063 // Fill1DHisto("gamma_cl_time_h", posClu.getTime()-calTimeOffset,weight);
1064 //Fill1DHisto("gamma_cl_ene_h",posClu.getEnergy(),weight);
1065
1066}
1067
1068/*
1069 * fill WAB specific histos
1070 * already have track and cluster plots from above methods
1071 */
1072void TridentHistos::FillWABHistos(std::pair<CalCluster*, Track*> ele, CalCluster* gamma,double weight){
1073 CalCluster* eleClu=ele.first;
1074 Track* eleTrk=ele.second;
1075 double esum=0;
1076 if(eleClu)
1077 esum=eleClu->getEnergy()+gamma->getEnergy();
1078 double psum=sqrt(eleTrk->getMomentum()[0]*eleTrk->getMomentum()[0]+
1079 eleTrk->getMomentum()[1]*eleTrk->getMomentum()[1]+
1080 eleTrk->getMomentum()[2]*eleTrk->getMomentum()[2])+ gamma->getEnergy();
1081 if(eleClu)
1082 Fill1DHisto("wab_esum_h",esum,weight);
1083 Fill1DHisto("wab_psum_h",psum,weight);
1084}
1085
1086void TridentHistos::Fill1DVertex(Vertex* vtx, float weight){
1087 //kept for now...
1088}
1089
1090
1091/*
1092 * overwrite the saveHistos routine to make layer-based directories....
1093 */
1094
1095void TridentHistos::saveHistos(TFile* outF,std::string folder) {
1096 if (outF) outF->cd();
1097 TDirectory* dir{nullptr};
1098 if (!folder.empty()) {
1099 dir = outF->mkdir(folder.c_str());
1100 dir->cd();
1101 }
1102
1103 std::regex pat;
1104
1105 for (it3d it = histos3d.begin(); it!=histos3d.end(); ++it) {
1106 if (!it->second){
1107 continue;
1108 }
1109 setOutputDir(outF,folder,it->first);
1110 it->second->Write();
1111 }
1112 for (it2d it = histos2d.begin(); it!=histos2d.end(); ++it) {
1113 if (!(it->second)) {
1114 continue;
1115 }
1116 setOutputDir(outF,folder,it->first);
1117 it->second->Write();
1118 }
1119
1120 for (it1d it = histos1d.begin(); it!=histos1d.end(); ++it) {
1121 if (!it->second){
1122 continue;
1123 }
1124 setOutputDir(outF,folder,it->first);
1125 it->second->Write();
1126 }
1127
1128 //dir->Write();
1129 //if (dir) {delete dir; dir=0;}
1130
1131 Clear();
1132
1133}
1134
1135/*
1136 * Figure out if we want this histo in a layer-seperated directory, find layer code,
1137 * make directory if needed, an cd to appropriate dir
1138 */
1139void TridentHistos::setOutputDir(TFile* outF, std::string folder, std::string histoName){
1140 std::string fullDirPath=folder;
1141 std::string layerCode="";
1142 //first, check if this histo was split by both ele/pos layers (e.g. for V0 quantities)
1143 std::string dualLayerCode=splitByElePosLayerCombos(histoName);
1144 if(dualLayerCode.length()>0)
1145 layerCode=dualLayerCode;
1146 else
1147 layerCode=getLayerCodeFromHistoName(histoName);
1148 if(layerCode.length()>0){
1149 fullDirPath=folder+"/"+layerCode;
1150 if(outF->GetDirectory(fullDirPath.c_str())==0){
1151 TDirectory* ldir{nullptr};
1152 ldir=outF->mkdir(std::string(fullDirPath).c_str());
1153 }
1154 }
1155 outF->GetDirectory(fullDirPath.c_str())->cd();
1156}
1157
1158
1159/*
1160 *check if the histo name has a layer code in it; if it does, return the layer code
1161 *else return empty string
1162 *....the easy way to do this is to use regex but it's broken for GCC 4.8.5??
1163 */
1164
1165std::string TridentHistos::getLayerCodeFromHistoName(std::string histoName){
1166 std::string layerCode="";
1167 int nLayers=4;
1168 int padding=1;//the "L"
1169 for (int i=0; i<histoName.length()-(padding+nLayers); i++)
1170 if(isLayerCode(histoName, nLayers,i)){
1171 layerCode=histoName.substr(i,padding+nLayers);
1172 break;
1173 }
1174
1175 return layerCode;
1176}
1177
1178bool TridentHistos::isLayerCode(std::string histoName, int nLayers, int ptr ){
1179 if(histoName.substr(ptr,1)!="L")//my code starts with capital L
1180 return false;
1181 for(int i=0; i<nLayers;i++){
1182 if(! (histoName.substr(ptr+i+1,1)=="0" ||histoName.substr(ptr+i+1,1)=="1"))
1183 return false;
1184 }
1185 //if we make it this far, it must be true
1186 return true;
1187
1188}
1189
1190
1191std::string TridentHistos::getLayerCodeFromTrack(Track* trk){
1192 std::string s1="0";
1193 std::string s2="0";
1194 std::string s3="0";
1195 std::string s4="0";
1196 for (auto & layer: trk->getHitLayers()) {
1197 // TrackerHit* hit = (TrackerHit*) trk->getSvtHits().At(ihit);
1198 //RawSvtHit* rhit=(RawSvtHit*)(hit->getRawHits()).At(0);
1199 //int layer=rhit->getLayer();
1200 // if(layer == 0 ){
1201 //std::cout<<"I didn't think you could have layer 0???"<<std::endl;
1202 //}
1203 if (layer == 0 ) {
1204 s1="1";
1205 }
1206 if (layer == 1) {
1207 s2="1";
1208 }
1209 if (layer == 2) {
1210 s3="1";
1211 }
1212 if (layer == 3) {
1213 s4="1";
1214 }
1215 }
1216
1217 return std::string(s1+s2+s3+s4);
1218
1219}
1220
1221/*
1222 * look for the electron/positron layer codes for splitting V0-related quantities
1223 * code should look like posXXXX_eleXXXX
1224 */
1225std::string TridentHistos::splitByElePosLayerCombos(std::string histoName){
1226 std::string layerCode="";
1227 int nLayers=4;
1228 int padding=3;//the "ele" or "pos"
1229 std::string firstSt="pos";
1230 std::string secSt="ele";
1231 bool foundFirst=false;
1232 for (int i=0; i<histoName.length()-(padding+nLayers); i++){
1233 foundFirst=false;
1234 layerCode="";
1235 if(histoName.substr(i,padding)==firstSt && (histoName.substr(i+padding,1)=="0"||histoName.substr(i+padding,1)=="1")){
1236 foundFirst=true;
1237 layerCode=histoName.substr(i,padding+nLayers);
1238 }
1239 if(foundFirst){
1240 if(histoName.substr(i+padding+nLayers+1,padding)==secSt && (histoName.substr(i+2*padding+nLayers+1,1)=="0"||histoName.substr(i+2*padding+nLayers+1,1)=="1")){
1241 layerCode+="_"+histoName.substr(i+padding+nLayers+1,padding+nLayers);
1242 return layerCode;
1243 }
1244 }
1245 }
1246 return layerCode;
1247
1248}
1249
1250/*************************************************
1251*
1252* Makes histograms based on template in json file
1253* like DefineHistos but _only_ makes the histos
1254* that are split by the makeCopyJsonTag.
1255*
1256*************************************************/
1257
1258void TridentHistos::DefineHistosFromTemplateOnly(std::vector<std::string> histoCopyNames, std::string makeCopyJsonTag){
1259 if (debug_ ) std::cout << "[HistoManager] DefineHistosFromTemplateOnly" << std::endl;
1260 std::string h_name = "";
1261 if (debug_ ){
1262 for (auto hist : _h_configs.items()){
1263 std::cout<<hist.key()<<std::endl;
1264 }
1265 }
1266 for (auto hist : _h_configs.items()) {
1267 bool singleCopy = true;
1268 for(int i = 0; i < histoCopyNames.size(); i++){
1269 h_name = m_name+"_" + hist.key() ;
1270 std::size_t found = (hist.key()).find_last_of("_");
1271 std::string extension = hist.key().substr(found+1);
1272 if (histoCopyNames.size() > 1 && std::string(hist.key()).find(makeCopyJsonTag) != std::string::npos){
1273 h_name = m_name+"_"+ histoCopyNames.at(i) + "_" + hist.key() ;
1274 if (debug_ )std::cout <<"DefineHistosFromTemplateOnly: base hist name: "<<m_name<<"_"<< hist.key()<<" copy tag: "<<makeCopyJsonTag<<" hist copy list size " << histoCopyNames.size() << std::endl;
1275 if (debug_ )std::cout<<"new name = "<<h_name<<std::endl;
1276
1277 if(histoCopyNames.at(i)=="")
1278 h_name = m_name+"_"+ hist.key();//this lets you use and empty string and not screw up spacing
1279 singleCopy = false;
1280 }
1282 if(singleCopy)
1283 break;
1284
1285 if(debug_){
1286 std::cout << "DefineHistosFromTemplateOnly: " << h_name << std::endl;
1287 std::cout << extension << hist.value().at("xtitle") << std::endl;
1288 }
1289 if (extension == "h") {
1290 histos1d[h_name] = plot1D(h_name,hist.value().at("xtitle"),
1291 hist.value().at("bins"),
1292 hist.value().at("minX"),
1293 hist.value().at("maxX"));
1294
1295 std::string ytitle = hist.value().at("ytitle");
1296
1297 histos1d[h_name]->GetYaxis()->SetTitle(ytitle.c_str());
1298
1299 if (hist.value().contains("labels")) {
1300 std::vector<std::string> labels = hist.value().at("labels").get<std::vector<std::string> >();
1301
1302 if (labels.size() < hist.value().at("bins")) {
1303 std::cout<<"Cannot apply labels to histogram:"<<h_name<<std::endl;
1304 }
1305 else {
1306 for (int i = 1; i<=hist.value().at("bins");++i)
1307 histos1d[h_name]->GetXaxis()->SetBinLabel(i,labels[i-1].c_str());
1308 }//bins
1309 }//labels
1310 }//1D histo
1311
1312 else if (extension == "hh") {
1313 histos2d[h_name] = plot2D(h_name,
1314 hist.value().at("xtitle"),hist.value().at("binsX"),hist.value().at("minX"),hist.value().at("maxX"),
1315 hist.value().at("ytitle"),hist.value().at("binsY"),hist.value().at("minY"),hist.value().at("maxY"));
1316 }
1317
1318
1319 }
1320 }//loop on config
1321 if(debug_)std::cout<<"DefineHistosFromTemplateOnly: Done with copies"<<std::endl;
1322}
1323
1324/*
1325 * call this after you've made all of the split histograms via DefineHistosFromTemplateOnly
1326 * to make the histograms with no split
1327 */
1328void TridentHistos::DefineOneTimeHistos(){
1329 if(debug_)std::cout << "[HistoManager] DefineOneTimeHistos" << std::endl;
1330 std::string h_name = "";
1331
1332 for (auto hist : _h_configs.items()) {
1333
1334 h_name =hist.key() ;
1335 bool foundMatch=false;
1336 //check 1d histo list
1337 for (auto i=histos1d.begin(); i!=histos1d.end(); i++){
1338 if( std::string(i->first).find(h_name) != std::string::npos && std::string(i->first).find(m_name)!= std::string::npos){
1339 foundMatch=true;
1340 }
1341 }
1342 if(foundMatch)
1343 continue;
1344 //check 2d histo list
1345 for (auto i=histos2d.begin(); i!=histos2d.end(); i++){
1346 if( std::string(i->first).find(h_name) != std::string::npos && std::string(i->first).find(m_name)!= std::string::npos){
1347 foundMatch=true;
1348 }
1349 }
1350
1351 if(foundMatch)
1352 continue;
1353
1354 //if we get here, haven't found any existing histos with a name like stub=hname
1355 //so, make it.
1356 h_name=m_name+"_"+hist.key();
1357 std::size_t found = (hist.key()).find_last_of("_");
1358 std::string extension = hist.key().substr(found+1);
1359 if(debug_){
1360 std::cout << "DefineOneTimeHistos: " << h_name << std::endl;
1361 std::cout << extension << hist.value().at("xtitle") << std::endl;
1362 }
1363 if (extension == "h") {
1364 histos1d[h_name] = plot1D(h_name,hist.value().at("xtitle"),
1365 hist.value().at("bins"),
1366 hist.value().at("minX"),
1367 hist.value().at("maxX"));
1368
1369 std::string ytitle = hist.value().at("ytitle");
1370
1371 histos1d[h_name]->GetYaxis()->SetTitle(ytitle.c_str());
1372
1373 if (hist.value().contains("labels")) {
1374 std::vector<std::string> labels = hist.value().at("labels").get<std::vector<std::string> >();
1375
1376 if (labels.size() < hist.value().at("bins")) {
1377 std::cout<<"Cannot apply labels to histogram:"<<h_name<<std::endl;
1378 }
1379 else {
1380 for (int i = 1; i<=hist.value().at("bins");++i)
1381 histos1d[h_name]->GetXaxis()->SetBinLabel(i,labels[i-1].c_str());
1382 }//bins
1383 }//labels
1384 }
1385 else if (extension == "hh") {
1386 histos2d[h_name] = plot2D(h_name,
1387 hist.value().at("xtitle"),hist.value().at("binsX"),hist.value().at("minX"),hist.value().at("maxX"),
1388 hist.value().at("ytitle"),hist.value().at("binsY"),hist.value().at("minY"),hist.value().at("maxY"));
1389 }
1390
1391
1392 }
1393 if(debug_)std::cout<<"DefineOneTimeHistos: Done with copies"<<std::endl;
1394}
1395
CalCluster::getEnergy
double getEnergy() const
Definition CalCluster.h:73
CalCluster::getPosition
std::vector< double > getPosition() const
Definition CalCluster.h:63
CalCluster::getTime
double getTime() const
Definition CalCluster.h:83
HistoManager::Fill2DHisto
void Fill2DHisto(const std::string &histoName, float valuex, float valuey, float weight=1.)
description
Definition HistoManager.cxx:328
HistoManager::histos2d
std::map< std::string, TH2F * > histos2d
description
Definition HistoManager.h:335
HistoManager::plot2D
TH2F * plot2D(std::string name, std::string xtitle, int nbinsX, float xmin, float xmax, std::string ytitle, int nbinsY, float ymin, float ymax)
description
Definition HistoManager.cxx:203
HistoManager::it3d
std::map< std::string, TH3F * >::iterator it3d
description
Definition HistoManager.h:338
HistoManager::it1d
std::map< std::string, TH1F * >::iterator it1d
description
Definition HistoManager.h:334
HistoManager::histos3d
std::map< std::string, TH3F * > histos3d
description
Definition HistoManager.h:337
HistoManager::histos1d
std::map< std::string, TH1F * > histos1d
description
Definition HistoManager.h:333
HistoManager::it2d
std::map< std::string, TH2F * >::iterator it2d
description
Definition HistoManager.h:336
HistoManager::Clear
virtual void Clear()
description
Definition HistoManager.cxx:18
HistoManager::Fill1DHisto
void Fill1DHisto(const std::string &histoName, float value, float weight=1.)
description
Definition HistoManager.cxx:363
HistoManager::plot1D
TH1F * plot1D(const std::string &name, const std::string &xtitle, int nbinsX, float xmin, float xmax)
description
Definition HistoManager.cxx:189
HistoManager::_h_configs
json _h_configs
description
Definition HistoManager.h:341
HistoManager::m_name
std::string m_name
description
Definition HistoManager.h:331
HistoManager::debug_
bool debug_
description
Definition HistoManager.h:340
Particle::getEnergy
double getEnergy() const
Definition Particle.h:143
Track
Definition Track.h:32
Track::getHitLayers
std::vector< int > getHitLayers()
Definition Track.h:99
Track::getPhi
double getPhi() const
Definition Track.h:88
Track::getOmega
double getOmega() const
Definition Track.h:90
Track::getTrackTime
double getTrackTime() const
Definition Track.h:160
Track::getPositionAtEcal
std::vector< double > getPositionAtEcal()
Definition Track.cxx:53
Track::getD0
double getD0() const
Definition Track.h:87
Track::getCharge
int getCharge() const
Definition Track.h:247
Track::getTanLambda
double getTanLambda() const
Definition Track.h:91
Track::getMomentum
std::vector< double > getMomentum()
Definition Track.h:285
Track::getZ0
double getZ0() const
Definition Track.h:92
Track::getPt
double getPt() const
Definition Track.h:293
Track::getP
double getP() const
Definition Track.h:291
TrackerHit::getLayer
int getLayer() const
Definition TrackerHit.h:109
TrackerHit::getPosition
std::vector< double > getPosition() const
Definition TrackerHit.h:47
TridentHistos::getTrackClusterPair
std::pair< CalCluster *, Track * > getTrackClusterPair(Track *trk, std::vector< CalCluster * > &clusters, float weight)
Definition TridentHistos.cxx:664
TridentHistos::DefineHistos
void DefineHistos()
Definition of histograms from json config.
Definition TridentHistos.cxx:9
TridentHistos::Fill1DTrack
void Fill1DTrack(Track *track, double trkTimeOffset, float weight=1., const std::string &trkname="")
Definition TridentHistos.cxx:393
TridentHistos::saveHistos
void saveHistos(TFile *outF, std::string folder)
save histograms
Definition TridentHistos.cxx:1095
TridentHistos::FillTrackComparisonHistograms
void FillTrackComparisonHistograms(Track *track_x, Track *track_y, float weight=1.)
Definition TridentHistos.cxx:605
TridentHistos::layerCode
std::string layerCode
Definition TridentHistos.h:67
TridentHistos::Fill2DTrack
void Fill2DTrack(Track *track, float weight=1., const std::string &trkname="")
Definition TridentHistos.cxx:340
TridentHistos::stdBeamEnergy_
double stdBeamEnergy_
Definition TridentHistos.h:76
TridentHistos::FillResidualHistograms
void FillResidualHistograms(Track *track, int ly, double res, double sigma)
Definition TridentHistos.cxx:627
TridentHistos::DefineHistosFromTemplateOnly
virtual void DefineHistosFromTemplateOnly(std::vector< std::string > histoCopyNames, std::string makeCopyJsonTag="default=single_copy")
Definition TridentHistos.cxx:1258
TridentHistos::Fill1DHistograms
void Fill1DHistograms(Track *track=nullptr, Vertex *vtx=nullptr, float weight=1.)
Definition TridentHistos.cxx:506
TridentHistos::Fill1DVertex
void Fill1DVertex(Vertex *vtx, float weight=1.)
Definition TridentHistos.cxx:1086
TridentHistos::FillWABHistos
void FillWABHistos(std::pair< CalCluster *, Track * > ele, CalCluster *gamma, double weight)
Definition TridentHistos.cxx:1072
TridentHistos::Define2DHistos
virtual void Define2DHistos()
define additional 2D histo by hand
Definition TridentHistos.cxx:124
TridentHistos::Fill1DTrackTruth
void Fill1DTrackTruth(Track *track, Track *truth_track, float weight=1., const std::string &="")
Definition TridentHistos.cxx:520
TridentHistos::DefineOneTimeHistos
virtual void DefineOneTimeHistos()
Definition TridentHistos.cxx:1328
TridentHistos::getLayerCodeFromTrack
std::string getLayerCodeFromTrack(Track *trk)
Definition TridentHistos.cxx:1191
TridentHistos::Fill2DHistograms
void Fill2DHistograms(Vertex *vtx=nullptr, float weight=1.)
Definition TridentHistos.cxx:576
TridentHistos::FillTrackClusterHistos
void FillTrackClusterHistos(std::pair< CalCluster, Track * > ele, std::pair< CalCluster, Track * > posOrGamma, double calTimeOffset, double trkTimeOffset, double weight)
Definition TridentHistos.cxx:713
TridentHistos::setOutputDir
void setOutputDir(TFile *outF, std::string folder, std::string histoName)
Definition TridentHistos.cxx:1139
TridentHistos::AssignLayerCode
void AssignLayerCode(Track *ele_trk, Track *pos_trk)
Definition TridentHistos.cxx:128
TridentHistos::isLayerCode
bool isLayerCode(std::string histoName, int nLayers, int ptr)
Definition TridentHistos.cxx:1178
TridentHistos::splitByElePosLayerCombos
std::string splitByElePosLayerCombos(std::string histoName)
Definition TridentHistos.cxx:1225
TridentHistos::getLayerCodeFromHistoName
std::string getLayerCodeFromHistoName(std::string name)
Definition TridentHistos.cxx:1165
Generated by doxygen 1.9.8