Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AnalyzeResiduals.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file AnalyzeResiduals.C
1 #include <trackbase/ActsGeometry.h>
2 
3 bool writer = true;
4 
5 const int nlayers = 7;
6 
7 namespace
8 {
9  template<class T> inline constexpr T square(const T& x) { return x*x;}
10  template<class T> inline constexpr T get_r(const T&x ,const T&y)
11  { return std::sqrt(square(x)+square(y)); }
12  inline const double get_angle(const Acts::Vector3&x, const Acts::Vector3& y)
13  { return std::acos(x.dot(y)/(x.norm() * y.norm())); }
14 }
15 
16 void AnalyzeResiduals(std::string infile, const bool use_helical = false)
17 {
18  TFile *file = TFile::Open(infile.c_str());
19  if(!file)
20  {
21  std::cout << "wrong file name" << std::endl;
22  return;
23  }
24  auto tree = (TTree*)file->Get("residualtree");
25 
26  float px, py, pz, pt, eta, phi, deltapt, quality, pcax, pcay, pcaz;
27  int nhits, nmaps, nintt, ntpc, nmms, charge;
28 
29  std::vector<float>* cluslx=0;
30  std::vector<float>* cluslz=0;
31  std::vector<float>* cluselx=0;
32  std::vector<float>* cluselz=0;
33  std::vector<float>* clusgx=0;
34  std::vector<float>* clusgy=0;
35  std::vector<float>* clusgz=0;
36  std::vector<int>* cluslayer=0;
37  std::vector<int>* clussize=0;
38  std::vector<uint32_t>* clushitsetkey = 0;
39 
40  std::vector<float>* idealsurfcenterx=0;
41  std::vector<float>* idealsurfcentery=0;
42  std::vector<float>* idealsurfcenterz=0;
43  std::vector<float>* idealsurfnormx=0;
44  std::vector<float>* idealsurfnormy=0;
45  std::vector<float>* idealsurfnormz=0;
46  std::vector<float>* missurfcenterx=0;
47  std::vector<float>* missurfcentery=0;
48  std::vector<float>* missurfcenterz=0;
49  std::vector<float>* missurfnormx=0;
50  std::vector<float>* missurfnormy=0;
51  std::vector<float>* missurfnormz=0;
52  std::vector<float>* clusgxideal=0;
53  std::vector<float>* clusgyideal=0;
54  std::vector<float>* clusgzideal=0;
55  std::vector<float>* missurfalpha=0;
56  std::vector<float>* missurfbeta=0;
57  std::vector<float>* missurfgamma=0;
58  std::vector<float>* idealsurfalpha=0;
59  std::vector<float>* idealsurfbeta=0;
60  std::vector<float>* idealsurfgamma=0;
61 
62  std::vector<float>* statelx=0;
63  std::vector<float>* statelz=0;
64  std::vector<float>* stateelx=0;
65  std::vector<float>* stateelz=0;
66  std::vector<float>* stategx=0;
67  std::vector<float>* stategy=0;
68  std::vector<float>* stategz=0;
69  std::vector<float>* statepx=0;
70  std::vector<float>* statepy=0;
71  std::vector<float>* statepz=0;
72  std::vector<float>* statepl=0;
73 
74  std::vector<float>* statelxglobderivdx=0;
75  std::vector<float> *statelxglobderivdy=0;
76  std::vector<float> *statelxglobderivdz=0;
77  std::vector<float> *statelxglobderivdalpha=0;
78  std::vector<float> *statelxglobderivdbeta=0;
79  std::vector<float> *statelxglobderivdgamma=0;
80 
81  std::vector<float> *statelxlocderivd0=0;
82  std::vector<float> *statelxlocderivz0=0;
83  std::vector<float> *statelxlocderivphi=0;
84  std::vector<float> *statelxlocderivtheta=0;
85  std::vector<float> *statelxlocderivqop=0;
86 
87  std::vector<float> *statelzglobderivdx=0;
88  std::vector<float> *statelzglobderivdy=0;
89  std::vector<float> *statelzglobderivdz=0;
90  std::vector<float> *statelzglobderivdalpha=0;
91  std::vector<float> *statelzglobderivdbeta=0;
92  std::vector<float> *statelzglobderivdgamma=0;
93 
94  std::vector<float> *statelzlocderivd0=0;
95  std::vector<float> *statelzlocderivz0=0;
96  std::vector<float> *statelzlocderivphi=0;
97  std::vector<float> *statelzlocderivtheta=0;
98  std::vector<float> *statelzlocderivqop=0;
99 
100  tree->SetBranchAddress("px",&px);
101  tree->SetBranchAddress("py",&py);
102  tree->SetBranchAddress("pz",&pz);
103  tree->SetBranchAddress("pt",&pt);
104  tree->SetBranchAddress("eta",&eta);
105  tree->SetBranchAddress("phi",&phi);
106  tree->SetBranchAddress("deltapt",&deltapt);
107  tree->SetBranchAddress("quality",&quality);
108  tree->SetBranchAddress("nhits",&nhits);
109  tree->SetBranchAddress("nmaps",&nmaps);
110  tree->SetBranchAddress("nintt",&nintt);
111  tree->SetBranchAddress("ntpc",&ntpc);
112  tree->SetBranchAddress("nmms",&nmms);
113  tree->SetBranchAddress("quality",&quality);
114  tree->SetBranchAddress("pcax", &pcax);
115  tree->SetBranchAddress("pcay",&pcay);
116  tree->SetBranchAddress("pcaz", &pcaz);
117 
118  tree->SetBranchAddress("clushitsetkey",&clushitsetkey);
119  tree->SetBranchAddress("cluslx",&cluslx);
120  tree->SetBranchAddress("cluslz",&cluslz);
121  tree->SetBranchAddress("cluselx",&cluselx);
122  tree->SetBranchAddress("cluselz",&cluselz);
123  tree->SetBranchAddress("clusgx",&clusgx);
124  tree->SetBranchAddress("clusgy",&clusgy);
125  tree->SetBranchAddress("clusgz",&clusgz);
126  tree->SetBranchAddress("cluslayer",&cluslayer);
127  tree->SetBranchAddress("clussize",&clussize);
128  tree->SetBranchAddress("idealsurfcenterx",&idealsurfcenterx);
129  tree->SetBranchAddress("idealsurfcentery",&idealsurfcentery);
130  tree->SetBranchAddress("idealsurfcenterz",&idealsurfcenterz);
131  tree->SetBranchAddress("missurfcenterx",&missurfcenterx);
132  tree->SetBranchAddress("missurfcentery",&missurfcentery);
133  tree->SetBranchAddress("missurfcenterz",&missurfcenterz);
134  tree->SetBranchAddress("idealsurfnormx",&idealsurfnormx);
135  tree->SetBranchAddress("idealsurfnormy",&idealsurfnormy);
136  tree->SetBranchAddress("idealsurfnormz",&idealsurfnormz);
137  tree->SetBranchAddress("missurfnormx",&missurfnormx);
138  tree->SetBranchAddress("missurfnormy",&missurfnormy);
139  tree->SetBranchAddress("missurfnormz",&missurfnormz);
140  tree->SetBranchAddress("clusgxideal",&clusgxideal);
141  tree->SetBranchAddress("clusgyideal", &clusgyideal);
142  tree->SetBranchAddress("clusgzideal",&clusgzideal);
143  tree->SetBranchAddress("missurfalpha",&missurfalpha);
144  tree->SetBranchAddress("missurfbeta",&missurfbeta);
145  tree->SetBranchAddress("missurfgamma",&missurfgamma);
146  tree->SetBranchAddress("idealsurfalpha",&idealsurfalpha);
147  tree->SetBranchAddress("idealsurfbeta",&idealsurfbeta);
148  tree->SetBranchAddress("idealsurfgamma",&idealsurfgamma);
149 
150  tree->SetBranchAddress("statelx",&statelx);
151  tree->SetBranchAddress("statelz",&statelz);
152  tree->SetBranchAddress("stateelx",&stateelx);
153  tree->SetBranchAddress("stateelz",&stateelz);
154  tree->SetBranchAddress("stategx",&stategx);
155  tree->SetBranchAddress("stategy",&stategy);
156  tree->SetBranchAddress("stategz",&stategz);
157  tree->SetBranchAddress("statepx",&statepx);
158  tree->SetBranchAddress("statepy",&statepy);
159  tree->SetBranchAddress("statepz",&statepz);
160  tree->SetBranchAddress("statepl",&statepl);
161 
162  tree->SetBranchAddress("statelxglobderivdx", &statelxglobderivdx);
163  tree->SetBranchAddress("statelxglobderivdy", &statelxglobderivdy);
164  tree->SetBranchAddress("statelxglobderivdz", &statelxglobderivdz);
165  tree->SetBranchAddress("statelxglobderivdalpha", &statelxglobderivdalpha);
166  tree->SetBranchAddress("statelxglobderivdbeta", &statelxglobderivdbeta);
167  tree->SetBranchAddress("statelxglobderivdgamma", &statelxglobderivdgamma);
168 
169  tree->SetBranchAddress("statelxlocderivd0",&statelxlocderivd0);
170  tree->SetBranchAddress("statelxlocderivz0",&statelxlocderivz0);
171  tree->SetBranchAddress("statelxlocderivphi",&statelxlocderivphi);
172  tree->SetBranchAddress("statelxlocderivtheta",&statelxlocderivtheta);
173  tree->SetBranchAddress("statelxlocderivqop",&statelxlocderivqop);
174 
175  tree->SetBranchAddress("statelzglobderivdx", &statelzglobderivdx);
176  tree->SetBranchAddress("statelzglobderivdy", &statelzglobderivdy);
177  tree->SetBranchAddress("statelzglobderivdz", &statelzglobderivdz);
178  tree->SetBranchAddress("statelzglobderivdalpha", &statelzglobderivdalpha);
179  tree->SetBranchAddress("statelzglobderivdbeta", &statelzglobderivdbeta);
180  tree->SetBranchAddress("statelzglobderivdgamma", &statelzglobderivdgamma);
181 
182  tree->SetBranchAddress("statelzlocderivd0",&statelzlocderivd0);
183  tree->SetBranchAddress("statelzlocderivz0",&statelzlocderivz0);
184  tree->SetBranchAddress("statelzlocderivphi",&statelzlocderivphi);
185  tree->SetBranchAddress("statelzlocderivtheta",&statelzlocderivtheta);
186  tree->SetBranchAddress("statelzlocderivqop",&statelzlocderivqop);
187 
188  TFile *outfile = new TFile("residualoutfile.root","recreate");
189 
190  TH2 *h_residuallayerx = new TH2F("residuallayerx",";x residual [mm];layer",1000,-1,1,60,0,60);
191  TH2 *h_residuallayerz = new TH2F("residuallayerz",";z residual [mm];layer",1000,-1,1,60,0,60);
192  TH2 *h_residualcluspulllayerx = new TH2F("residualcluspulllayerx",";x residual / clus error;layer",400,-10,10,60,0,60);
193  TH2 *h_residualcluspulllayerz = new TH2F("residualcluspulllayerz",";z residual / clus error;layer",400,-10,10,60,0,60);
194  TH2 *h_residualstatepulllayerx = new TH2F("residualstatepulllayerx",";x residual / state error;layer",400,-10,10,60,0,60);
195  TH2 *h_residualstatepulllayerz = new TH2F("residualstatepulllayerz",";z residual / state error;layer",400,-10,10,60,0,60);
196 
197  // as a function of layer
198  TH2 *h_residualxclusr[nlayers];
199  TH2 *h_residualxclusphi[nlayers];
200  TH2 *h_residualgxclusphi[nlayers];
201  TH2 *h_residualgyclusphi[nlayers];
202  TH2 *h_residualgzclusphi[nlayers];
203  TH2 *h_residualxclusz[nlayers];
204  TH2 *h_residualzclusr[nlayers];
205  TH2 *h_residualzclusphi[nlayers];
206  TH2 *h_residualzclusz[nlayers];
207  TH2 *h_surfradiusdiffvsphi[nlayers];
208  TH2 *h_surfxdiffvsphi[nlayers];
209  TH2 *h_surfydiffvsphi[nlayers];
210  TH2 *h_surfzdiffvsphi[nlayers];
211  TH2 *h_surfphidiffvsphi[nlayers];
212  TH2 *h_surfadiffvsphi[nlayers];
213  TH2 *h_surfbdiffvsphi[nlayers];
214  TH2 *h_surfgdiffvsphi[nlayers];
215  TH2 *h_surfadiffvsz[nlayers];
216  TH2 *h_surfbdiffvsz[nlayers];
217  TH2 *h_surfgdiffvsz[nlayers];
218  TH3 *h_residualxclusphiz[nlayers];
219  TH3 *h_residualzclusphiz[nlayers];
220 
221  for(int i=0; i<nlayers; i++)
222  {
223  int zlow = -105;
224  int zhigh = 105;
225  int rlow = 30;
226  int rhigh = 85;
227 
228  if(i<3){
229  zlow = -10; zhigh = 10; rlow = 0; rhigh = 4;
230  } else if (i < 7) {
231  zlow = -30; zhigh = 30; rlow = 4; rhigh = 14;
232  }
233 
234  ostringstream name;
235  name.str("");
236  name << "surfphidiffvsphi_"<<i;
237  h_surfphidiffvsphi[i] = new TH2F(name.str().c_str(),";Ideal #phi_{surf} [deg]; Ideal - misaligned phi [deg]",120,-180,180,1080,-2,2);
238 
239  name.str("");
240  name << "surfadiffvsphi_"<<i;
241  h_surfadiffvsphi[i] = new TH2F(name.str().c_str(),";Ideal #phi_{surf} [deg]; Ideal - misaligned angle to x axis [mrad]",120,-180,180,100,-10,10);
242 
243  name.str("");
244  name << "surfbdiffvsphi_"<<i;
245  h_surfbdiffvsphi[i] = new TH2F(name.str().c_str(),";Ideal #phi_{surf} [deg]; Ideal - misaligned angle to y axis [mrad]",120,-180,180,100,-10,10);
246 
247 name.str("");
248  name << "surfgdiffvsphi_"<<i;
249  h_surfgdiffvsphi[i] = new TH2F(name.str().c_str(),";Ideal #phi_{surf} [deg]; Ideal - misaligned angle to z axis [mrad]",120,-180,180,100,-10,10);
250 
251  name.str("");
252  name << "surfadiffvsz_"<<i;
253  h_surfadiffvsz[i] = new TH2F(name.str().c_str(),";Ideal z_{surf} [mm]; Ideal - misaligned angle to x axis [mrad]",120,-180,180,100,-10,10);
254 
255  name.str("");
256  name << "surfbdiffvsz_"<<i;
257  h_surfbdiffvsz[i] = new TH2F(name.str().c_str(),";Ideal z_{surf} [mm]; Ideal - misaligned angle to y axis [mrad]",120,-180,180,100,-10,10);
258 
259  name.str("");
260  name << "surfgdiffvsz_"<<i;
261  h_surfgdiffvsz[i] = new TH2F(name.str().c_str(),";Ideal z_{surf} [deg]; Ideal - misaligned angle to z axis [mrad]",120,-180,180,100,-10,10);
262 
263  name.str("");
264  name << "surfradiusdiffvsphi_"<<i;
265  h_surfradiusdiffvsphi[i] = new TH2F(name.str().c_str(),";Ideal #phi_{surf} [deg]; Ideal radius - misalign radius [mm]",120,-180,180,1000,-1,1);
266 
267  name.str("");
268  name << "surfxdiffvsphi_"<<i;
269  h_surfxdiffvsphi[i] = new TH2F(name.str().c_str(),";Ideal #phi_{surf} [deg]; Misalign x - ideal x [mm]",120,-180,180,1000,-1,1);
270 
271  name.str("");
272  name << "surfydiffvsphi_"<<i;
273  h_surfydiffvsphi[i] = new TH2F(name.str().c_str(),";Ideal #phi_{surf} [deg]; Misalign y - ideal y [mm]", 120,-180,180,1000,-1,1);
274 
275  name.str("");
276  name << "surfzdiffvsphi_"<<i;
277  h_surfzdiffvsphi[i] = new TH2F(name.str().c_str(), ";Ideal #phi_{surf} [deg]; Misalign z - ideal z [mm]",120,-180,180,2500,-10,10);
278 
279  name.str("");
280  name << "residualxclusphiz_"<<i;
281  //h_residualxclusphiz[i] = new TH3F(name.str().c_str(),";#phi [deg]; z [mm]; x residual [mm]",180,-180,180,100,-200,200,1000,-1,1);
282  name.str("");
283  name << "residualzclusphiz_"<<i;
284  //h_residualzclusphiz[i] = new TH3F(name.str().c_str(),";#phi [deg]; z[mm]; z residual [mm]",180,-180,180,100,-200,200,1000,-1,1);
285 
286  name.str("");
287  name<<"residualgxclusphi_"<<i;
288  h_residualgxclusphi[i] = new TH2F(name.str().c_str(),";#phi [deg]; x_{glob} residual [mm]", 180,-180,180,2000,-1,1);
289 
290  name.str("");
291  name << "residualgyclusphi_"<<i;
292  h_residualgyclusphi[i] = new TH2F(name.str().c_str(),";#phi [deg]; y_{glob} residual [mm]",180,-180,180,2000,-1,1);
293 
294  name.str("");
295  name << "residualgzclusphi_"<<i;
296  h_residualgzclusphi[i] = new TH2F(name.str().c_str(),";#phi [deg]; z_{glob} residual [mm]",180,-180,180,2000,-10,10);
297 
298  name.str("");
299  name << "residualxclusr_"<<i;
300  h_residualxclusr[i] = new TH2F(name.str().c_str(),";r [cm]; x residual [mm]",
301  170,rlow,rhigh,2000,-1,1);
302  name.str("");
303  name << "residualxclusphi_"<<i;
304  h_residualxclusphi[i] = new TH2F(name.str().c_str(),";phi [deg]; x residual [mm]",
305  180,-180,180,2000,-1,1);
306  name.str("");
307  name << "residualxclusz_"<<i;
308  h_residualxclusz[i] = new TH2F(name.str().c_str(),";z [cm]; x residual [mm]",
309  210,zlow,zhigh,2000,-1,1);
310  name.str("");
311  name<<"residualzclusr_"<<i;
312  h_residualzclusr[i] = new TH2F(name.str().c_str(),";r [cm]; z residual [mm]",
313  170,rlow,rhigh,2000,-1,1);
314  name.str("");
315  name << "residualzclusphi_"<<i;
316  h_residualzclusphi[i] = new TH2F(name.str().c_str(),";phi [deg]; z residual [mm]",
317  180,-180,180,2000,-1,1);
318  name.str("");
319  name << "residualzclusz_"<<i;
320  h_residualzclusz[i] = new TH2F(name.str().c_str(),";z [cm]; z residual [mm]",
321  210,zlow,zhigh,2000,-1,1);
322  }
323 
324 
325 
326  TH2F *h_residualhitsetkeyx = new TH2F("residualhitsetkeyx",";x residual [mm];hitsetkey",1000,-1,1,1200,0,1200);
327  TH2F *h_residualhitsetkeyz = new TH2F("residualhitsetkeyz",";z residual [mm];hitsetkey",1000,-1,1,1200,0,1200);
328  TH2F *h_residualcluspullhitsetkeyx = new TH2F("residualcluspullhitsetkeyx",";x residual / clus error;hitsetkey",400,-10,10,1200,0,1200);
329  TH2F *h_residualcluspullhitsetkeyz = new TH2F("residualcluspullhitsetkeyz",";z residual / clus error;hitsetkey",400,-10,10,1200,0,1200);
330  TH2F *h_residualstatepullhitsetkeyx = new TH2F("residualstatepullhitsetkeyx",";x residual / state error;hitsetkey",400,-10,10,1200,0,1200);
331  TH2F *h_residualstatepullhitsetkeyz = new TH2F("residualstatepullhitsetkeyz",";z residual / state error;hitsetkey",400,-10,10,1200,0,1200);
332 
333  TH2F *h_xresidualgdx = new TH2F("xresidualgdx",";dx_{resid}/dx_{glob};x_{resid} [mm]",1000,-1.3,1.3,1000,-1,1);
334  TH2F *h_xresidualgdy = new TH2F("xresidualgdy",";dx_{resid}/dy_{glob};x_{resid} [mm]",1000,-1.3,1.3,1000,-1,1);
335  TH2F *h_xresidualgdz = new TH2F("xresidualgdz",";dx_{resid}/dz_{glob};x_{resid} [mm]",1000,-0.001,0.001,1000,-1,1);
336  TH2F *h_xresidualdalpha = new TH2F("xresidualdalpha",";dx_{resid}/d#alpha;x_{resid} [mm]",1000,-600,600,1000,-1,1);
337  TH2F *h_xresidualdbeta = new TH2F("xresidualdbeta",";dx_{resid}/d#beta;x_{resid} [mm]",1000,-600,600,1000,-1,1);
338  TH2F *h_xresidualdgamma = new TH2F("xresidualdgamma",";dx_{resid}/d#gamma;x_{resid} [mm]",1000,-20,20,1000,-1,1);
339  TH2F *h_xresidualdd0 = new TH2F("xresidualdd0",";dx_{resid}/dd_{0}; x_{resid} [mm]",1000,0.9,1.1,1000,-1,1);
340  TH2F *h_xresidualdz0 = new TH2F("xresidualdz0",";dx_{resid}/dz_{0}; x_{resid} [mm]",1000,-0.002,0.002,1000,-1,1);
341  TH2F *h_xresidualdphi = new TH2F("xresidualdphi",";dx_{resid}/d#phi; x_{resid} [mm]", 1000,-100,100,1000,-1,1);
342  TH2F *h_xresidualdtheta = new TH2F("xresidualdtheta",";dx_{resid}/d#theta; x_{resid} [mm]",1000,-10,10,1000,-1,1);
343  TH2F *h_xresidualdqop = new TH2F("xresidualdqop",";dx_{resid}/dqop; x_{resid} [mm]",1000,-20,10,1000,-1,1);
344 
345 
346  TH2F *h_zresidualgdx = new TH2F("zresidualgdx",";dz_{resid}/dx_{glob};z_{resid} [mm]",1000,-1.3,1.3,1000,-1,1);
347  TH2F *h_zresidualgdy = new TH2F("zresidualgdy",";dz_{resid}/dy_{glob};z_{resid} [mm]",1000,-1.3,1.3,1000,-1,1);
348  TH2F *h_zresidualgdz = new TH2F("zresidualgdz",";dz_{resid}/dz_{glob};z_{resid} [mm]",1000,-1.3,1.3,1000,-1,1);
349  TH2F *h_zresidualdalpha = new TH2F("zresidualdalpha",";dz_{resid}/d#alpha;z_{resid} [mm]",1000,-600,600,1000,-1,1);
350  TH2F *h_zresidualdbeta = new TH2F("zresidualdbeta",";dz_{resid}/d#beta;z_{resid} [mm]",1000,-600,600,1000,-1,1);
351  TH2F *h_zresidualdgamma = new TH2F("zresidualdgamma",";dz_{resid}/d#gamma;z_{resid} [mm]",1000,-20,20,1000,-1,1);
352  TH2F *h_zresidualdd0 = new TH2F("zresidualdd0",";dz_{resid}/dd_{0}; z_{resid} [mm]",1000,-1,1,1000,-1,1);
353  TH2F *h_zresidualdz0 = new TH2F("zresidualdz0",";dz_{resid}/dz_{0}; z_{resid} [mm]",1000,0.9,1.1,1000,-1,1);
354  TH2F *h_zresidualdphi = new TH2F("zresidualdphi",";dz_{resid}/d#phi; z_{resid} [mm]", 1000,-100,100,1000,-1,1);
355  TH2F *h_zresidualdtheta = new TH2F("zresidualdtheta",";dz_{resid}/d#theta; z_{resid} [mm]",5000,-600,10,1000,-1,1);
356  TH2F *h_zresidualdqop = new TH2F("zresidualdqop",";dz_{resid}/dqop; z_{resid} [mm]",1000,-2,2,1000,-1,1);
357 
358  TH2F *h_dxdxvsphi = new TH2F("dxdxvsphi",";#phi [deg];dx_{resid}/dx_{glob}",120,-180,180,1000,-1.3,1.3);
359  TH2F *h_dxdyvsphi = new TH2F("dxdyvsphi",";#phi [deg];dx_{resid}/dy_{glob}",120,-180,180,1000,-1.3,1.3);
360  TH2F *h_dxdzvsphi = new TH2F("dxdzvsphi",";#phi [deg];dx_{resid}/dz_{glob}",120,-180,180,1000,-0.002,0.002);
361  TH2F *h_dxdalphavsphi = new TH2F("dxdalphavsphi",";#phi [deg];dx_{resid}/d#alpha",120,-180,180,1000,-600,600);
362  TH2F *h_dxdbetavsphi = new TH2F("dxdbetavsphi",";#phi [deg];dx_{resid}/d#beta",120,-180,180,1000,-600,600);
363  TH2F *h_dxdgammavsphi = new TH2F("dxdgammavsphi",";#phi [deg];dx_{resid}/d#gamma",120,-180,180,1000,-20,20);
364  TH2F *h_dxdd0vsphi = new TH2F("dxdd0vsphi",";#phi [deg];dx_{resid}/dd_{0}",120,-180,180,1000,0.9,1.1);
365  TH2F *h_dxdz0vsphi = new TH2F("dxdz0vsphi",";#phi [deg];dx_{resid}/dz_{0}",120,-180,180,1000,-0.002,0.002);
366  TH2F *h_dxdphivsphi = new TH2F("dxdphivsphi",";#phi [deg];dx_{resid}/d#phi",120,-180,180,1000,-100,100);
367  TH2F *h_dxdthetavsphi = new TH2F("dxdthetavsphi",";#phi [deg];dx_{resid}/d#theta",120,-180,180,1000,-10,10);
368  TH2F *h_dxdqopvsphi = new TH2F("dxdqopvsphi",";#phi [deg];dx_{resid}/dqop",120,-180,180,1000,-1,1);
369 
370  TH2F *h_dzdxvsphi = new TH2F("dzdxvsphi",";#phi [deg];dz_{resid}/dx_{glob}",120,-180,180,1000,-1.3,1.3);
371  TH2F *h_dzdyvsphi = new TH2F("dzdyvsphi",";#phi [deg];dz_{resid}/dy_{glob}",120,-180,180,1000,-1.3,1.3);
372  TH2F *h_dzdzvsphi = new TH2F("dzdzvsphi",";#phi [deg];dz_{resid}/dz_{glob}",120,-180,180,1000,-1.3,1.3);
373  TH2F *h_dzdalphavsphi = new TH2F("dzdalphavsphi",";#phi [deg];dz_{resid}/d#alpha",120,-180,180,1000,-600,600);
374  TH2F *h_dzdbetavsphi = new TH2F("dzdbetavsphi",";#phi [deg];dz_{resid}/d#beta",120,-180,180,1000,-600,600);
375  TH2F *h_dzdgammavsphi = new TH2F("dzdgammavsphi",";#phi [deg];dz_{resid}/d#gamma",120,-180,180,1000,-20,20);
376  TH2F *h_dzdd0vsphi = new TH2F("dzdd0vsphi",";#phi [deg];dz_{resid}/dd_{0}",120,-180,180,1000,-1,1);
377  TH2F *h_dzdz0vsphi = new TH2F("dzdz0vsphi",";#phi [deg];dz_{resid}/dz_{0}",120,-180,180,1000,0.9,1.1);
378  TH2F *h_dzdphivsphi = new TH2F("dzdphivsphi",";#phi [deg];dz_{resid}/d#phi",120,-180,180,1000,-100,100);
379  TH2F *h_dzdthetavsphi = new TH2F("dzdthetavsphi",";#phi [deg];dz_{resid}/d#theta",120,-180,180,1000,-10,10);
380  TH2F *h_dzdqopvsphi = new TH2F("dzdqopvsphi",";#phi [deg];dz_{resid}/dqop",120,-180,180,1000,-1,1);
381 
382  TH2F *h_dxdxvseta = new TH2F("dxdxvseta",";#eta ;dx_{resid}/dx_{glob}",100,-1,1,1000,-1.3,1.3);
383  TH2F *h_dxdyvseta = new TH2F("dxdyvseta",";#eta ;dx_{resid}/dy_{glob}",100,-1,1,1000,-1.3,1.3);
384  TH2F *h_dxdzvseta = new TH2F("dxdzvseta",";#eta ;dx_{resid}/dz_{glob}",100,-1,1,1000,-1.3,1.3);
385  TH2F *h_dxdalphavseta = new TH2F("dxdalphavseta",";#eta ;dx_{resid}/d#alpha",100,-1,1,1000,-600,600);
386  TH2F *h_dxdbetavseta = new TH2F("dxdbetavseta",";#eta ;dx_{resid}/d#beta",100,-1,1,1000,-600,600);
387  TH2F *h_dxdgammavseta = new TH2F("dxdgammavseta",";#eta ;dx_{resid}/d#gamma",100,-1,1,1000,-20,20);
388  TH2F *h_dxdd0vseta = new TH2F("dxdd0vseta",";#eta ;dx_{resid}/dd_{0}",100,-1,1,1000,0.9,1.1);
389  TH2F *h_dxdz0vseta = new TH2F("dxdz0vseta",";#eta ;dx_{resid}/dz_{0}",100,-1,1,1000,-0.002,0.002);
390  TH2F *h_dxdphivseta = new TH2F("dxdphivseta",";#eta ;dx_{resid}/d#phi",100,-1,1,1000,-100,100);
391  TH2F *h_dxdthetavseta = new TH2F("dxdthetavseta",";#eta ;dx_{resid}/d#theta",100,-1,1,1000,-10,10);
392  TH2F *h_dxdqopvseta = new TH2F("dxdqopvseta",";#eta ;dx_{resid}/dqop",100,-1,1,1000,-1,1);
393 
394  TH2F *h_dzdxvseta = new TH2F("dzdxvseta",";#eta ;dz_{resid}/dx_{glob}",100,-1,1,1000,-1.3,1.3);
395  TH2F *h_dzdyvseta = new TH2F("dzdyvseta",";#eta ;dz_{resid}/dy_{glob}",100,-1,1,1000,-1.3,1.3);
396  TH2F *h_dzdzvseta = new TH2F("dzdzvseta",";#eta ;dz_{resid}/dz_{glob}",100,-1,1,1000,-1.3,1.3);
397  TH2F *h_dzdalphavseta = new TH2F("dzdalphavseta",";#eta ;dz_{resid}/d#alpha",100,-1,1,1000,-600,600);
398  TH2F *h_dzdbetavseta = new TH2F("dzdbetavseta",";#eta ;dz_{resid}/d#beta",100,-1,1,1000,-600,600);
399  TH2F *h_dzdgammavseta = new TH2F("dzdgammavseta",";#eta ;dz_{resid}/d#gamma",100,-1,1,1000,-20,20);
400  TH2F *h_dzdd0vseta = new TH2F("dzdd0vseta",";#eta ;dz_{resid}/dd_{0}",100,-1,1,1000,-1,1);
401  TH2F *h_dzdz0vseta = new TH2F("dzdz0vseta",";#eta ;dz_{resid}/dz_{0}",100,-1,1,1000,0.9,1.1);
402  TH2F *h_dzdphivseta = new TH2F("dzdphivseta",";#eta ;dz_{resid}/d#phi",100,-1,1,1000,-100,100);
403  TH2F *h_dzdthetavseta = new TH2F("dzdthetavseta",";#eta ;dz_{resid}/d#theta",100,-1,1,1000,-10,10);
404  TH2F *h_dzdqopvseta = new TH2F("dzdqopvseta",";#eta ;dz_{resid}/dqop",100,-1,1,1000,-1,1);
405  TH2F *h_sensorlayercounts = new TH2F("sensorlayercounts",";sensor number; layer",1200,0,1200,58,0,58);
406 
407  TH2F *h_lowcountsxy = new TH2F("lowcountsxy",";x_{clus} [cm];y_{clus} [cm];",
408  10000,-20,20,10000,-20,20);
409 
410  TH2F *h_lowcountsrz = new TH2F("lowcountsrz",";z_{clus} [cm];r_{clus} [cm];",
411  10000,-40,40,100,0,20);
412 
413  int badent = 0;
414  int sensornum = 0;
415  std::map<uint32_t, int> hitsetkey_sensornum_map, hitsetkey_numcounts_map;
416  std::cout << "iterating " << tree->GetEntries() << std::endl;
417  for(int i=0; i<tree->GetEntries(); i++)
418  {
419  if(i > 25000)
420  break;
421  tree->GetEntry(i);
422  for(int j=0; j<cluslx->size(); j++)
423  {
424  unsigned int layer = cluslayer->at(j);
425 
426  if(layer > 6) continue;
427  float xresidual = cluslx->at(j) - statelx->at(j);
428  float zresidual = cluslz->at(j) - statelz->at(j);
429  //convert to mm
430  xresidual *= 10;
431  zresidual *= 10;
432  float celx = cluselx->at(j);
433  float celz = cluselz->at(j);
434  float cgz = clusgz->at(j);
435  celx *= 10;
436  celz *= 10;
437  cgz *= 10;
438  uint32_t hitsetkey = clushitsetkey->at(j);
439 
440  auto iter = hitsetkey_sensornum_map.find(hitsetkey);
441  if(iter == hitsetkey_sensornum_map.end())
442  {
443  hitsetkey_sensornum_map.insert(std::make_pair(hitsetkey,sensornum));
444  hitsetkey_numcounts_map.insert(std::make_pair(hitsetkey, 1));
445  sensornum++;
446  }
447  else
448  {
449  hitsetkey_numcounts_map.find(hitsetkey)->second++;
450  }
451 
452  int sensornum = hitsetkey_sensornum_map.find(hitsetkey)->second;
453  if(sensornum > 130)
454  {
455  h_lowcountsxy->Fill(clusgx->at(j), clusgy->at(j));
456  h_lowcountsrz->Fill(clusgz->at(j),sqrt(clusgx->at(j)*clusgx->at(j)+clusgy->at(j)*clusgy->at(j)));
457  }
458  h_sensorlayercounts->Fill(sensornum, layer);
459  TVector3 clusvec;
460  clusvec.SetXYZ(clusgx->at(j), clusgy->at(j), clusgz->at(j));
461  float clusphi = clusvec.Phi() * 180 / 3.14159;
462  float cluseta = clusvec.Eta();
463  float statelxglobderivdxv = statelxglobderivdx->at(j);
464 
465  float xresidualcluspull = xresidual / celx;
466  float zresidualcluspull = zresidual / celz;
467  float xresidualstatepull = xresidual / (10*stateelx->at(j));
468  float zresidualstatepull = zresidual / (10*stateelz->at(j));
469 
470  float statelxglobdx = statelxglobderivdx->at(j);
471  float statelxglobdy = statelxglobderivdy->at(j);
472  float statelxglobdz = statelxglobderivdz->at(j);
473  float statelxglobdalpha = statelxglobderivdalpha->at(j);
474  float statelxglobdbeta = statelxglobderivdbeta->at(j);
475  float statelxglobdgamma = statelxglobderivdgamma->at(j);
476  float statelzglobdx = statelxglobderivdx->at(j);
477  float statelzglobdy = statelxglobderivdy->at(j);
478  float statelzglobdz = statelxglobderivdz->at(j);
479  float statelzglobdalpha = statelzglobderivdalpha->at(j);
480  float statelzglobdbeta = statelzglobderivdbeta->at(j);
481  float statelzglobdgamma = statelzglobderivdgamma->at(j);
482 
483  if(use_helical)
484  {
485  statelxglobdalpha *= 10;
486  statelxglobdbeta *= 10;
487  statelxglobdgamma *= 10;
488  statelzglobdalpha *= 10;
489  statelzglobdbeta *= 10;
490  statelzglobdgamma *= 10;
491  }
492 
493  Acts::Vector3 missurfnorm(missurfnormx->at(j),
494  missurfnormy->at(j),
495  missurfnormz->at(j));
496 
497  Acts::Vector3 idealsurfnorm(idealsurfnormx->at(j),
498  idealsurfnormy->at(j),
499  idealsurfnormz->at(j));
500 
501  float alpha = idealsurfalpha->at(j)-missurfalpha->at(j);
502  float beta = idealsurfbeta->at(j)-missurfbeta->at(j);
503  float gamma = idealsurfgamma->at(j)-missurfgamma->at(j);
504 
506  alpha *= 1000;
507  beta *=1000;
508  gamma *=1000;
509  float idealsurfphi = 180/M_PI*atan2(idealsurfcentery->at(j),
510  idealsurfcenterx->at(j));
511  float missurfphi = 180/M_PI*atan2(missurfcentery->at(j),
512  missurfcenterx->at(j));
513  h_surfphidiffvsphi[layer]->Fill(idealsurfphi, missurfphi-idealsurfphi);
514  h_surfadiffvsphi[layer]->Fill(idealsurfphi, alpha);
515  h_surfbdiffvsphi[layer]->Fill(idealsurfphi, beta);
516  h_surfgdiffvsphi[layer]->Fill(idealsurfphi, gamma);
517  h_surfadiffvsz[layer]->Fill(idealsurfcenterz->at(j)*10, alpha);
518  h_surfbdiffvsz[layer]->Fill(idealsurfcenterz->at(j)*10, beta);
519  h_surfgdiffvsz[layer]->Fill(idealsurfcenterz->at(j)*10, gamma);
520 
521  h_surfradiusdiffvsphi[layer]->Fill(idealsurfphi,
522  get_r(missurfcenterx->at(j),
523  missurfcentery->at(j))*10 -
524  get_r(idealsurfcenterx->at(j),
525  idealsurfcentery->at(j))*10);
526  h_surfxdiffvsphi[layer]->Fill(idealsurfphi,
527  missurfcenterx->at(j)*10 -
528  idealsurfcenterx->at(j)*10);
529  h_surfydiffvsphi[layer]->Fill(idealsurfphi,
530  missurfcentery->at(j)*10 -
531  idealsurfcentery->at(j)*10);
532 
533  h_surfzdiffvsphi[layer]->Fill(idealsurfphi,
534  missurfcenterz->at(j)*10 -
535  idealsurfcenterz->at(j)*10);
536  h_residualxclusr[layer]->Fill(clusvec.Perp(), xresidual);
537  h_residualgxclusphi[layer]->Fill(clusphi, (clusgx->at(j) - stategx->at(j))*10);
538  h_residualgyclusphi[layer]->Fill(clusphi, (clusgy->at(j) - stategy->at(j))*10);
539  //std::cout << clusgz->at(j) << ", " << stategz->at(j) << std::endl;
540  h_residualgzclusphi[layer]->Fill(clusphi, (clusgz->at(j) - stategz->at(j))*10);
541 
542  h_residualxclusphi[layer]->Fill(clusphi, xresidual);
543  h_residualxclusz[layer]->Fill(cgz, xresidual);
544  h_residualzclusr[layer]->Fill(clusvec.Perp(), zresidual);
545  h_residualzclusphi[layer]->Fill(clusphi, zresidual);
546  h_residualzclusz[layer]->Fill(cgz, zresidual);
547  //h_residualxclusphiz[layer]->Fill(clusphi,cgz,xresidual);
548  //h_residualzclusphiz[layer]->Fill(clusphi,cgz,zresidual);
549 
550  h_dxdxvsphi->Fill(clusphi, statelxglobdx);
551  h_dxdyvsphi->Fill(clusphi, statelxglobdy);
552  h_dxdzvsphi->Fill(clusphi, statelxglobdz);
553  h_dxdalphavsphi->Fill(clusphi, statelxglobdalpha);
554  h_dxdbetavsphi->Fill(clusphi, statelxglobdbeta);
555  h_dxdgammavsphi->Fill(clusphi, statelxglobdgamma);
556  h_dxdd0vsphi->Fill(clusphi,statelxlocderivd0->at(j));
557  h_dxdz0vsphi->Fill(clusphi,statelxlocderivz0->at(j));
558  h_dxdphivsphi->Fill(clusphi,statelxlocderivphi->at(j));
559  h_dxdthetavsphi->Fill(clusphi,statelxlocderivtheta->at(j));
560  h_dxdqopvsphi->Fill(clusphi,statelxlocderivqop->at(j));
561 
562  h_dzdxvsphi->Fill(clusphi, statelzglobdx);
563  h_dzdyvsphi->Fill(clusphi, statelzglobdy);
564  h_dzdzvsphi->Fill(clusphi, statelzglobdz);
565  h_dzdalphavsphi->Fill(clusphi, statelzglobdalpha);
566  h_dzdbetavsphi->Fill(clusphi, statelzglobdbeta);
567  h_dzdgammavsphi->Fill(clusphi, statelzglobdgamma);
568  h_dzdd0vsphi->Fill(clusphi,statelzlocderivd0->at(j));
569  h_dzdz0vsphi->Fill(clusphi,statelzlocderivz0->at(j));
570  h_dzdphivsphi->Fill(clusphi,statelzlocderivphi->at(j));
571  h_dzdthetavsphi->Fill(clusphi,statelzlocderivtheta->at(j));
572  h_dzdqopvsphi->Fill(clusphi,statelzlocderivqop->at(j));
573 
574  h_dxdxvseta->Fill(cluseta, statelxglobdx);
575  h_dxdyvseta->Fill(cluseta, statelxglobdy);
576  h_dxdzvseta->Fill(cluseta, statelxglobdz);
577  h_dxdalphavseta->Fill(cluseta, statelxglobdalpha);
578  h_dxdbetavseta->Fill(cluseta, statelxglobdbeta);
579  h_dxdgammavseta->Fill(cluseta, statelxglobdgamma);
580  h_dxdd0vseta->Fill(cluseta,statelxlocderivd0->at(j));
581  h_dxdz0vseta->Fill(cluseta,statelxlocderivz0->at(j));
582  h_dxdphivseta->Fill(cluseta,statelxlocderivphi->at(j));
583  h_dxdthetavseta->Fill(cluseta,statelxlocderivtheta->at(j));
584  h_dxdqopvseta->Fill(cluseta,statelxlocderivqop->at(j));
585 
586  h_dzdxvseta->Fill(cluseta, statelzglobdx);
587  h_dzdyvseta->Fill(cluseta, statelzglobdy);
588  h_dzdzvseta->Fill(cluseta, statelzglobdz);
589  h_dzdalphavseta->Fill(cluseta, statelzglobdalpha);
590  h_dzdbetavseta->Fill(cluseta, statelzglobdbeta);
591  h_dzdgammavseta->Fill(cluseta, statelzglobdgamma);
592  h_dzdd0vseta->Fill(cluseta,statelzlocderivd0->at(j));
593  h_dzdz0vseta->Fill(cluseta,statelzlocderivz0->at(j));
594  h_dzdphivseta->Fill(cluseta,statelzlocderivphi->at(j));
595  h_dzdthetavseta->Fill(cluseta,statelzlocderivtheta->at(j));
596  h_dzdqopvseta->Fill(cluseta,statelzlocderivqop->at(j));
597 
598  h_xresidualgdx->Fill(statelxglobdx,xresidual);
599  h_xresidualgdy->Fill(statelxglobdy,xresidual);
600  h_xresidualgdz->Fill(statelxglobdz,xresidual);
601  h_xresidualdalpha->Fill(statelxglobdalpha,xresidual);
602  h_xresidualdbeta->Fill(statelxglobdbeta,xresidual);
603  h_xresidualdgamma->Fill(statelxglobdgamma,xresidual);
604  h_xresidualdd0->Fill(statelxlocderivd0->at(j),xresidual);
605  h_xresidualdz0->Fill(statelxlocderivz0->at(j),xresidual);
606  h_xresidualdphi->Fill(statelxlocderivphi->at(j),xresidual);
607  h_xresidualdtheta->Fill(statelxlocderivtheta->at(j),xresidual);
608  h_xresidualdqop->Fill(statelxlocderivqop->at(j),xresidual);
609 
610  h_zresidualgdx->Fill(statelzglobdx,zresidual);
611  h_zresidualgdy->Fill(statelzglobdy,zresidual);
612  h_zresidualgdz->Fill(statelzglobdz,zresidual);
613  h_zresidualdalpha->Fill(statelzglobdalpha,zresidual);
614  h_zresidualdbeta->Fill(statelzglobdbeta,zresidual);
615  h_zresidualdgamma->Fill(statelzglobdgamma,zresidual);
616  h_zresidualdd0->Fill(statelzlocderivd0->at(j),zresidual);
617  h_zresidualdz0->Fill(statelzlocderivz0->at(j),zresidual);
618  h_zresidualdphi->Fill(statelzlocderivphi->at(j),zresidual);
619  h_zresidualdtheta->Fill(statelzlocderivtheta->at(j),zresidual);
620  h_zresidualdqop->Fill(statelzlocderivqop->at(j),zresidual);
621 
622  h_residuallayerx->Fill(xresidual, layer);
623  h_residuallayerz->Fill(zresidual, layer);
624  h_residualcluspulllayerx->Fill(xresidualcluspull, layer);
625  h_residualcluspulllayerz->Fill(zresidualcluspull, layer);
626  h_residualstatepulllayerx->Fill(xresidualstatepull, layer);
627  h_residualstatepulllayerz->Fill(zresidualstatepull, layer);
628 
629  h_residualhitsetkeyx->Fill(xresidual, sensornum);
630  h_residualhitsetkeyz->Fill(zresidual, sensornum);
631  h_residualcluspullhitsetkeyx->Fill(xresidualcluspull, sensornum);
632  h_residualcluspullhitsetkeyz->Fill(zresidualcluspull, sensornum);
633  h_residualstatepullhitsetkeyx->Fill(xresidualstatepull, sensornum);
634  h_residualstatepullhitsetkeyz->Fill(zresidualstatepull, sensornum);
635 
636  }
637  }
638  std::cout << "badent " << badent<<std::endl;
639  for(const auto& [key, counts] : hitsetkey_numcounts_map)
640  {
641  std::cout << "hitsetkey " << key << " has " << counts << std::endl;
642  }
643 
644  TCanvas *can1 = new TCanvas("can1","can1",300,300,1000,700);
645  can1->Divide(3,2);
646  can1->cd(1);
647  h_xresidualdalpha->Draw();
648  can1->cd(2);
649  h_xresidualdbeta->Draw();
650  can1->cd(3);
651  h_xresidualdgamma->Draw();
652  can1->cd(4);
653  h_xresidualgdx->Draw();
654  can1->cd(5);
655  h_xresidualgdy->Draw();
656  can1->cd(6);
657  h_xresidualgdz->Draw();
658 
659 
660  TCanvas *can2 = new TCanvas("can2","can2",300,300,1000,700);
661  can2->Divide(3,2);
662  can2->cd(1);
663  h_xresidualdd0->Draw();
664  can2->cd(2);
665  h_xresidualdz0->Draw();
666  can2->cd(3);
667  h_xresidualdphi->Draw();
668  can2->cd(4);
669  h_xresidualdtheta->Draw();
670  can2->cd(5);
671  h_xresidualdqop->Draw();
672 
673 
674  TCanvas *can3 = new TCanvas("can3","can3",300,300,1000,700);
675  can3->Divide(3,2);
676  can3->cd(1);
677  h_zresidualdalpha->Draw();
678  can3->cd(2);
679  h_zresidualdbeta->Draw();
680  can3->cd(3);
681  h_zresidualdgamma->Draw();
682  can3->cd(4);
683  h_zresidualgdx->Draw();
684  can3->cd(5);
685  h_zresidualgdy->Draw();
686  can3->cd(6);
687  h_zresidualgdz->Draw();
688 
689 
690  TCanvas *can4 = new TCanvas("can4","can4",300,300,1000,700);
691  can4->Divide(3,2);
692  can4->cd(1);
693  h_zresidualdd0->Draw();
694  can4->cd(2);
695  h_zresidualdz0->Draw();
696  can4->cd(3);
697  h_zresidualdphi->Draw();
698  can4->cd(4);
699  h_zresidualdtheta->Draw();
700  can4->cd(5);
701  h_zresidualdqop->Draw();
702 
703 
704  TCanvas *can5 = new TCanvas("can5","can5",300,300,1000,700);
705  can5->Divide(3,2);
706  can5->cd(1);
707 
708  h_dxdxvsphi->Draw();
709  can5->cd(2);
710  h_dxdyvsphi->Draw();
711  can5->cd(3);
712  h_dxdzvsphi->Draw();
713  can5->cd(4);
714  h_dxdalphavsphi->Draw();
715  can5->cd(5);
716  h_dxdbetavsphi->Draw();
717  can5->cd(6);
718  h_dxdgammavsphi->Draw();
719 
720  TCanvas *can6 = new TCanvas("can6","can6",300,300,1000,700);
721  can6->Divide(3,2);
722  can6->cd(1);
723 
724  h_dxdd0vsphi->Draw();
725  can6->cd(2);
726  h_dxdz0vsphi->Draw();
727  can6->cd(3);
728  h_dxdphivsphi->Draw();
729  can6->cd(4);
730  h_dxdthetavsphi->Draw();
731  can6->cd(5);
732  h_dxdqopvsphi->Draw();
733 
734  TCanvas *can7 = new TCanvas("can7","can7",300,300,1000,700);
735  can7->Divide(3,2);
736  can7->cd(1);
737 
738  h_dzdxvsphi->Draw();
739  can7->cd(2);
740  h_dzdyvsphi->Draw();
741  can7->cd(3);
742  h_dzdzvsphi->Draw();
743  can7->cd(4);
744  h_dzdalphavsphi->Draw();
745  can7->cd(5);
746  h_dzdbetavsphi->Draw();
747  can7->cd(6);
748  h_dzdgammavsphi->Draw();
749 
750  TCanvas *can8 = new TCanvas("can8","can8",300,300,1000,700);
751  can8->Divide(3,2);
752  can8->cd(1);
753 
754  h_dzdd0vsphi->Draw();
755  can8->cd(2);
756  h_dzdz0vsphi->Draw();
757  can8->cd(3);
758  h_dzdphivsphi->Draw();
759  can8->cd(4);
760  h_dzdthetavsphi->Draw();
761  can8->cd(5);
762  h_dzdqopvsphi->Draw();
763 
764 
765  TCanvas *can9 = new TCanvas("can9","can9",300,300,1000,700);
766  can9->Divide(3,2);
767  can9->cd(1);
768 
769  h_dxdxvseta->Draw();
770  can9->cd(2);
771  h_dxdyvseta->Draw();
772  can9->cd(3);
773  h_dxdzvseta->Draw();
774  can9->cd(4);
775  h_dxdalphavseta->Draw();
776  can9->cd(5);
777  h_dxdbetavseta->Draw();
778  can9->cd(6);
779  h_dxdgammavseta->Draw();
780 
781  TCanvas *can10 = new TCanvas("can10","can10",300,300,1000,700);
782  can10->Divide(3,2);
783  can10->cd(1);
784 
785  h_dxdd0vseta->Draw();
786  can10->cd(2);
787  h_dxdz0vseta->Draw();
788  can10->cd(3);
789  h_dxdphivseta->Draw();
790  can10->cd(4);
791  h_dxdthetavseta->Draw();
792  can10->cd(5);
793  h_dxdqopvseta->Draw();
794 
795  TCanvas *can11 = new TCanvas("can11","can11",300,300,1000,700);
796  can11->Divide(3,2);
797  can11->cd(1);
798 
799  h_dzdxvseta->Draw();
800  can11->cd(2);
801  h_dzdyvseta->Draw();
802  can11->cd(3);
803  h_dzdzvseta->Draw();
804  can11->cd(4);
805  h_dzdalphavseta->Draw();
806  can11->cd(5);
807  h_dzdbetavseta->Draw();
808  can11->cd(6);
809  h_dzdgammavseta->Draw();
810 
811  TCanvas *can12 = new TCanvas("can12","can12",300,300,1000,700);
812  can12->Divide(3,2);
813  can12->cd(1);
814 
815  h_dzdd0vseta->Draw();
816  can12->cd(2);
817  h_dzdz0vseta->Draw();
818  can12->cd(3);
819  h_dzdphivseta->Draw();
820  can12->cd(4);
821  h_dzdthetavseta->Draw();
822  can12->cd(5);
823  h_dzdqopvseta->Draw();
824 
825  if(writer)
826  {
827  outfile->cd();
828  h_lowcountsxy->Write();
829  h_lowcountsrz->Write();
830  h_residuallayerx->Write();
831  h_residuallayerz->Write();
832  h_residualcluspulllayerx->Write();
833  h_residualcluspulllayerz->Write();
834  h_residualstatepulllayerx->Write();
835  h_residualstatepulllayerz->Write();
836  h_sensorlayercounts->Write();
837 
838  h_residualhitsetkeyx->Write();
839  h_residualhitsetkeyz->Write();
840  h_residualcluspullhitsetkeyx->Write();
841  h_residualcluspullhitsetkeyz->Write();
842  h_residualstatepullhitsetkeyx->Write();
843  h_residualstatepullhitsetkeyz->Write();
844 
845  h_xresidualgdx->Write();
846  h_xresidualgdy->Write();
847  h_xresidualgdz->Write();
848  h_xresidualdalpha->Write();
849  h_xresidualdbeta->Write();
850  h_xresidualdgamma->Write();
851  h_xresidualdd0->Write();
852  h_xresidualdz0->Write();
853  h_xresidualdphi->Write();
854  h_xresidualdtheta->Write();
855  h_xresidualdqop->Write();
856 
857  h_zresidualgdx->Write();
858  h_zresidualgdy->Write();
859  h_zresidualgdz->Write();
860  h_zresidualdalpha->Write();
861  h_zresidualdbeta->Write();
862  h_zresidualdgamma->Write();
863  h_zresidualdd0->Write();
864  h_zresidualdz0->Write();
865  h_zresidualdphi->Write();
866  h_zresidualdtheta->Write();
867  h_zresidualdqop->Write();
868 
869  h_dxdxvsphi->Write();
870  h_dxdyvsphi->Write();
871  h_dxdzvsphi->Write();
872  h_dxdalphavsphi->Write();
873  h_dxdbetavsphi->Write();
874  h_dxdgammavsphi->Write();
875  h_dxdd0vsphi->Write();
876  h_dxdz0vsphi->Write();
877  h_dxdphivsphi->Write();
878  h_dxdthetavsphi->Write();
879  h_dxdqopvsphi->Write();
880 
881  h_dzdxvsphi->Write();
882  h_dzdyvsphi->Write();
883  h_dzdzvsphi->Write();
884  h_dzdalphavsphi->Write();
885  h_dzdbetavsphi->Write();
886  h_dzdgammavsphi->Write();
887  h_dzdd0vsphi->Write();
888  h_dzdz0vsphi->Write();
889  h_dzdphivsphi->Write();
890  h_dzdthetavsphi->Write();
891  h_dzdqopvsphi->Write();
892 
893  h_dxdxvseta->Write();
894  h_dxdyvseta->Write();
895  h_dxdzvseta->Write();
896  h_dxdalphavseta->Write();
897  h_dxdbetavseta->Write();
898  h_dxdgammavseta->Write();
899  h_dxdd0vseta->Write();
900  h_dxdz0vseta->Write();
901  h_dxdphivseta->Write();
902  h_dxdthetavseta->Write();
903  h_dxdqopvseta->Write();
904 
905  h_dzdxvseta->Write();
906  h_dzdyvseta->Write();
907  h_dzdzvseta->Write();
908  h_dzdalphavseta->Write();
909  h_dzdbetavseta->Write();
910  h_dzdgammavseta->Write();
911  h_dzdd0vseta->Write();
912  h_dzdz0vseta->Write();
913  h_dzdphivseta->Write();
914  h_dzdthetavseta->Write();
915  h_dzdqopvseta->Write();
916 
917 
918  for(int i=0; i<nlayers; i++)
919  {
920  h_residualxclusr[i]->Write();
921  h_residualxclusz[i]->Write();
922  h_residualgxclusphi[i]->Write();
923  h_residualgyclusphi[i]->Write();
924  h_residualgzclusphi[i]->Write();
925  h_residualzclusr[i]->Write();
926  h_residualzclusz[i]->Write();
927  h_residualxclusphi[i]->Write();
928  h_residualzclusphi[i]->Write();
929 
930  //h_residualxclusphiz[i]->Write();
931  //h_residualzclusphiz[i]->Write();
932  h_surfradiusdiffvsphi[i]->Write();
933  h_surfxdiffvsphi[i]->Write();
934  h_surfphidiffvsphi[i]->Write();
935  h_surfydiffvsphi[i]->Write();
936  h_surfzdiffvsphi[i]->Write();
937  h_surfadiffvsphi[i]->Write();
938  h_surfbdiffvsphi[i]->Write();
939  h_surfgdiffvsphi[i]->Write();
940 
941  h_surfadiffvsz[i]->Write();
942  h_surfbdiffvsz[i]->Write();
943  h_surfgdiffvsz[i]->Write();
944  }
945 
946  outfile->Write();
947  outfile->Close();
948 
949 
950  }
951 
952 
953 
954 }