-
Notifications
You must be signed in to change notification settings - Fork 1
/
holdcon
142 lines (123 loc) · 6.51 KB
/
holdcon
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
void MCSAnalysis::ConvolveWithInputDistribution(std::string distname){
int isfirst = 0;
if (distname.find(modelname1.c_str()) != std::string::npos)
isfirst = 1;
TFile* infile = new TFile(modelfile.c_str());
std::string tmpname = "thetaX_refconv_";
tmpname += distname;
TH1D* thetaX_refconv =
new TH1D(tmpname.c_str(),"Change in Projected Angle (X);#Delta#theta_{X}; Events per 4 mrad",
_histlimits["NbinsXY"], _histlimits["minXY"], _histlimits["maxXY"]);
tmpname = "thetaY_refconv_";
tmpname += distname;
TH1D* thetaY_refconv =
new TH1D(tmpname.c_str(),"Change in Projected Angle (Y);#Delta#theta_{Y}; Events per 4 mrad",
_histlimits["NbinsXY"], _histlimits["minXY"], _histlimits["maxXY"]);
tmpname = "thetaScatt_refconv_";
tmpname += distname;
TH1D* thetaScatt_refconv =
new TH1D(tmpname.c_str(),"Scattering Angle between Momentum Vectors;#theta_{Scatt}; Events per mrad",
_histlimits["NbinsTh"], _histlimits["minTh"], _histlimits["maxTh"]);
tmpname = "theta2Scatt_refconv_";
tmpname += distname;
TH1D* theta2Scatt_refconv =
new TH1D(tmpname.c_str(),"Scattering Angle between Momentum Vectors;#theta^{2}_{Scatt}; Events per mrad",
_histlimits["NbinsTh2"], _histlimits["minTh2"], _histlimits["maxTh2"]);
tmpname = "thetaScatt_refconv_vp_";
tmpname += distname;
TH2D* thetaScatt_refconv_vp =
new TH2D(tmpname.c_str(),"Scattering Angle between Momentum Vectors;Momentum (MeV/c); #theta_{Scatt}",
200, 100, 300, _histlimits["NbinsTh"], _histlimits["minTh"], _histlimits["maxTh"]);
TH2D* thetaXUS_thetaXDS =
new TH2D("thetaXUS_thetaXDS","Upstream vs. Downstream Angle;#theta_{X}^{US}; #theta_{X}^{DS}",
_histlimits["NbinsXY"], _histlimits["minXY"], _histlimits["maxXY"],
_histlimits["NbinsXY"], _histlimits["minXY"], _histlimits["maxXY"]);
TH2D* thetaYUS_thetaYDS =
new TH2D("thetaYUS_thetaYDS","Upstream vs. Downstream Angle;#theta_{X}^{US}; #theta_{X}^{DS}",
_histlimits["NbinsXY"], _histlimits["minXY"], _histlimits["maxXY"],
_histlimits["NbinsXY"], _histlimits["minXY"], _histlimits["maxXY"]);
tmpname = "thetaX_";
tmpname += distname;
std::cout<<"Convolution with "<<tmpname<<" from "<<modelfile<<std::endl;
TH1D* hx = (TH1D*)infile->Get(tmpname.c_str());
TH1D* hy = (TH1D*)infile->Get(tmpname.c_str());
// Collection DSConvSet;
// for (int l=0; l<10; l++){
for (int i=0; i<_USMCset.N(); i++){
for (int j=0; j<20; j++){
double dthetaX = hx->GetRandom() * _sys["resX"];
double dthetaY = hy->GetRandom() * _sys["resY"];
// First project the upstream track to the absorber
double zabspos = _sys["abspos"] + 0.0;
Vars projvarAbs = PropagateVarsMu(_USMCset.E(i), zabspos);
double xabs = projvarAbs.X; /// _USMCset.E(i).X + _USMCset.E(i).dXdz * dzabsUS;
double yabs = projvarAbs.Y; /// _USMCset.E(i).Y + _USMCset.E(i).dYdz * dzabsUS;
// Now add the angle from the model to the downstream measurement.
double dXdz_abs = _DSMCset.E(i).dXdz + tan(dthetaY);
double dYdz_abs = _DSMCset.E(i).dYdz + tan(dthetaX);
// double d_thetaY = atan(dXdz_abs) - atan(_DSMCset.E(i).dXdz);
// double d_thetaX = atan(dYdz_abs) - atan(_DSMCset.E(i).dYdz);
// Project the track into the downstream reference plane
// double xref = xabs + dXdz_abs * dzabsDS;
// double yref = yabs + dYdz_abs * dzabsDS;
Vars tmpvar = _DSMCset.E(i);
tmpvar.X = xabs;
tmpvar.Y = yabs;
tmpvar.Z = _sys["abspos"];
tmpvar.dXdz = dXdz_abs;
tmpvar.dYdz = dYdz_abs;
tmpvar.px = dXdz_abs * _USMCset.E(i).pz;
tmpvar.py = dYdz_abs * _USMCset.E(i).pz;
tmpvar.pz = _USMCset.E(i).pz;
tmpvar.TOF12= _USMCset.E(i).TOF12;
tmpvar.TOF01= _USMCset.E(i).TOF01;
std::vector<double> projDTheta = DefineProjectionAngles(tmpvar, _DSMCset.E(i));
double d_thetaX = projDTheta[0];
double d_thetaY = projDTheta[1];
//double cosdthetaScatt = ( (1 + dXdz_abs * _DSMCset.E(i).dXdz + dYdz_abs*_DSMCset.E(i).dYdz) /
// sqrt( 1 + dXdz_abs*dXdz_abs + dYdz_abs*dYdz_abs) /
// sqrt(1 + _DSMCset.E(i).dXdz*_DSMCset.E(i).dXdz + _DSMCset.E(i).dYdz * _DSMCset.E(i).dYdz) );
double dthetaScatt = projDTheta[2]; /// acos(cosdthetaScatt);
Vars projvar = PropagateVarsMu(tmpvar, _sys["abspos"] + 2993.05);
// Remove events that do not pass through the tracker at its centre
// double xtracker = xabs + dXdz_abs * (_sys["abspos"] + 549.95);
// double ytracker = yabs + dYdz_abs * (dzabsDS + 549.95);
if( sqrt(projvar.X*projvar.X + projvar.Y*projvar.Y) > meanp ){
tmpvar.X = 0.0;
tmpvar.Y = 0.0;
tmpvar.Z = 0.0;
tmpvar.dXdz = 1./2.;
tmpvar.dYdz = 1./2.;
}
std::vector<double> projTheta = DefineProjectionAngles(_USMCset.E(i), tmpvar);
double thetaY = projTheta[1]; /// atan(tmpvar.dXdz) - atan(_USMCset.E(i).dXdz);
double thetaX = projTheta[0]; /// atan(tmpvar.dYdz) - atan(_USMCset.E(i).dYdz);
// double cosScatt = ( (1 + _USMCset.E(i).dXdz * tmpvar.dXdz +
// _USMCset.E(i).dYdz * tmpvar.dYdz )/
// sqrt(1 + _USMCset.E(i).dXdz*_USMCset.E(i).dXdz +
// _USMCset.E(i).dYdz*_USMCset.E(i).dYdz)/
// sqrt(1 + tmpvar.dXdz*tmpvar.dXdz +
// tmpvar.dYdz*tmpvar.dYdz));
double thetaScatt = projTheta[2]; /// acos(cosScatt);
thetaXUS_thetaXDS->Fill(atan(_USMCset.E(i).dXdz), atan(tmpvar.dXdz));
thetaYUS_thetaYDS->Fill(atan(_USMCset.E(i).dYdz), atan(tmpvar.dYdz));
thetaX_refconv->Fill(thetaX);
thetaY_refconv->Fill(thetaY);
thetaScatt_refconv->Fill(thetaScatt);
theta2Scatt_refconv->Fill(thetaScatt*thetaScatt);
thetaScatt_refconv_vp->Fill(_USMCset.E(i).pz, thetaScatt);
isfirst == 1 ? tresp_thetaX.Fill(thetaX, d_thetaX) : resp_thetaX.Fill(thetaX, d_thetaX);
isfirst == 1 ? tresp_thetaY.Fill(thetaY, d_thetaY) : resp_thetaY.Fill(thetaY, d_thetaY);
isfirst == 1 ? tresp_thetaScatt.Fill(thetaScatt, dthetaScatt) : resp_thetaScatt.Fill(thetaScatt, dthetaScatt);
isfirst == 1 ? tresp_theta2Scatt.Fill(thetaScatt*thetaScatt, dthetaScatt*dthetaScatt) : resp_theta2Scatt.Fill(thetaScatt*thetaScatt, dthetaScatt*dthetaScatt);
}
}
outfile->cd();
thetaXUS_thetaXDS->Write();
thetaYUS_thetaYDS->Write();
thetaX_refconv->Write();
thetaY_refconv->Write();
thetaScatt_refconv->Write();
theta2Scatt_refconv->Write();
thetaScatt_refconv_vp->Write();
}