StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
CheckMirrors.C
1 //#define DEBUG
2 #if !defined(__CINT__) || defined(__MAKECINT__)
3 #include "Riostream.h"
4 #include <stdio.h>
5 #include "TSystem.h"
6 #include "TMath.h"
7 #include "TH1.h"
8 #include "TH2.h"
9 #include "TH3.h"
10 #include "TProfile.h"
11 #include "TStyle.h"
12 #include "TF1.h"
13 #include "TTree.h"
14 #include "TChain.h"
15 #include "TFile.h"
16 #include "TNtuple.h"
17 #include "TCanvas.h"
18 #include "TMinuit.h"
19 #include "TSpectrum.h"
20 #include "TString.h"
21 #include "TLine.h"
22 #include "TText.h"
23 #include "TROOT.h"
24 #include "TList.h"
25 #include "TPolyMarker.h"
26 #include "StBichsel/Bichsel.h"
27 #include "BetheBloch.h"
28 #include "TDirIter.h"
29 #include "TTreeIter.h"
30 #else
31 class TSystem;
32 class TMath;
33 class TH1;
34 class TH2;
35 class TH3;
36 class TProfile;
37 class TStyle;
38 class TF1;
39 class TTree;
40 class TChain;
41 class TFile;
42 class TNtuple;
43 class TCanvas;
44 class TMinuit;
45 class TSpectrum;
46 class TString;
47 class TLine;
48 class TText;
49 class TROOT;
50 class TList;
51 class TPolyMarker;
52 class Bichsel;
53 class BetheBloch;
54 class TDirIter;
55 class TTreeIter;
56 #endif
57 #include "StLaserAnalysisMaker/laserino.h"
58 void CheckMirrors(Int_t sector, Int_t var = 1) {
59  static const Char_t *Names[3] = {"X","Y","Z"};
60  TTree *laser = (TTree *) gDirectory->Get("laser");
61  if (! laser ) return;
62  TCanvas *c = new TCanvas(Form("Sector%02i%s",sector,Names[var-1]),Form("Sector%02i%s",sector,Names[var-1]),600,800);
63  c->Divide(2,3);
64  for (Int_t Bundle = 1; Bundle <= 6; Bundle++) {
65  c->cd(Bundle);
66  laser->Draw(Form("fTracks.Laser.XyzB.mX%i:fTracks.Laser.Mirror",var),
67  Form("fTracks.Laser.Sector==%i&&fTracks.Laser.Bundle==%i",sector,Bundle),"colz",10000);
68  }
69 }
70 //________________________________________________________________________________
71 void CheckMirrors(const Char_t *files = "./laser_8102110.root") {
72  TDirIter Dir(files);
73  TTreeIter iter("laser");
74  iter.AddFile(files);
75  const Int_t& fNtrack = iter("fNtrack");
76  const UShort_t*& fTracks_mNumberOfFitPointsTpc = iter("fTracks.mNumberOfFitPointsTpc");
77  const Double32_t*& fTracks_fpTInv = iter("fTracks.fpTInv");
78  const Int_t*& fTracks_Flag = iter("fTracks.Flag");
79  const Double_t *&fTracks_dirPU_mX1 = iter("fTracks.dirPU.mX1");
80  const Double_t*& fTracks_dirPM_mX1 = iter("fTracks.dirPM.mX1");
81  const Double_t*& fTracks_dirPM_mX2 = iter("fTracks.dirPM.mX2");
82  const Int_t *&fTracks_Laser_Sector = iter("fTracks.Laser.Sector");
83  const Int_t *&fTracks_Laser_Bundle = iter("fTracks.Laser.Bundle");
84  const Int_t *&fTracks_Laser_Mirror = iter("fTracks.Laser.Mirror");
85  const Double_t*& fTracks_XyzPB_mX1 = iter("fTracks.XyzPB.mX1");
86  const Double_t*& fTracks_XyzPB_mX2 = iter("fTracks.XyzPB.mX2");
87  const Double_t*& fTracks_XyzPB_mX3 = iter("fTracks.XyzPB.mX3");
88  const Double_t*& fTracks_XyzPL_mX3 = iter("fTracks.XyzPL.mX3");
89  const Double_t*& fTracks_Laser_XyzB_mX1 = iter("fTracks.Laser.XyzB.mX1");
90  const Double_t*& fTracks_Laser_XyzB_mX2 = iter("fTracks.Laser.XyzB.mX2");
91  const Double_t*& fTracks_Laser_XyzB_mX3 = iter("fTracks.Laser.XyzB.mX3");
92  const Int_t NS = 12;
93  const Int_t NB = 6;
94  const Int_t NM = 7;
95  TH2D *X = new TH2D("X","dX versus mirror, bundle, sector",NS*NB*NM,0,NS*NB*NM,500,-.5,0.5);
96  TH2D *Y = new TH2D("Y","dY versus mirror, bundle, sector",NS*NB*NM,0,NS*NB*NM,500,-.5,0.5);
97  TH2D *Z = new TH2D("Z","dZ versus mirror, bundle, sector",NS*NB*NM,0,NS*NB*NM,500,-.5,0.5);
98  while (iter.Next()) {
99  for (Int_t k = 0; k < fNtrack; k++) {
100  if (fTracks_Flag[k] < 2) continue;
101  if (fTracks_dirPU_mX1[k] > -0.5) continue;
102  if (fTracks_mNumberOfFitPointsTpc[k] < 35) continue;
103  static const Double_t pTInv0 = 4.78815e-03;
104  static const Double_t DpTInv0 = 9.75313e-03;
105  if (TMath::Abs(fTracks_fpTInv[k] - pTInv0) > 3*DpTInv0) continue;
106  Int_t s = fTracks_Laser_Sector[k]/2 - 1;
107  if (s < 0 || s >= NS) continue;
108  Int_t b = fTracks_Laser_Bundle[k] - 1;
109  if (b < 0 || b >= NB) continue;
110  Int_t m = fTracks_Laser_Mirror[k] - 1;
111  if (m < 0 || m >= NM) continue;
112  Double_t x = fTracks_XyzPB_mX1[k] - fTracks_Laser_XyzB_mX1[k];
113  Double_t y = fTracks_XyzPB_mX2[k] - fTracks_Laser_XyzB_mX2[k];
114  Double_t z = fTracks_XyzPB_mX3[k] - fTracks_Laser_XyzB_mX3[k];
115  //laser->Draw("fTracks.XyzPB.mX3-fTracks.Laser.XyzB.mX3:fTracks.XyzPL.mX3>>WE6","fTracks.Flag>1&&abs(fTracks.XyzPM.mX3-6.5)<1&&fTracks.Laser.Sector!=6","colz")
116 #if 0
117  if (s < 6) z -= 6.39057;
118  else z -= 6.59919;
119  z -= 3.29319000000000009e-04*fTracks_XyzPL_mX3[k];
120 #else
121  if (s < 6) z -= 6.45043e+00+1.26835e-04*fTracks_XyzPL_mX3[k];
122  else z -= 6.60302e+00+3.93078e-04*fTracks_XyzPL_mX3[k];
123 #endif
124  Double_t indx = m + NM*(b + NB*s) + 0.5;
125  X->Fill(indx,x);
126  Y->Fill(indx,y);
127  Z->Fill(indx,z);
128  }
129  }
130  X->FitSlicesY();
131  Y->FitSlicesY();
132  Z->FitSlicesY();
133  TH1D *Fit[3][3];
134  const Char_t *xyz[3] = {"X","Y","Z"};
135  const Char_t *res[3] = {"1","2","chi2"};
136  for (Int_t i = 0; i < 3; i++)
137  for (Int_t j = 0; j < 3; j++) Fit[i][j] = (TH1D*) gDirectory->Get(Form("%s_%s",xyz[i],res[j]));
138  ofstream out;
139  TString fOut(files);
140  fOut,ReplaceAll("*","");
141  fOut.ReplaceAll(".root","");
142  fOut,ReplaceAll(".","");
143  fOut.ReplaceAll("/","");
144  fOut += ".data";
145  cout << "Create " << fOut << endl;
146  out.open(fOut.Data());
147  for (Int_t s = 0; s < 12; s++) {
148  for (Int_t b = 0; b < 6; b++) {
149  for (Int_t m = 0; m < 7; m++) {
150  Int_t bin = m + NM*(b + NB*s) + 1;
151  Double_t dxyz[3] = {0,0,0};
152  Double_t ddxyz[3] = {0,0,0};
153  Int_t iflag[3] = {0,0,0};
154  for (Int_t i = 0; i < 3; i++) {
155  if (Fit[i][0]->GetBinError(bin) > 0 && Fit[i][0]->GetBinError(bin) < 0.1) {
156  if (TMath::Abs(Fit[i][0]->GetBinContent(bin)) > Fit[i][0]->GetBinError(bin)) iflag[i] = 1;
157  dxyz[i] = Fit[i][0]->GetBinContent(bin);
158  ddxyz[i] = Fit[i][0]->GetBinError(bin);
159  }
160  }
161  if (iflag[0] + iflag[1] + iflag[2] || ddxyz[0] > 0 || dxyz[1] > 0 || dxyz[2] > 0)
162  cout << Form("Sector: %2i",2*s + 2) << " Bundle: " << b+1 << " Mirror: " << m+1
163  << " <x> = " << Form("%7.4f",dxyz[0]) << " +/- " << Form("%7.4f",ddxyz[0])
164  << " <y> = " << Form("%7.4f",dxyz[1]) << " +/- " << Form("%7.4f",ddxyz[1])
165  << " <z> = " << Form("%7.4f",dxyz[2]) << " +/- " << Form("%7.4f",ddxyz[2])
166  << endl;
167  out << Form(",%7.4f,%7.4f",Mirrors[s][b][m].dX+dxyz[0],ddxyz[0])
168  << Form(",%7.4f,%7.4f",Mirrors[s][b][m].dY+dxyz[1],ddxyz[1])
169  << Form(",%7.4f,%7.4f",Mirrors[s][b][m].dZ+dxyz[2],ddxyz[2])
170  << "}";
171  if (m == 6) out << "}";
172  if (b == 5 && m == 6) out << "}";
173  out << ",";
174  out << "// S/B/M " << 2*(s+1) << "/" << b+1 << "/" << m+1 << endl;
175  }
176  }
177  }
178  out.close();
179 }
180 #if 0
181 //________________________________________________________________________________
182 void CheckMirrors(const Char_t *files = "./laser_8102110.root") {
183  TDirIter Dir(files);
184  TTreeIter iter("laser");
185  iter.AddFile(files);
186  const Int_t& fNtrack = iter("fNtrack");
187  const UShort_t*& fTracks_mNumberOfFitPointsTpc = iter("fTracks.mNumberOfFitPointsTpc");
188  const Double32_t*& fTracks_fpTInv = iter("fTracks.fpTInv");
189  const Int_t*& fTracks_Flag = iter("fTracks.Flag");
190  const Double_t *&fTracks_dirPU_mX1 = iter("fTracks.dirPU.mX1");
191  const Double_t*& fTracks_dirPM_mX1 = iter("fTracks.dirPM.mX1");
192  const Double_t*& fTracks_dirPM_mX2 = iter("fTracks.dirPM.mX2");
193  const Int_t *&fTracks_Laser_Sector = iter("fTracks.Laser.Sector");
194  const Int_t *&fTracks_Laser_Bundle = iter("fTracks.Laser.Bundle");
195  const Int_t *&fTracks_Laser_Mirror = iter("fTracks.Laser.Mirror");
196  const Double_t*& fTracks_XyzPB_mX1 = iter("fTracks.XyzPB.mX1");
197  const Double_t*& fTracks_XyzPB_mX2 = iter("fTracks.XyzPB.mX2");
198  const Double_t*& fTracks_XyzPB_mX3 = iter("fTracks.XyzPB.mX3");
199  const Double_t*& fTracks_Laser_XyzB_mX1 = iter("fTracks.Laser.XyzB.mX1");
200  const Double_t*& fTracks_Laser_XyzB_mX2 = iter("fTracks.Laser.XyzB.mX2");
201  const Double_t*& fTracks_Laser_XyzB_mX3 = iter("fTracks.Laser.XyzB.mX3");
202  Int_t N = 0;
203  struct Kxy_t {
204  Int_t N;
205  Double_t K, Kx, K2, K2x, Ky, x, y, x2, xy, y2, del, Kdel, z, z2;
206  };
207  Kxy_t Offset[12][6][8];
208  const Int_t NW = (sizeof(Kxy_t) - sizeof(Int_t))/sizeof(Double_t);
209  memset(Offset, 0, 12*6*8*sizeof(Kxy_t));
210  while (iter.Next()) {
211  // cout << " fEvtHdr_fRun " << fEvtHdr_fRun << " fEvtHdr_fEvtNum " << fEvtHdr_fEvtNum << " fNtrack " << fNtrack << endl;
212  N++;
213  // if (N > 10000) break;
214  for (Int_t k = 0; k < fNtrack; k++) {
215  if (fTracks_Flag[k] < 2) continue;
216  if (fTracks_dirPU_mX1[k] > -0.5) continue;
217  if (fTracks_mNumberOfFitPointsTpc[k] < 35) continue;
218  static const Double_t pTInv0 = 4.78815e-03;
219  static const Double_t DpTInv0 = 9.75313e-03;
220  if (TMath::Abs(fTracks_fpTInv[k] - pTInv0) > 3*DpTInv0) continue;
221  Int_t s = fTracks_Laser_Sector[k]/2 - 1;
222  if (s < 0 || s > 11) continue;
223  Int_t b = fTracks_Laser_Bundle[k] - 1;
224  if (b < 0 || b > 5) continue;
225  Int_t m = fTracks_Laser_Mirror[k] - 1;
226  if (m < 0 || m > 6) continue;
227  LOOP:
228  Kxy_t *mirror = &Offset[s][b][m];
229  Double_t K = fTracks_dirPM_mX2[k]/fTracks_dirPM_mX1[k];
230  Double_t x = fTracks_XyzPB_mX1[k] - fTracks_Laser_XyzB_mX1[k];
231  Double_t y = fTracks_XyzPB_mX2[k] - fTracks_Laser_XyzB_mX2[k];
232  Double_t z = fTracks_XyzPB_mX3[k] - fTracks_Laser_XyzB_mX3[k];
233  if (s < 6) z -= 6.46746;
234  else z -= 6.56439;
235  if (TMath::Abs(z) > 0.5) continue;
236  if (TMath::Abs(x) > 1 || TMath::Abs(y) > 1) continue;
237  Double_t del = y - K*x;
238  mirror->N++;
239  mirror->K += K;
240  mirror->Kx += x*K;
241  mirror->K2 += K*K;
242  mirror->K2x += K*K*x;
243  mirror->Ky += K*y;
244  mirror->x += x;
245  mirror->y += y;
246  mirror->x2 += x*x;
247  mirror->xy += x*y;
248  mirror->y2 += y*y;
249  mirror->del += del;
250  mirror->Kdel+= K*del;
251  mirror->z += z;
252  mirror->z2 += z*z;
253  if (m != 7) {
254  m = 7; goto LOOP;
255  }
256  }
257  }
258  ofstream out;
259  TString fOut = "LaserCorrection.data";
260  cout << "Create " << fOut << endl;
261  out.open(fOut.Data());
262  for (Int_t s = 0; s < 12; s++) {
263  for (Int_t b = 0; b < 6; b++) {
264  for (Int_t m = 0; m < 8; m++) {
265  Kxy_t *mirror = &Offset[s][b][m];
266  Double_t *xx = &mirror->K;
267  if (mirror->N > 2) for (Int_t i = 0; i < NW; i++) xx[i] /= mirror->N;
268  else for (Int_t i = 0; i < NW; i++) xx[i] = 0;
269  Double_t dx, dy, dz;
270  dx = dy = dz = 0;
271  if ( mirror->N > 2 ) {
272  dx = TMath::Sqrt(mirror->x2 - mirror->x* mirror->x);
273  dy = TMath::Sqrt(mirror->y2 - mirror->y* mirror->y);
274  dz = TMath::Sqrt(mirror->z2 - mirror->z* mirror->z);
275  cout << "Sector: " << 2*s + 2 << " Bundle: " << b+1 << " Mirror: " << m+1 << " N: " << Form("%6i",mirror->N)
276  << " <x> = " << Form("%7.4f",mirror->x) << " +/- " << Form("%7.4f",dx)
277  << " <y> = " << Form("%7.4f",mirror->y) << " +/- " << Form("%7.4f",dy)
278  << " <z> = " << Form("%7.4f",mirror->z) << " +/- " << Form("%7.4f",dz)
279  << endl;
280  }
281  if (m < 7) {
282  out << Form(",%7.4f,%7.4f",Mirrors[s][b][m].dX+mirror->x,dx)
283  << Form(",%7.4f,%7.4f",Mirrors[s][b][m].dY+mirror->y,dy)
284  << Form(",%7.4f,%7.4f",Mirrors[s][b][m].dZ+mirror->z,dz)
285  << "}";
286  if (m == 6) out << "}";
287  if (b == 5 && m == 6) out << "}";
288  out << ",";
289  out << "// S/B/M " << 2*(s+1) << "/" << b+1 << "/" << m+1 << " => " << mirror->N << endl;
290  }
291 #if 0
292  else {
293  Double_t K = mirror->K;
294  Double_t K2 = mirror->K2;
295  Double_t del = mirror->del;
296  Double_t Kdel = mirror->Kdel;
297  Double_t det = K2 - K*K;
298  cout << "======================================================" << endl;
299  cout << "Sector: " << 2*s + 2 << " Bundle: " << b+1 << " Mirror: " << m+1 << " N: " << mirror->N << " det: " << det << endl;
300  if (TMath::Abs(det) < 1e-7) continue;
301  Double_t x0 = - (Kdel - K*del)/det;
302  Double_t y0 = K*x0 + del;
303  cout << "Sector: " << 2*s + 2 << " Bundle: " << b+1 << " Mirror: " << m+1
304  << " x0/y0 = " << x0 << " / " << y0 << endl;
305  cout << "======================================================" << endl;
306  }
307 #endif
308  }
309  }
310  }
311  out.close();
312 }
313 //________________________________________________________________________________
314 void dZ(const Char_t *files = "./laser_8102110.root") {
315  TDirIter Dir(files);
316  TTreeIter iter("laser");
317  iter.AddFile(files);
318  const Int_t &fNtrack = iter("fNtrack");
319  const Int_t *&fTracks_Flag = iter("fTracks.Flag");
320  const Double_t *&fTracks_dU_mX1 = iter("fTracks.dU.mX1");
321  const Double_t *&fTracks_dU_mX2 = iter("fTracks.dU.mX2");
322  const Double_t *&fTracks_dU_mX3 = iter("fTracks.dU.mX3");
323  const Int_t *&fTracks_Laser_Sector = iter("fTracks.Laser.Sector");
324  const Int_t *&fTracks_Laser_Bundle = iter("fTracks.Laser.Bundle");
325  const Int_t *&fTracks_Laser_Mirror = iter("fTracks.Laser.Mirror");
326  const Double_t *&fTracks_Laser_XyzU_mX3 = iter("fTracks.Laser.XyzU.mX3");
327  //laser->Draw("fTracks.XyzPU.mX3:fTracks.dU.mX3","fTracks.Flag>1&&abs(fTracks.dU.mX1)<.1&&abs(fTracks.dU.mX2)<.1&&fTracks.Laser.Sector==2&&fTracks.Laser.Bundle==4&&fTracks.Laser.Mirror==3","colz")
328  while (iter.Next()) {
329  for (Int_t k = 0; k < fNtrack; k++) {
330  if (fTracks_Flag[k] < 2) continue;
331  if (TMath::Abs(fTracks_dU_mX1[k]) > 0.1 || TMath::Abs(fTracks_dU_mX2[k]) > 0.1) continue;
332  }
333 }
334 #endif
Generated by   doxygen 1.8.5