StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
T.C
1 #define T_xxx
2 #include <assert.h>
3 #include "TH2.h"
4 #include "TStyle.h"
5 #include "TCanvas.h"
6 #include "TF1.h"
7 #include "TRVector.h"
8 #include "TRMatrix.h"
9 #include "TRSymMatrix.h"
10 #include "T.h"
11 #include "Riostream.h"
12 ClassImp(TT);
13 struct Geometry_t {
14  Int_t Barrel;
15  Int_t Layer;
16  Int_t NoLadders;
17  Int_t NoWafers;
18 };
19 const Int_t NoLayers = 7;
20 // Barrel, Layer ladder wafer
21 const Geometry_t SvtSsdConfig[] =
22  { {1, 1, 8, 4}, // even
23  {1, 2, 8, 4}, // odd
24  {2, 3, 12, 6}, // event
25  {2, 4, 12, 6}, // odd
26  {3, 5, 16, 7}, // even
27  {3, 6, 16, 7}, // odd
28  {4, 7, 20, 16} // Ssd
29  };
30 const Int_t BL[4] = {8, 12, 16, 20}; // ladders in barrel
31 struct HybridFit_t {
32  HybridFit_t() {noentries=0; AmX = TRVector(6); S = TRSymMatrix(6);}
33  Int_t noentries;
34  TRVector AmX;
35  TRSymMatrix S;
36 };
37 //________________________________________________________________________________
38 Double_t STcheb(Int_t N, Double_t *par, Double_t x) {// N polynome degree, dimension is par[N+1]
39  if (N < 0 || N > 12) return 0;
40  Double_t T0 = 1;
41  Double_t T1 = 2*x - 1;
42  Double_t T2;
43  Double_t Sum = par[0]*T0;
44  if (N >= 1) {
45  T1 = 2*x - 1;
46  Sum += par[1]*T1;
47  for (int n = 2; n <= N; n++) {
48  T2 = 2*(2*x - 1)*T1 - T0;
49  Sum += par[n]*T2;
50  T0 = T1;
51  T1 = T2;
52  }
53  }
54  return Sum;
55 }
56 #if 0
57 //________________________________________________________________________________
58 Double_t DriftCorHack(Int_t barrel, Int_t ladder, Int_t wafer, Double_t u) {
59  struct data_t {
60  Int_t barrel, layer, ladder, wafer, hybrid, Npar;
61  Double_t param[12];
62  Double_t dparam[12];
63  Char_t *Comment;
64  };
65  static Int_t N = 0;
66  static data_t *pointers[3][16][7][2];
67  if (N == 0) {
68  N = sizeof(Data)/sizeof(data_t);
69  memset (pointers,0, 3*16*7*2*sizeof(data_t *));
70  }
71  Int_t h = 1;
72  if (u > 0) h = 2;
73  assert(barrel >= 1 && barrel <= 3);
74  assert(ladder >= 1 && ladder <= 16);
75  assert(wafer >= 1 && wafer <= 7);
76  data_t *p = pointers[barrel-1][ladder-1][wafer-1][h-1];
77  if (! p) {
78  for (Int_t i = 0; i < N; i++) {
79  if (Data[i].barrel == barrel &&
80  Data[i].ladder == ladder &&
81  Data[i].wafer == wafer &&
82  Data[i].hybrid == h) {
83  p = Data + i;
84  pointers[barrel-1][ladder-1][wafer-1][h-1] = p;
85  break;
86  }
87  }
88  }
89 
90  return p ? STcheb(p->Npar, p->param, TMath::Abs(u/3.)) : 0;
91 }
92 #endif
93 //________________________________________________________________________________
94 void TT::Loop(Int_t Nevents) {
95 // In a ROOT session, you can do:
96 // Root > .L T.C
97 // Root > T t
98 // Root > t.GetEntry(12); // Fill t data members with entry number 12
99 // Root > t.Show(); // Show values of entry 12
100 // Root > t.Show(16); // Read and show values of entry 16
101 // Root > t.Loop(); // Loop on all entries
102 //
103 
104 // This is the loop skeleton where:
105 // jentry is the global entry number in the chain
106 // ientry is the entry number in the current Tree
107 // Note that the argument to GetEntry must be:
108 // jentry for TChain::GetEntry
109 // ientry for TTree::GetEntry and TBranch::GetEntry
110 //
111 // To read only selected branches, Insert statements like:
112 // METHOD1:
113 // fChain->SetBranchStatus("*",0); // disable all branches
114 // fChain->SetBranchStatus("branchname",1); // activate branchname
115 // METHOD2: replace line
116 // fChain->GetEntry(jentry); //read all branches
117 //by b_branchname->GetEntry(ientry); //read only this branch
118  /* Local
119  --------
120  (1 0 0)
121  RotateX(alpha) = (0 1 -alpha)
122  (0 alpha 1)
123 
124  (1 0 beta)
125  RotateY(beta) = (0 1 0)
126  (-beta 0 1)
127 
128  (1 -gamma 0)
129  RotateZ(gamma) = (gamma 1 0)
130  (0 0 1)
131 
132  (1 -gamma beta)
133  Rx*Ry*Rz = (gamma 1 -alpha)
134  (-beta alpha 1)
135 
136 
137  T is transformation from "real" local (l) coordinate system to "known" as local (l') (local -> Master)
138  (u') ( 1 gamma -beta)(u) (du)
139  l'=(v') = (-gamma 1 alpha)(v) + (dv)
140  (w') ( beta -alpha 1)(w) (dw)
141 
142  du dv dw alpha beta gamma
143  (u - uP) = (-1, 0, tuP, tuP*vP, -tuP*uP, vP) =
144  (v - vP) ( 0, -1, tvP, tvP*vP, -tvP*uP, -uP)
145 
146  (u - uP) = -du +tuP*dw +tuP*vP*alpha -tuP*uP*beta +vP*gamma;
147  (v - vP) = -dv +tvP*dw +tvP*vP*alpha -tvP*uP*beta -uP*gamma;
148 or
149  (u - uP) = -du +tuP*(dw +vP*alpha -uP*beta) +vP*gamma;
150  (v - vP) = -dv +tvP*(dw +vP*alpha -uP*beta) -uP*gamma;
151  Assume uniform distribution over tuP, tvP, vP, uP then:
152  <u - uP> => -du;
153  <v - vP> => -dv
154  1. <u - uP> versus tuP => dw
155  2 <v - vP> versus tvP => dw
156  3. <u - uP> versus vP => gamma
157  4 <v - vP> versus -uP => gamma
158  5. <(u - uP)/tuP> versus vP => alpha
159  6 <(v - vP)/tvP> versus vP => alpha
160  7. <(u - uP)/tuP> versus -uP => beta;
161  8. <(v - vP)/tvP> versus -uP => beta;
162  ________________________________________________________________________________
163  Global
164  ------
165  (dx) ( 1 -gamma beta )(xG) (xP)
166  r' = dr + R*r = (dy) + ( gamma 1 -alpha)(yG) => (yP)
167  (dz) ( -beta alpha 1 )(zG) (zP)
168 
169  dX = xP - xG = dx -gamma*yG + beta*zG
170  dY = yP - yG = dy + gamma*xG -alpha*zG
171  dZ = zP - zG = dz -beta*xG +alpha*yG
172 
173  ( 1 gamma -beta )(xP-dx) (xG)
174  R^-1*(r-dr) = (-gamma 1 alpha)(yP-dy) => (yG)
175  ( beta -alpha 1 )(zP-dz) (zG)
176 
177 
178  dX = xG - xP = -dx +gamma*yP -beta*zP
179  dY = yG - yP = -dy -gamma*xP +alpha*zP
180  dZ = zG - zP = -dz +beta*xP -alpha*yP
181 
182  ________________________________________________________________________________
183  <dX> => -dx
184  <dY> => -dy
185  <dZ> => -dz
186  0. dX versus xP
187  1. dX versus yP => gamma
188  2. dX versus -zP => beta
189 
190  3. dY versus -xP => gamma
191  4. dY versus yP
192  5. dY versus zP => alpha
193 
194  6. dZ versus xP => beta
195  7. dZ versus -yP => alpha
196  8. dZ versus zP
197  ________________________________________________________________________________
198  Account that prediction should be on detector plane
199  (xP,yP,zP) == (x,y,z)
200 
201  dX = xG - x;
202  dY = yG - y;
203  dZ = zG - z;
204  (dxP,dyP,dPz) = track prediction direction in GCS
205  (xP,yP,zP) = track prediction position in GCS
206  (wx.wy.wz) = detector plane direction in GCS
207  (vx,vy,vz) = (wx,wy,wz)/dLz, dLz track direction in LCS
208 
209  dX = dx*(-1+dxP*vx) + dy*( dxP*vy) + dz*( dxP*vz) + alpha*( dxP*(-vy*z+vz*y)) + beta*(-z+dxP*(vx*z-vz*x)) + gamma*( y+dxP*(-vx*y+vy*x))
210  dY = dx*( dyP*vx) + dy*(-1+dyP*vy) + dz*( dyP*vz) + alpha*( z+dyP*(-vy*z+vz*y)) + beta*( dyP*(vx*z-vz*x)) + gamma*(-x+dyP*(-vx*y+vy*x))
211  dZ = dx*( dzP*vx) + dy*( dzP*vy) + dz*(-1+dzP*vz) + alpha*(-y+dzP*(-vy*z+vz*y)) + beta*( x+dzP*(vx*z-vz*x)) + gamma*( dzP*(-vx*y+vy*x))
212 
213 
214 ------------------
215  dRT = I + DT;
216  q = (dx,dy,dz,alpha,beta,gamma)
217  A = -I + j * vT; v = wG/(wG*dG)
218  Delta = (dX,dY,dZ) = DRT*(x_q - dt) = DRT*(p + h*j) = A ( dt - DT *p) = A * s;
219  r0 = DR*t + dt); w0 = DR*R*wL = DR*wG
220  w0T *( p + h*j - DR*t - dt) = 0;
221 
222  ( 1 gamma -beta )
223  DRT = (-gamma 1 alpha);
224  ( beta -alpha 1 )
225 
226  ( 0 gamma -beta )
227  DT = DRT - I = (-gamma 0 alpha);
228  ( beta -alpha 0 )
229 
230  ( 0 gamma -beta ) (x) ( gamma*y - beta *z)
231  DT*p = (-gamma 0 alpha)*(y) = (-gamma*x + alpha*z)
232  ( beta -alpha 0 ) (z) ( beta *x -alpha*y )
233 
234 
235 
236  ( dx - gamma*y + beta *z)
237  s = dt - DT * p = ( dy + gamma*x - alpha*z)
238  ( dz - beta *x + alpha*y)
239  q = (dx,dy,dZ,alpha,beta,gamma)
240  ( 1 0 0 0 z -y)
241  ds / d q = B = ( 0 1 0 -z 0 x)
242  ( 0 0 1 y -x 0)
243 
244 
245  (-1 0 0) (jx) ((-1 + jx*vx) jx*vy jx*vz )
246  A = -I + j*vT = ( 0 -1 0) + (jy) * (vx vy vz) = ( jy*vx (-1 + jy*vy) jy*vz )
247  ( 0 0 -1) (jz) ( jz*vx jz*vy (-1 + jz*vz))
248 
249  (-1+jx*vx jx*vy jx*vz) ( 1 0 0 0 z -y)
250  A * B = ( jy*vx -1+jy*vy jy*vz) * ( 0 1 0 -z 0 x) =
251  ( jz*vx jz*vy -1+jz*vz) ( 0 0 1 y -x 0)
252 
253  (-1+jx*vx jx*vy jx*vz jx*(-vy*z+vz*y) -z+jx*(vx*z-vz*x) y+jx*(-vx*y+vy*x))
254  = ( jy*vx -1+jy*vy jy*vz z+jy*(-vy*z+vz*y) jy*(vx*z-vz*x) -x+jy*(-vx*y+vy*x))
255  ( jz*vx jz*vy -1+jz*vz -y+jz*(-vy*z+vz*y) x+jz*(vx*z-vy*x) jz*(-vx*y+vy*x))
256 
257  */
258  if (fChain == 0) return;
259  struct PlotName_t {
260  Char_t *Name;
261  Char_t *Title;
262  Double_t xmax[2]; // svt ssd
263  };
264  // svt ssd
265  // static Double_t Dv[2] = {3.000, 2.10};
266  // Double_t Radii[7] = { 6.37, 7.38, 10.38, 11.27, 14.19, 15.13, 23.80};
267  static Double_t Du[2] = {3.000, 3.65};
268  static Double_t Sv[2] = {6.305, 4.35};
269  const PlotName_t plotNameD[10] = {// plots for drift
270  {"dutuP","<u - uP> versus tuP => dw for Drift", { 0.5, 0.5}}, // 0
271  {"dvtvP","<v - vP> versus tvP => dw for Drift", { 2.5, 2.5}}, // 1
272  {"duvP", "<u - uP> versus v => gamma for Drift", { -2, -2}}, // 2 z
273  {"dvuP", "<v - vP> versus u => -gamma for Drift", { -1, -1}}, // 3
274  {"duOvertuPvP","<(u - uP)/tuP> versus v => alpha for Drift", { -2, -2}}, // 4 z
275  {"dvOvertvPvP","<(v - vP)/tvP> versus v => alpha for Drift", { -2, -2}}, // 5 z
276  {"duOvertuPuP","<(u - uP)/tuP> versus u => -beta for Drift", { -1, -1}}, // 6
277  {"dvOvertvPuP","<(v - vP)/tvP> versus u => -beta for Drift", { -1, -1}}, // 7
278  {"duuH" , "<u - uP> versus uHat for Drift", {1.2, -1}}, // 8
279  {"duvH" , "<u - uP> versus vHat for Drift", {1.2, -2}} // 9 z
280  };
281  const PlotName_t plotName[37] = {
282  {"dutuP","<u - uP> versus tuP => dw", { 0.5, 0.5}}, // 0
283  {"dvtvP","<v - vP> versus tvP => dw", { 2.5, 2.5}}, // 1
284  {"duvP", "<u - uP> versus vP => gamma", { -2, -2}}, // 2 z
285  {"dvuP", "<v - vP> versus -uP => gamma", { -1, -1}}, // 3
286  {"duOvertuPvP","<(u - uP)/tuP> versus vP => alpha", { -2, -2}}, // 4 z
287  {"dvOvertvPvP","<(v - vP)/tvP> versus vP => alpha", { -2, -2}}, // 5 z
288  {"duOvertuPuP","<(u - uP)/tuP> versus -uP => beta", { -1, -1}}, // 6
289  {"dvOvertvPuP","<(v - vP)/tvP> versus -uP => beta", { -1, -1}}, // 7
290  {"duuP", "<u - uP> versus -uP", { -1, -1}}, // 8
291  {"dvvP", "<v - vP> versus vP", { -2, -2}}, // 9 z
292  {"dXvsX","dX versus x" , { 16, 24}}, // 10
293  {"dXvsY","dX versus y => gamma", { 16, 24}}, // 11
294  {"dXvsZ","dX versus -z => beta", { 24.,36}}, // 12
295  {"dYvsX","dY versus -x => gamma", { 16, 24}}, // 13
296  {"dYvsY","dY versus y" , { 16, 24}}, // 14
297  {"dYvsZ","dY versus z => alpha", { 24.,36}}, // 15
298  {"dZvsX","dZ versus x => beta", { 16, 24}}, // 16
299  {"dZvsY","dZ versus -y => alpha", { 16, 24}}, // 17
300  {"dZvsZ","dZ versus z", { 24.,36}}, // 18
301 
302  {"dX4dx","dX vs -1+jx*vx => dx", {2.2,2.2}}, // 19
303  {"dX4dy","dX vs jx*vy => dy", { 1, 1}}, // 20
304  {"dX4dz","dX vs jx*vz => dz", {.01,.01}}, // 21
305  {"dX4da","dX vs jx*(-vy*z+vz*y)=> alpha", {20,20}}, // 22
306  {"dX4db","dX vs -z+jx*(vx*z-vz*x)=> beta ", {40,40}}, // 23
307  {"dX4dg","dX vs y+jx*(-vx*y+vy*x)=> alpha", {25,25}}, // 24
308 
309  {"dY4dx","dY vs jy*vx => dx", { 1, 1}}, // 25
310  {"dY4dy","dY vs -1+jy*vy => dy", {2.2,2.2}}, // 26
311  {"dY4dz","dY vs jy*vz => dz", {.01,.01}}, // 27
312  {"dY4da","dY vs z+jy*(-vy*z+vz*y)=> alpha", {35,35}}, // 28
313  {"dY4db","dY vs jy*(vx*z-vz*x)=> beta ", {20,20}}, // 29
314  {"dY4dg","dY vs -x+jy*(-vx*y+vy*x)=> gamma", {25,25}}, // 30
315 
316  {"dZ4dx","dZ vs jz*vx => dx", {2.2,2.2}}, // 31
317  {"dZ4dy","dZ vs jz*vy => dy", {2.2,2.2}}, // 32
318  {"dZ4dz","dZ vs -1+jz*vz => dz", {2.2,2.2}}, // 33
319  {"dZ4da","dZ vs -y+jz*(-vy*z+vz*y)=> alpha", {80,80}}, // 34
320  {"dZ4db","dZ vs x+jz*( vx*z-vy*x)=> beta ", {80,80}}, // 35
321  {"dZ4dg","dZ vs jz*(-vx*y+vy*x)=> gamma", {10,10}}, // 36
322 
323  };
324 
325  const Int_t ssdSector[20] = {// 100*sector + ladder
326  101, 102,
327  203, 204, 205, 206, 207, 208, 209,
328  310, 311, 312,
329  413, 414, 415, 416, 417, 418, 419,
330  120
331  };
332  TFile *fOut = new TFile(fOutFileName,"recreate");
333  TString Name, Title;
334  TH1D *LSF = new TH1D("LSF","Matrix and right part for Least Squred Fit",6*28,0,6*28);
335  TH1D *LSFB[4];
336  for (Int_t barrel = 1; barrel <= 4; barrel++)
337  LSFB[barrel-1] = new TH1D(Form("LSFB%i",barrel),
338  Form("Matrix and right part for Least Squred Fit for barrel %i",barrel),
339  BL[barrel-1]*28,0,BL[barrel-1]*28);
340  // T B l W
341  TH2F *LocPlots[10][4][20][17];
342  memset(LocPlots,0,10*4*20*17*sizeof(TH2F *));
343  // TH2F *LocPlots[9][4][20][17];
344  // memset(LocPlots,0,9*4*20*17*sizeof(TH2F *));
345  for (Int_t L = 0; L < NoLayers; L++) {// over Layers
346  Int_t barrel = SvtSsdConfig[L].Barrel;
347  Int_t layer = SvtSsdConfig[L].Layer;
348  Int_t NoLadders = SvtSsdConfig[L].NoLadders;
349  Int_t NoWafers = SvtSsdConfig[L].NoWafers;
350  // if (! AllWafers) NoWafers = 2; // use wafer index for Positive / negatives
351  for (Int_t ladder = 1; ladder <= NoLadders; ladder++) {
352  if (barrel <= 3 && (ladder-1)%2 != layer%2) continue;
353  for (Int_t wafer = 0; wafer <= NoWafers; wafer++) {// wafer == 0 for whole ladder
354  Int_t Id = ladder + 100*(wafer + 10*layer);
355  for (Int_t t = 0; t < 10; t++) {
356  if (NoWafers > 2 && wafer != 0) {
357  Name = Form("%s%04i",plotNameD[t].Name,Id);
358  Title = Form("%s for barrel %i, layer %i ladder %i, ",plotNameD[t].Title,barrel,layer,ladder);
359  } else {
360  Name = Form("%s%04i",plotNameD[t].Name,Id);
361  Title = Form("%s for barrel %i, layer %i ladder %i, ",plotNameD[t].Title,barrel,layer,ladder);
362  }
363  if (AllWafers) {
364  if (wafer == 0) Title += "all wafers";
365  else Title += Form("wafer %i",wafer);
366  } else {
367  Title += "all wafers";
368  if (wafer == 1) Title += " Positive";
369  if (wafer == 2) Title += " Negative";
370  }
371  Int_t n = 100;
372  Int_t k = 0;
373  if (barrel > 3) k = 1;
374  Double_t xmax = plotNameD[t].xmax[k];
375 #if 0
376  cout << plotNameD[t].Name << "/" << plotNameD[t].Title
377  << "\txmax " << plotNameD[t].xmax[0] << "\t" << plotNameD[t].xmax[1]
378  << "\txmax = " << xmax << endl;
379 #endif
380  if (xmax > 0 && t >= 8 && ! (NoWafers > 2 && wafer != 0)) xmax = -1;
381  if (xmax < 0) {
382  Int_t m = - (Int_t) xmax;
383  if (m == 1) xmax = Du[k];
384  else xmax = Sv[k]/2.;
385  }
386  if ((wafer == 0 || ! AllWafers) && (t == 2 || t == 4 || t == 5 || t == 9)) {
387  switch (barrel) {
388  case 1: xmax = 12; break;
389  case 2: xmax = 18; break;
390  case 3: xmax = 21; break;
391  case 0:
392  case 4: xmax = 35; break;
393  default: xmax = 40; break;
394  }
395  xmax += 2;
396  n = (Int_t) (4*xmax);
397  }
398  Double_t ymax = rCut/2;
399  Int_t ny = 500;
400  Double_t dy = ymax/ny;
401  // if (t >= 4 && t <= 7) ymax *= 10;
402  LocPlots[t][barrel-1][ladder-1][wafer] = new TH2F(Name,Title,n,-xmax,xmax,ny+1,-ymax-dy,ymax+dy);
403  }
404  }
405  }
406  }
407  // T S
408  TH2F *GloPlots[27][6];
409  memset(GloPlots,0,27*6*sizeof(TH2F *));
410  for (Int_t s = 0; s < 6; s++) {
411  for (Int_t i = 0; i < 27; i++) {
412  Int_t t = i+10;
413  Name = Form("%s%i",plotName[t].Name,s);
414  if (s < 2)
415  Title = Form("%s for SVT Clam shell %i",plotName[t].Title,s);
416  else
417  Title = Form("%s for SSD Sector %i",plotName[t].Title,s-1);
418  Int_t m = 0;
419  if (s > 1) m = 1;
420  Double_t xmax = plotName[t].xmax[m];
421  Int_t n = (Int_t) (4.*xmax);
422  if (n < 100) n = 100;
423  Double_t ymax = rCut;
424  GloPlots[i][s] = new TH2F(Name,Title, n,-xmax,xmax,500,-ymax,ymax);
425  }
426  }
427  // T B L
428  TH2F *GloBLPlots[27][4][20];
429  memset(GloBLPlots,0,27*4*20*sizeof(TH2F *));
430  if (LaddersInGlobal) {
431  for (Int_t barrel = 0; barrel < 4; barrel++) {
432  for (Int_t ladder = 0; ladder < BL[barrel]; ladder++) {
433  for (Int_t i = 0; i < 27; i++) {
434  Int_t t = i+10;
435  Name = Form("%sB%iL%02i",plotName[t].Name,barrel+1,ladder+1);
436  Title = Form("%s for barrel %i ladder %02i",plotName[t].Title,barrel+1,ladder+1);
437  Int_t m = 0;
438  if (barrel > 3) m = 1;
439  Double_t xmax = plotName[t].xmax[m];
440  Int_t n = (Int_t) (4.*xmax);
441  if (n < 100) n = 100;
442  Double_t ymax = 2.50;
443  GloBLPlots[i][barrel][ladder] = new TH2F(Name,Title, n,-xmax,xmax,500,-ymax,ymax);
444  }
445  }
446  }
447  }
448  Long64_t nentries = fChain->GetEntriesFast();
449  if (Nevents > 0 && nentries > Nevents) nentries = Nevents;
450  Long64_t nbytes = 0, nb = 0;
451  Int_t TreeNo = -1;
452  TString currentFile("");
453  for (Long64_t jentry=0; jentry<nentries;jentry++) {
454  Long64_t ientry = LoadTree(jentry);
455  if (ientry < 0) break;
456  nb = fChain->GetEntry(jentry); nbytes += nb;
457  if (! jentry%1000 || TreeNo != fChain->GetTreeNumber()) {
458  if (jentry > 0) fOut->Flush();
459  cout << "Read event \t" << jentry
460  << " so far. switch to file " << fChain->GetCurrentFile()->GetName() << endl;
461  TreeNo = fChain->GetTreeNumber();
462  }
463  if (VertexZCut > 0 && TMath::Abs(fVertex[2]) > VertexZCut) continue;
464  UInt_t Ntrack = fNPTracks;
465  Int_t run = fEvtHdr_fRun;
466  for (UInt_t trk = 0; trk < Ntrack; trk++) {
467  Int_t Npoints = fTracks_fNpoint[trk];
468  if (minNoFitPoints > 0 && Npoints%100 < minNoFitPoints) continue;
469  if (UseSsd && Npoints < 1000) continue;
470  if (UseSvt && Npoints < 100) continue;
471  if (TpcLengthCut > 0 && fTracks_fLength[trk] < TpcLengthCut) continue;
472  if (dEdxCut > 0 && (fTracks_fdEdx[trk] <= 1e-7 || fTracks_fdEdx[trk] > dEdxCut)) continue;
473 
474  if (TMath::Abs(fTracks_fTanL[trk]) > 3) continue;
475  if (EastWest) {
476  Double_t zTPC = fVertex[2] + 60*fTracks_fTanL[trk];
477  Double_t zCut = 0;
478  if (EastWest > 2) zCut = 60.*TMath::SinH(0.5);
479  if (EastWest%2 == 1 && ! (zTPC < -zCut || fTracks_fTanL[trk] < 0)) continue;
480  if (EastWest%2 == 0 && ! (zTPC > zCut || fTracks_fTanL[trk] > 0)) continue;
481  }
482  Int_t Nsp = fTracks_fNsp[trk];
483  // if (Nsp <= 0 || Nsp >= 10) continue;
484  for (Int_t hit = 0; hit < Nsp; hit++) {
485  Int_t k = fTracks_fIdHitT[trk][hit]-1;
486  // for (UInt_t k = 0; k < fNhit; k++) {
487  Int_t barrel = fHits_barrel[k];
488  Int_t layer = fHits_layer[k];
489  Int_t ladder = fHits_ladder[k];
490  Int_t wafer = fHits_wafer[k];
491  Int_t hybrid = fHits_hybrid[k];
492  Double32_t anode = fHits_anode[k];
493  Int_t sector = -1;
494  if (fHits_isTrack[k]) continue; // prediction only
495  if (layer < 7 && fHits_hitFlag[k] > 3) continue;
496  //Run V if (layer < 7 && IsNotValidHybrid(barrel,ladder,wafer,hybrid,run,anode)) continue;
497  //Run VI
498  if (layer < 7 && IsNotValidHybrid(barrel,ladder,wafer,hybrid,run,anode)) continue;
499  if (layer < 7) {
500  sector = 0;
501  if (ladder > SvtSsdConfig[layer-1].NoLadders/2) sector = 1;
502  } else {sector = ssdSector[ladder-1]/100 + 1; barrel = 4;}
503  if (sector < 0 || sector > 5) {
504  cout << "Sector " << sector;
505  cout << " for barrel " << barrel
506  << ", ladder " << ladder
507  << ",wafer " << wafer << "\t";
508  cout << "is not been defined" << endl;
509  continue;
510  }
511  if (! GloPlots[0][sector]) {
512  cout << "GloPlots[0][" << sector << "]";
513  cout << " for barrel " << barrel
514  << ", ladder " << ladder
515  << ",wafer " << wafer << "\t";
516  cout << "is not been defined" << endl;
517  continue;
518  }
519  if (LaddersInGlobal && ! GloBLPlots[0][barrel-1][ladder-1]) {
520  cout << "GloBLPlots[0][" << sector << "]";
521  cout << " for barrel " << barrel
522  << ", ladder " << ladder
523  << ",wafer " << wafer << "\t";
524  cout << "is not been defined" << endl;
525  continue;
526  }
527  if (AllWafers) {
528  if (! LocPlots[0][barrel-1][ladder-1][wafer]) {
529  cout << "locPlots";
530  cout << " for barrel " << barrel
531  << ", ladder " << ladder
532  << ",wafer " << wafer << "\t";
533  cout << "is not been defined" << endl;
534  continue;
535  }
536  }
537  if (DipCut > 0 && TMath::Abs(fHits_pT[k]) < DipCut*fHits_pMom[k]) continue;
538  // Int_t NoWafers = SvtSsdConfig[layer-1].NoWafers;
539  Double32_t xPG = fHits_xPG[k];
540  Double32_t zPG = fHits_zPG[k];
541  Double32_t yPG = fHits_yPG[k];
542  Double32_t uP = fHits_uP[k];
543  Double32_t vP = fHits_vP[k];
544  Double32_t tuP = fHits_tuP[k];
545  Double32_t tvP = fHits_tvP[k];
546  Double32_t xPL = fHits_xPL[k];
547  Double32_t zPL = fHits_zPL[k];
548  Double_t dxP = fHits_cxPG[k];
549  Double_t dyP = fHits_cyPG[k];
550  Double_t dzP = fHits_czPG[k];
551 #ifdef __USE_GLOBAL__
552  if (fGlobal) {
553  xPG = fHits_xPGlG[k];
554  zPG = fHits_zPGlG[k];
555  yPG = fHits_yPGlG[k];
556  uP = fHits_uPGl[k];
557  vP = fHits_vPGl[k];
558  tuP = fHits_tuPGl[k];
559  tvP = fHits_tvPGl[k];
560  xPL = fHits_xPGlL[k];
561  zPL = fHits_zPGlL[k];
562  dxP = fHits_cxPGlG[k];
563  dyP = fHits_cyPGlG[k];
564  dzP = fHits_czPGlG[k];
565  }
566 #endif
567  Double32_t xG = fHits_xG[k];
568  Double32_t yG = fHits_yG[k];
569  Double32_t zG = fHits_zG[k];
570  Double32_t dX = xG - xPG;
571  Double32_t dY = yG - yPG;
572  Double32_t dZ = zG - zPG;
573  Double32_t u = fHits_u[k];
574  Double32_t v = fHits_v[k];
575  Double32_t uHat = fHits_uHat[k];
576  if (hybrid == 1) uHat -= 0.1;
577  if (hybrid == 2) uHat += 0.1;
578  // Double32_t vHat = fHits_vHat[k];
579  Double32_t vHat = 1. - anode/240.;
580  if (hybrid == 1) vHat = - vHat;
581  if (hybrid == 1) vHat -= 0.1;
582  if (hybrid == 2) vHat += 0.1;
583  Double32_t zL = fHits_zL[k];
584  if (barrel <= 3) {zPL -= 23.5250; zL -= 23.5250;}
585  Double32_t xL = fHits_xL[k];
586  Double_t DxL = xL - xPL;
587 #if 0
588  Double32_t yL = fHits_yL[k];
589  Double32_t yPL = fHits_yPL[k];
590  Double_t DyL = yL - yPL;
591 #endif
592  Double_t DzL = zL - zPL;
593  Double32_t du = u - uP;
594 #if 0
595  if (barrel <= 3) du -= DriftCorHack(barrel, ladder, wafer, u);
596 #endif
597  Double32_t dv = v - vP;
598  if (TMath::Abs(fHits_pT[k]) < 0.2) continue;
599  // if (TMath::Abs(du) > 0.5 || TMath::Abs(dv) > 0.5) continue;
600  if (TMath::Abs(du) > rCut || TMath::Abs(dv) > rCut) continue;
601  Int_t m = 0;
602  if (barrel > 3) m = 1;
603  if (TMath::Abs(xPG) > plotName[10].xmax[m] ||
604  TMath::Abs(yPG) > plotName[11].xmax[m] ||
605  TMath::Abs(zPG) > plotName[12].xmax[m]) continue;
606  if (TMath::Abs(uP) > Du[m] || TMath::Abs(vP) > Sv[m]/2) continue;
607  Double_t uA = TMath::Abs(u);
608  Double_t vA = TMath::Abs(v);
609  if (layer < 7) {
610  if (uMax > 0 && uA > uMax) continue;
611  if (uMin > 0 && uA < uMin) continue;
612  if (vMax > 0 && vA > vMax) continue;
613  if (vMin > 0 && vA < vMin) continue;
614  }
615  Double_t x = xPG;
616  Double_t y = yPG;
617  Double_t z = zPG;
618  Double_t jL2 = TMath::Sqrt(1. + tuP*tuP + tvP*tvP);
619  Double_t wG[3] = {fHits_wGu[k],fHits_wGv[k], fHits_wGw[k]};
620  Double_t vx = fHits_wGu[k]*jL2;
621  Double_t vy = fHits_wGv[k]*jL2;
622  Double_t vz = fHits_wGw[k]*jL2;
623  Double_t vars[27][2] = {
624  {dX , x}, // 0
625  {dX , y},
626  {dX , -z},
627  {dY , -x},
628  {dY , y},
629  {dY , z},
630  {dZ , x},
631  {dZ , -y},
632  {dZ , z}, // 8
633  {dX, (-1+dxP*vx) }, // 9
634  {dX, ( dxP*vy) }, // 10
635  {dX, ( dxP*vz) },
636  {dX, ( dxP*(-vy*z+vz*y))},
637  {dX, (-z+dxP*( vx*z-vz*x))},
638  {dX, ( y+dxP*(-vx*y+vy*x))},
639  {dY, ( dyP*vx) }, // 15
640  {dY, (-1+dyP*vy) },
641  {dY, ( dyP*vz) },
642  {dY, ( z+dyP*(-vy*z+vz*y))},
643  {dY, ( dyP*( vx*z-vz*x))},
644  {dY, (-x+dyP*(-vx*y+vy*x))},
645  {dZ, ( dzP*vx) }, // 21
646  {dZ, ( dzP*vy) },
647  {dZ, (-1+dzP*vz) },
648  {dZ, (-y+dzP*(-vy*z+vz*y))},
649  {dZ, ( x+dzP*( vx*z-vz*x))},
650  {dZ, ( dzP*(-vx*y+vy*x))}
651  };
652  for (Int_t l = 0; l < 9; l++) {
653  GloPlots[l][sector]->Fill(vars[l][1],vars[l][0]);
654  if (LaddersInGlobal)
655  GloBLPlots[l][barrel-1][ladder-1]->Fill(vars[l][1],vars[l][0]);
656  }
657  Double_t vxyz[3] = {vx, vy, vz};
658  TRMatrix vR(3,1,vxyz);
659  Double_t dxyzP[3] = {dxP, dyP, dzP};
660  TRMatrix dR(3,1,dxyzP);
661  static TRMatrix UR(TRArray::kUnit,3);
662  TRMatrix AR(dR,TRArray::kAxBT,vR);// cout << "AR\t" << AR;
663  AR -= UR;// cout << "AR\t" << AR;
664  TRMatrix BR(3,6,
665  1. , 0., 0., 0., z,-y,
666  0., 1., 0., -z, 0., x,
667  0., 0., 1., y, -x, 0.);// cout << "BR\t" << BR;
668  TRMatrix ABR(AR,TRArray::kAxB,BR);// cout << "ABR\t" << ABR;
669  Double_t dxyz[3] = {dX, dY, dZ};
670  TRVector mX;
671  TRMatrix A(0,6);
672  for (UInt_t l = 0; l < 3; l++) {
673  if (l == 0 && TMath::Abs(wG[0]) >= 0.999 ||
674  l == 1 && TMath::Abs(wG[1]) >= 0.999) continue;
675  mX.AddRow(&dxyz[l]);
676  A.AddRow(ABR.GetRow(l));
677  for (UInt_t k = 0; k < 6; k++) {
678  UInt_t lk = k + 6*l + 9;
679  GloPlots[lk][sector]->Fill(ABR(l,k),dxyz[l]);
680  if (LaddersInGlobal)
681  GloBLPlots[lk][barrel-1][ladder-1]->Fill(ABR(l,k),dxyz[l]);
682 #if 0
683  if (TMath::Abs(vars[l][0]) < 1.e-3 || TMath::Abs(vars[l][1]) < 1.e-3) {
684  cout << "l\t" << l << "\tdX/dY/dZ = " << dX << "\t" << dY << "\t" << dZ << endl;
685  cout << "x/y/z(PG) = " << x << "\t" << y << "\t" << z << endl;
686  cout << "x/y/z( G) = " << xG << "\t" << yG << "\t" << zG << endl;
687  cout << "dirXYZ = " << dxP << "\t" << dyP << "\t" << dzP << endl;
688  cout << "v xyz = " << vx << "\t" << vy << "\t" << vz << endl;
689  }
690 #endif
691  }
692  }
693  TRVector AmX(A,TRArray::kATxB,mX);// cout << "AmX\t" << AmX << endl;
694  TRSymMatrix SX(A,TRArray::kATxA);// cout << "SX\t" << SX << endl;
695  Double_t *array = LSF->GetArray();
696  Double_t *amX = AmX.GetArray();
697  Double_t *sX = SX.GetArray();
698  Int_t im = 1 + 28*sector;
699  Int_t is = im + 6;
700  TCL::vadd(amX,array+im,array+im,6);
701  TCL::vadd(sX,array+is,array+is,21);
702 
703  TRVector duv(2,du,dv);
704  TRMatrix P(2,6,
705  -1., 0., tuP, tuP*vP, -tuP*uP, vP,
706  0., -1., tvP, tvP*vP, -tvP*uP, -uP);
707  TRVector pm(P,TRArray::kATxB,duv);
708  TRSymMatrix PX(P,TRArray::kATxA);
709  array = LSFB[barrel-1]->GetArray();
710  amX = pm.GetArray();
711  sX = PX.GetArray();
712  im = 1 + 28*(ladder-1);
713  is = im + 6;
714  TCL::vadd(amX,array+im,array+im,6);
715  TCL::vadd(sX,array+is,array+is,21);
716 
717 
718  Double32_t duOvertuP = du/tuP;
719  Double32_t dvOvertvP = dv/tvP;
720 #if 1
721  if (AllWafers) {
722  if (wafer == 0) {
723  LocPlots[0][barrel-1][ladder-1][wafer]->Fill(tuP,du);
724  LocPlots[1][barrel-1][ladder-1][wafer]->Fill(tvP,dv);
725  LocPlots[2][barrel-1][ladder-1][wafer]->Fill( vP,du);
726  LocPlots[3][barrel-1][ladder-1][wafer]->Fill(-uP,dv);
727  LocPlots[4][barrel-1][ladder-1][wafer]->Fill( vP,duOvertuP);
728  LocPlots[5][barrel-1][ladder-1][wafer]->Fill( vP,dvOvertvP);
729  LocPlots[6][barrel-1][ladder-1][wafer]->Fill(-uP,duOvertuP);
730  LocPlots[7][barrel-1][ladder-1][wafer]->Fill(-uP,dvOvertvP);
731  LocPlots[8][barrel-1][ladder-1][wafer]->Fill(-uP,du);
732  LocPlots[9][barrel-1][ladder-1][wafer]->Fill( vP,dv);
733  } else {
734  LocPlots[0][barrel-1][ladder-1][wafer]->Fill(tuP,du);
735  LocPlots[1][barrel-1][ladder-1][wafer]->Fill(tvP,dv);
736  LocPlots[2][barrel-1][ladder-1][wafer]->Fill( v ,du);
737  LocPlots[3][barrel-1][ladder-1][wafer]->Fill( u ,dv);
738  LocPlots[4][barrel-1][ladder-1][wafer]->Fill( v ,duOvertuP);
739  LocPlots[5][barrel-1][ladder-1][wafer]->Fill( v ,dvOvertvP);
740  LocPlots[6][barrel-1][ladder-1][wafer]->Fill( u ,duOvertuP);
741  LocPlots[7][barrel-1][ladder-1][wafer]->Fill( u ,dvOvertvP);
742  LocPlots[8][barrel-1][ladder-1][wafer]->Fill( uHat ,du);
743  LocPlots[9][barrel-1][ladder-1][wafer]->Fill( vHat ,du);
744  }
745  } else {
746  wafer = 1;
747  if (fHits_pT[k] < 0) wafer = 2;
748  LocPlots[0][barrel-1][ladder-1][wafer]->Fill(tuP,du);
749  LocPlots[1][barrel-1][ladder-1][wafer]->Fill(tvP,dv);
750  LocPlots[2][barrel-1][ladder-1][wafer]->Fill( zPL,du);
751  LocPlots[3][barrel-1][ladder-1][wafer]->Fill(-uP,dv);
752  LocPlots[4][barrel-1][ladder-1][wafer]->Fill( zPL,duOvertuP);
753  LocPlots[5][barrel-1][ladder-1][wafer]->Fill( zPL,dvOvertvP);
754  LocPlots[6][barrel-1][ladder-1][wafer]->Fill(-uP,duOvertuP);
755  LocPlots[7][barrel-1][ladder-1][wafer]->Fill(-uP,dvOvertvP);
756  LocPlots[8][barrel-1][ladder-1][wafer]->Fill(-uP,du);
757  LocPlots[9][barrel-1][ladder-1][wafer]->Fill( zPL,dv);
758  }
759  // vP = vP + Sv[m]*(wafer - NoWafers/2 - 0.5);
760 #endif
761  wafer = 0;
762  Double32_t DxLOvertuP = DxL/tuP;
763  Double32_t DzLOvertvP = DzL/tvP;
764  LocPlots[0][barrel-1][ladder-1][wafer]->Fill(tuP,DxL); // => dw
765  LocPlots[1][barrel-1][ladder-1][wafer]->Fill(tvP,DzL); // => dw
766  LocPlots[2][barrel-1][ladder-1][wafer]->Fill( zPL,DxL); // => gamma
767  LocPlots[3][barrel-1][ladder-1][wafer]->Fill(-xPL,DzL); // => gamma (-)
768  LocPlots[4][barrel-1][ladder-1][wafer]->Fill( zPL,DxLOvertuP); //=> alpha
769  LocPlots[5][barrel-1][ladder-1][wafer]->Fill( zPL,DzLOvertvP); //=> alpha
770  LocPlots[6][barrel-1][ladder-1][wafer]->Fill(-xPL,DxLOvertuP); //=> beta (-)
771  LocPlots[7][barrel-1][ladder-1][wafer]->Fill(-xPL,DzLOvertvP); //=> beta (-)
772  LocPlots[8][barrel-1][ladder-1][wafer]->Fill(-xPL,DxL); //
773  LocPlots[9][barrel-1][ladder-1][wafer]->Fill( zPL,DzL); //
774  }
775  }
776  }
777 #if 0
778  if (LSF) {
779  for (Int_t s = 0; s < 6; s++) {
780  Double_t *array = LSF->GetArray();
781  Int_t im = 1 + 28*s;
782  Int_t is = im + 6;
783  cout << "sector " << s << "================================" << endl;
784  TRVector AmX(6,array+im); cout << "AmX " << AmX << endl;
785  TRSymMatrix S(6,array+is); cout << "S " << S << endl;
786  TRSymMatrix SInv(S,TRArray::kInverted); cout << "SInv " << SInv << endl;
787  TRVector X(SInv,TRArray::kSxA,AmX); cout << "X " << X << endl;
788  }
789  }
790 #endif
791  fOut->Write();
792 }
793 #ifndef __CINT__
794 //________________________________________________________________________________
795 void FillNt(HybridFit_t *HFit[4][20][16][2]) {
796  if (! HFit) return;
797  for (Int_t barrel = 1; barrel <= 4; barrel++) {
798  for (Int_t ladder = 1; ladder <= 20; ladder++) {
799  for (Int_t wafer = 1; wafer <= 16; wafer++) {
800  for (Int_t hybrid = 1; hybrid <= 2; hybrid++) {
801  HybridFit_t *fit = HFit[barrel-1][ladder-1][wafer-1][hybrid-1];
802  if (! fit) continue;
803  if (fit->noentries < 100) continue;
804  cout << "B/L/W/H\t" << barrel << "/" << ladder << "/" << wafer << "/" << hybrid << endl;
805  cout << "AmX " << fit->AmX << endl;
806  cout << "S " << fit->S << endl;
807  TRSymMatrix SInv(fit->S,TRArray::kInverted); cout << "SInv " << SInv << endl;
808  TRVector X(SInv,TRArray::kSxA,fit->AmX); cout << "X " << X << endl;
809  for (Int_t i = 0; i < 6; i++) {
810  cout << X(i) << " +/- " << TMath::Sqrt(SInv(i,i)) << endl;
811  }
812  }
813  }
814  }
815  }
816 }
817 //________________________________________________________________________________
818 void TT::MakeNt() {
819  if (fChain == 0) return;
820  // Book
821  Long64_t nentries = fChain->GetEntriesFast();
822  Long64_t nbytes = 0, nb = 0;
823  Int_t TreeNo = -1;
824  TString currentFile("");
825  // B l W H
826  HybridFit_t *HFit[4][20][16][2];
827  memset(HFit, 0, 4*20*16*2*sizeof(HybridFit_t *));
828  for (Long64_t jentry=0; jentry<nentries;jentry++) {
829  Long64_t ientry = LoadTree(jentry);
830  if (ientry < 0) break;
831  nb = fChain->GetEntry(jentry); nbytes += nb;
832  if (! jentry%1000 || TreeNo != fChain->GetTreeNumber()) {
833  cout << "Read event \t" << jentry
834  << " so far. switch to file " << fChain->GetCurrentFile()->GetName() << endl;
835  if (TreeNo > -1) FillNt(HFit);
836  TreeNo = fChain->GetTreeNumber();
837  }
838  if (VertexZCut > 0 && TMath::Abs(fVertex[2]) > VertexZCut) continue;
839  UInt_t Ntrack = fNPTracks;
840  for (UInt_t trk = 0; trk < Ntrack; trk++) {
841  Int_t Npoints = fTracks_fNpoint[trk];
842  if (minNoFitPoints > 0 && Npoints%100 < minNoFitPoints) continue;
843  if (UseSsd && Npoints < 1000) continue;
844  if (UseSvt && Npoints < 100) continue;
845  Int_t Nsp = fTracks_fNsp[trk];
846  // if (Nsp <= 0 || Nsp >= 10) continue;
847  for (Int_t hit = 0; hit < Nsp; hit++) {
848  Int_t k = fTracks_fIdHitT[trk][hit];
849  // for (UInt_t k = 0; k < fNhit; k++) {
850  Int_t barrel = fHits_barrel[k];
851  Int_t layer = fHits_layer[k];
852  if (layer < 7 && fHits_hitFlag[k] > 3) continue;
853  Int_t ladder = fHits_ladder[k];
854  Int_t wafer = fHits_wafer[k];
855  Int_t hybrid = fHits_hybrid[k];
856  Double32_t u = fHits_u[k];
857  Double32_t v = fHits_v[k];
858  Double32_t uP = fHits_uP[k];
859  Double32_t vP = fHits_vP[k];
860  Double32_t du = u - uP;
861  Double32_t dv = v - vP;
862  if (TMath::Abs(fHits_pT[k]) < 0.2) continue;
863  // if (TMath::Abs(du) > 0.5 || TMath::Abs(dv) > 0.5) continue;
864  if (TMath::Abs(du) > rCut || TMath::Abs(dv) > rCut) continue;
865  Double_t uA = TMath::Abs(u);
866  Double_t vA = TMath::Abs(v);
867  if (layer < 7) {
868  if (uMax > 0 && uA > uMax) continue;
869  if (uMin > 0 && uA < uMin) continue;
870  if (vMax > 0 && vA > vMax) continue;
871  if (vMin > 0 && vA < vMin) continue;
872  }
873  HybridFit_t *fit = HFit[barrel-1][ladder-1][wafer-1][hybrid-1];
874  if (! fit) {
875  HFit[barrel-1][ladder-1][wafer-1][hybrid-1] = fit = new HybridFit_t();
876  fit->noentries = 0;
877  }
878  fit->noentries++;
879  Double_t a[12] = {
880  1, u, v, 0, 0, 0,
881  0, 0, 0, 1, u, v
882  };
883  TRMatrix A(2,6,a); // cout << "A\t" << A << endl;
884  TRVector duv(2,du,dv);// cout << "duv\t" << duv << endl;
885  // cout << "fit->AmX\t" << fit->AmX << endl;
886  fit->AmX += TRVector(A,TRArray::kATxB,duv);
887  // cout << "fit->AmX\t" << fit->AmX << endl;
888  // cout << "fit->S\t" << fit->S << endl;
889  fit->S += TRSymMatrix(A,TRArray::kATxA);
890  // cout << "fit->S\t" << fit->S << endl;
891  }
892  }
893  }
894  FillNt(HFit);
895 }
896 #endif
897 //________________________________________________________________________________
898 void TT::Loop4BadAnodes(Int_t Nevents) {
899  // B l W
900  const Int_t B = 3;
901  const Int_t L = 16;
902  const Int_t W = 7;
903  const Int_t H = 2;
904  TFile *fOut = new TFile(fOutFileName,"recreate");
905  TH1F *hists[B][L][W][H];
906  memset(hists, 0, B*L*W*H*sizeof(TH1F *));
907  for (Int_t barrel = 1; barrel <= B; barrel++) {
908  Int_t nl = SvtSsdConfig[2*barrel-1].NoLadders;
909  Int_t nw = SvtSsdConfig[2*barrel-1].NoWafers;
910  for (Int_t ladder = 1; ladder <= nl; ladder++)
911  for (Int_t wafer = 1; wafer <= nw; wafer++)
912  for (Int_t hybrid = 1; hybrid <= H; hybrid++) {
913  hists[barrel-1][ladder-1][wafer-1][hybrid-1] =
914  new TH1F(Form("B%iL%02iW%iH%i",barrel,ladder,wafer,hybrid),
915  Form("Anode for Barrel %i Ladder %i, Wafer %i Hybrid %i",
916  barrel,ladder,wafer,hybrid),
917  240,0,240);
918  }
919  }
920  Long64_t nentries = fChain->GetEntriesFast();
921  if (Nevents > 0 && nentries > Nevents) nentries = Nevents;
922  Long64_t nbytes = 0, nb = 0;
923  Int_t TreeNo = -1;
924  TString currentFile("");
925  for (Long64_t jentry=0; jentry<nentries;jentry++) {
926  Long64_t ientry = LoadTree(jentry);
927  if (ientry < 0) break;
928  nb = fChain->GetEntry(jentry); nbytes += nb;
929  if (! jentry%1000 || TreeNo != fChain->GetTreeNumber()) {
930  if (jentry > 0) fOut->Flush();
931  cout << "Read event \t" << jentry
932  << " so far. switch to file " << fChain->GetCurrentFile()->GetName() << endl;
933  TreeNo = fChain->GetTreeNumber();
934  }
935  for (UInt_t k = 0; k < fNhit; k++) {
936  Int_t barrel = fHits_barrel[k];
937  if (barrel < 1 || barrel > 3) continue;
938  Int_t ladder = fHits_ladder[k];
939  Int_t wafer = fHits_wafer[k];
940  Int_t hybrid = fHits_hybrid[k];
941  TH1F *hist = hists[barrel-1][ladder-1][wafer-1][hybrid-1];
942  if (! hist) continue;
943  Double32_t anode = fHits_anode[k];
944  hist->Fill(anode);
945  }
946  }
947  fOut->Write();
948 }
949 #if 0
950 //________________________________________________________________________________
951 void T::LoopTB(Int_t Nevents) {
952  struct PlotName_t {
953  Char_t *Name;
954  Char_t *Title;
955  Int_t nx;
956  Int_t ny;
957  Double_t xmin;
958  Double_t xmax;
959  Double_t ymin;
960  Double_t ymax;
961  Double_t zmin;
962  Double_t zmax;
963  };
964  const PlotName_t plotNameTB = // plots for time bins and anodes
965  { "timeB","time for 80 anodes", 128, 3, 0,128, 0,3, 0,0 };
966 
967  TFile *fOut = new TFile(fOutFileName,"recreate");
968 
969  TString Name, Title;
970  const Int_t NB = 3;
971  const Int_t NL = 16;
972  const Int_t NW = 7;
973  const Int_t NH = 2;
974  const Int_t NA = 3;
975  // B L W H A
976  TH1F *LocPlots[NB][NL][NW][NH][NA];
977  TH1F * LocAll = new TH1F("All","all", plotNameTB.nx, plotNameTB.xmin, plotNameTB.xmax);
978  TH1F * uAll = new TH1F("Uall","ua", 200, -5., 5.);
979  TH1F * uCut = new TH1F("Ucut","uc", 200, -3., 3.);
980  TH1F * vCut = new TH1F("Vcut","vc", 200, -3., 3.);
981  memset(LocPlots,0,NB*NL*NW*NH*sizeof(TH2F *));
982  for (Int_t L = 0; L < 6; L++) {// over Layers
983  Int_t barrel = SvtSsdConfig[L].Barrel;
984  Int_t layer = SvtSsdConfig[L].Layer;
985  Int_t NoLadders = SvtSsdConfig[L].NoLadders;
986  Int_t NoWafers = SvtSsdConfig[L].NoWafers;
987  for (Int_t ladder = 1; ladder <= NoLadders; ladder++) {
988  if (barrel <= 3 && (ladder-1)%2 != layer%2) continue;
989  for (Int_t wafer = 1; wafer <= NoWafers; wafer++) {// wafer == 0 for whole ladder
990  for (Int_t hybrid = 1; hybrid <= 2; hybrid++) {
991  for (Int_t anode =1; anode <=3; anode++) {
992  Name = plotNameTB.Name;
993  Name += Form("L%02iB%iW%02iH%iA%i", ladder, barrel, wafer, hybrid, anode);
994  Title = Form("%s for layer %i B %i L %i W %i H %i G %i",
995  plotNameTB.Title, layer, barrel, ladder, wafer, hybrid, anode);
996  LocPlots[barrel-1][ladder-1][wafer-1][hybrid-1][anode-1] =
997  new TH1F(Name, Title, plotNameTB.nx, plotNameTB.xmin, plotNameTB.xmax );
998  }
999  }
1000  }
1001  }
1002  }
1003  Long64_t nentries = fChain->GetEntriesFast();
1004  if (Nevents > 0 && nentries > Nevents) nentries = Nevents;
1005  Long64_t nbytes = 0, nb = 0;
1006  Int_t TreeNo = -1;
1007  TString currentFile("");
1008 
1009  for (Long64_t jentry=0; jentry<nentries;jentry++) {
1010  Long64_t ientry = LoadTree(jentry);
1011  if (ientry < 0) break;
1012  nb = fChain->GetEntry(jentry); nbytes += nb;
1013  if (! jentry%1000 || TreeNo != fChain->GetTreeNumber()) {
1014  cout << "Read event \t" << jentry
1015  << " so far, switch to file " << fChain->GetCurrentFile()->GetName() << endl;
1016  TreeNo = fChain->GetTreeNumber();
1017  }
1018 // if (VertexZCut > 0 && TMath::Abs(fVertex[2]) > VertexZCut) continue;
1019  UInt_t Ntrack = fNPTracks;
1020  for (UInt_t trk = 0; trk < Ntrack; trk++) {
1021  Int_t Npoints = fTracks_fNpoint[trk];
1022  if (minNoFitPoints > 0 && Npoints%100 < minNoFitPoints) continue;
1023  if (UseSsd && Npoints < 1000) continue;
1024  if (UseSvt && Npoints < 100) continue;
1025  Int_t Nsp = fTracks_fNsp[trk];
1026  for (Int_t hit = 0; hit < Nsp; hit++) {
1027  Int_t k = fTracks_fIdHitT[trk][hit] - 1;
1028  assert(k>=0);
1029  Int_t barrel = fHits_barrel[k];
1030  Int_t layer = fHits_layer[k];
1031  Int_t ladder = fHits_ladder[k];
1032  Int_t wafer = fHits_wafer[k];
1033  Int_t hybrid = fHits_hybrid[k];
1034 
1035  if (layer <= 6 ){
1036  Double32_t u = fHits_u[k];
1037  Double32_t v = fHits_v[k];
1038  Double32_t uP = fHits_uP[k];
1039  Double32_t vP = fHits_vP[k];
1040  Double32_t du = u - uP;
1041  Double32_t dv = v - vP;
1042  if (TMath::Abs(fHits_pT[k]) < 0.2) continue;
1043  uCut->Fill(du);
1044  if (TMath::Abs(du) > rCut ) continue;
1045  vCut->Fill(dv);
1046  if (TMath::Abs(dv) > rCut) continue;
1047  Double32_t anode = fHits_anode[k];
1048  Double32_t timeb = fHits_timeb[k];
1049  Int_t group = ((Int_t)anode)/80;
1050  if (group >= 3) group = 2;
1051  if (group < 0) group = 0;
1052  LocPlots[barrel-1][ladder-1][wafer-1][hybrid-1][group]->Fill( timeb );
1053  LocAll->Fill( timeb );
1054  uAll->Fill( u );
1055  }
1056  }
1057  }
1058  }
1059  fOut->Write();
1060 }
1061 #endif
1062 
1063 //________________________________________________________________________________
1064 Int_t TT::IsNotValidHybrid(Int_t barrel, Int_t ladder, Int_t wafer, Int_t hybrid, Int_t run, Double_t anode) {
1065  Int_t iok = 0;
1066  if (anode <= 10 || anode >= 230) {iok = 1; goto ENDL;}
1067 #if 0
1068 #if !defined(__CINT__)
1069  struct Hybrids_t {
1070  Char_t hybrid[10];
1071  Char_t run[5];
1072  Int_t npeaks; // -88 noise only, -99 dead
1073  Double_t mu0, sigma0;
1074  Double_t mu1, sigma1;
1075  Double_t mu2, sigma2;
1076  Double_t mu3, sigma3;
1077  Double_t mu4, sigma4;
1078  Double_t mu5, sigma5;
1079  Double_t mu6, sigma6;
1080  Double_t mu7, sigma7;
1081  Double_t mu8, sigma8;
1082  Double_t mu9, sigma9;
1083  };
1084  static const Hybrids_t Hybrids[] = {
1085  };
1086  static const UInt_t NoHybrids = sizeof(Hybrids)/sizeof(Hybrids_t);
1087  static TString oldHybrid("");
1088  static TString oldRun("");
1089  static Bool_t InitDone = kFALSE;
1090  static const Int_t NB = 3;
1091  static const Int_t NL =12;
1092  static const Int_t NW = 7;
1093  static const Int_t NH = 2;
1094  static const Int_t ND = 3;
1095  static Hybrids_t *ptr[NB][NL][NW][NH][ND];
1096  if (! InitDone) {
1097  InitDone = kTRUE;
1098  memset(ptr, 0, NB*NL*NW*NH*ND*sizeof(Hybrids_t *));
1099  Int_t b, l, w, h, d, day;
1100  for (UInt_t k = 0; k < NoHybrids; k++) {
1101  sscanf(&Hybrids[k].hybrid[0],"B%iL%02iW%iH%i",&b,&l,&w,&h);
1102  if (b > 0 && b <= NB && l > 0 && l <= NL && w > 0 && w <= NW && h > 0 && h <= NH) {
1103  sscanf(&Hybrids[k].run[0],"0%2iD",&day);
1104  d = 0;
1105  if (day == 49) d = 1;
1106  if (day == 69) d = 2;
1107  ptr[b-1][l-1][w-1][h-1][d-1] = (Hybrids_t *) &Hybrids[k];
1108  }
1109  }
1110  }
1111  Int_t day = (run/1000)%1000;
1112  Int_t d = 0;
1113  Hybrids_t *p = 0;
1114  Int_t k;
1115  Double_t *peaks = 0;
1116  if (day == 49) d = 1;
1117  if (day == 69) d = 2;
1118  if (barrel > 0 && barrel <= NB && ladder > 0 && ladder <= NL && wafer > 0 && wafer <= NW && hybrid > 0 && hybrid <= NH && d < ND) {
1119  p = ptr[barrel-1][ladder-1][wafer-1][hybrid-1][d-1];
1120  if (! p) {iok = 4; goto ENDL;}
1121  if (p->npeaks < 0) {iok = 2; goto ENDL;}
1122  peaks = (Double_t *) &p->mu0;
1123  for (k = 0; k < p->npeaks; k++) {
1124  if (TMath::Abs(anode-peaks[2*k]) < 3*peaks[2*k+1]) {iok = 3; goto ENDL;}
1125  }
1126  }
1127 #endif
1128 #endif
1129  ENDL:
1130  return iok;
1131 }
Geometry_t
Definition: T.C:13
TT
Definition: T.h:18
HybridFit_t
Definition: T.C:31
data_t
Definition: MakeShellSectorOnGlobal.C:1
Generated by   doxygen 1.8.5