StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StiDebug.cxx
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <math.h>
4 #include <assert.h>
5 #include <map>
6 #include <string>
7 #include "TClass.h"
8 #include "TH1F.h"
9 #include "TH2F.h"
10 #include "TCanvas.h"
11 #include "TGraph.h"
12 #include "TSystem.h"
13 #include "TString.h"
14 
15 #include "StiDebug.h"
16 #include "TMath.h"
17 #include "TROOT.h"
18 #include "TObjArray.h"
19 #include "Sti/StiKalmanTrack.h"
20 #include "Sti/StiKalmanTrackNode.h"
21 static int myReady=0;
22 typedef std::map<std::string, int> myTallyMap_t;
23 typedef myTallyMap_t::const_iterator myTallyIter_t;
24 static myTallyMap_t mgTally;
25 
26 
27 int StiDebug::mgDebug=1; //0=no debug, 1=Normal, 2=count is on
28 int StiDebug::mgRecov=1;
29 int StiDebug::mgCheck=1;
30 
31 typedef std::map<std::string, TH1*> myDebMap_t;
32 typedef myDebMap_t::const_iterator myDebIter_t;
33 static myDebMap_t myDebMap;
34 
35 typedef std::map<std::string, TCanvas*> myCanMap_t;
36 typedef myCanMap_t::const_iterator myCanIter_t;
37 myCanMap_t myCanMap;
38 //______________________________________________________________________________
39 int StiDebug::Debug()
40 {
41 return mgDebug;
42 }
43 //______________________________________________________________________________
44 void StiDebug::Init()
45 {
46  if (gROOT->IsBatch()) return;
47 }
48 //______________________________________________________________________________
49 int StiDebug::iFlag(const char *flagName,int dflt)
50 {
51  const char *val = gSystem->Getenv(flagName);
52  if (!val) return dflt;
53  printf("\nStiDebug::iFlag: %s = %s\n\n",flagName,val);
54 
55  return atoi(val);
56 }
57 
58 //______________________________________________________________________________
59 double StiDebug::dFlag(const char *flagName, double dflt)
60 {
61  const char *val = gSystem->Getenv(flagName);
62  if (!val) return dflt;
63  printf("\nStiDebug::dFlag: %s = %s\n\n",flagName,val);
64 
65  return atof(val);
66 }
67 //______________________________________________________________________________
68 void StiDebug::tally(const char *name,int val)
69 {
70  if (mgDebug<2) return;
71  mgTally[name] += val;
72 }
73 //______________________________________________________________________________
74 void StiDebug::Finish()
75 {
76 //VP Sumary();
77 }
78 //______________________________________________________________________________
79 void StiDebug::Break(int kase)
80 {
81 static int myBreak=-2005;
82 if (kase!=myBreak) return;
83 if (kase!=-2005)
84  printf("*** Break(%d) ***\n",kase);
85 }
86 //______________________________________________________________________________
87 void StiDebug::FpeOn()
88 {
89  gSystem->SetFPEMask(kInvalid | kDivByZero | kOverflow );
90  printf("*** Float Point Exception is ON ***\n");
91 
92 }
93 //______________________________________________________________________________
94 void StiDebug::show(StiKalmanTrack *kt)
95 {
96 // lev=0 draw all nodes
97 // lev=1 draw hit nodes
98 // lev=2 draw hits only
99 
100 
101 
102 static TCanvas *myCanvas=0;
103 static TGraph *graph[3][3] = {{0,0,0},{0,0,0},{0,0,0}};
104 // P G P G
105  float X[3][3][100],Y[3][3][100];
106  int N[3][3];
107 
108  if (!myCanvas) {myCanvas = new TCanvas("C1","",600,800);
109  myCanvas->Divide(1,3);}
110  if(kt) {
111  double curv = 0,xPrev=0,yPrev=0; int nCurv = 0;
112 
113  for (int i=0;i<9;i++) {delete graph[0][i];graph[0][i]=0;}
114 
115  StiKTNBidirectionalIterator it;
116  for (int ig=0;ig<3;ig++) {
117  int n=0;
118  double s = 0;
119  xPrev = 3e33;
120  for (it=kt->begin();it!=kt->end();++it) {
121  StiKalmanTrackNode *node = &(*it);
122  if (!node->isValid()) continue;
123 // S calculation common based on node x,y for both hit and node
124  double xNode = node->x_g();
125  double yNode = node->y_g();
126  if (xPrev<3e33) {
127  double ds = sqrt(pow(xNode-xPrev,2)+pow(yNode-yPrev,2));
128  double si = 0.5*ds*curv; if (si>0.99) si=0.99;
129  if (si>0.01) ds = 2*asin(si)/curv;
130  s += ds;
131  }
132  xPrev = xNode;
133  yPrev = yNode;
134 
135  StiHit *hit = node->getHit();
136  if (hit && node->getChi2()>1000) hit=0;
137  if (ig==2 && !hit) continue;
138  if (ig==0) { curv += node->getCurvature(); nCurv++; continue;}
139  if (ig==1) {//draw nodes
140  X[0][ig][n] = node->x_g();
141  Y[0][ig][n] = node->y_g();
142  Y[2][ig][n] = node->z_g();
143  } else {//draw hits only
144  X[0][ig][n] = hit->x_g();
145  Y[0][ig][n] = hit->y_g();
146  Y[2][ig][n] = hit->z_g();
147  }
148 
149  if (n) {
150  float xh = X[0][ig][n]-X[0][ig][0];
151  float yh = Y[0][ig][n]-Y[0][ig][0];
152  float rh = xh*xh+yh*yh+1E-10;
153  X[1][ig][n-1] = xh/rh;
154  Y[1][ig][n-1] = yh/rh;
155  }
156  X[2][ig][n]=s;
157  n++;
158  }//end for nodes
159  if (ig==0) { curv=fabs(curv)/nCurv; continue;}
160  N[0][ig] = n;
161  N[1][ig] = n-1;
162  N[2][ig] = n;
163  }//end for ig
164 
165  for (int ip=0;ip<3;ip++) {
166  double xMin=999,xMax=-999,yMin=999,yMax=-999;
167  for (int ig=1;ig<3;ig++) {
168  for (int j=0;j<N[ip][ig];j++) {
169  double x = X[ip][ig][j];
170  if (xMin> x) xMin = x;
171  if (xMax< x) xMax = x;
172  double y = Y[ip][ig][j];
173  if (yMin> y) yMin = y;
174  if (yMax< y) yMax = y;
175  }
176  }
177  X[ip][0][0] = xMin; Y[ip][0][0] = yMin;
178  X[ip][0][1] = xMin; Y[ip][0][1] = yMax;
179  X[ip][0][2] = xMax; Y[ip][0][2] = yMin;
180  X[ip][0][3] = xMax; Y[ip][0][3] = yMax;
181  N[ip][0] = 4;
182  }
183 static const char *opt[]={"AP","Same CP","Same *"};
184  for (int ip=0;ip<3;ip++) {
185  for (int ig =0;ig<3;ig++) {
186  graph[ip][ig] = new TGraph(N[ip][ig] , X[ip][ig], Y[ip][ig]);
187  if(ig==2) graph[ip][ig]->SetMarkerColor(kRed);
188  myCanvas->cd(ip+1); graph[ip][ig]->Draw(opt[ig]);
189  }//end for ig
190  }//end ip
191 
192  }//end if
193 
194 
195  if (!myCanvas) return;
196  myCanvas->Modified();
197  myCanvas->Update();
198  myReady = 2005;
199  while(!gSystem->ProcessEvents()){};
200 }
201 
202 //______________________________________________________________________________
203 StiAux::StiAux():TDataSet("StiAux")
204 {
205  fN=0;
206 }
207 //______________________________________________________________________________
208 StiAux_t* StiAux::Get(int id) const
209 {
210  if (id>fN) return 0;
211  return ((StiAux_t*)(fArr.GetArray())) + id-1;
212 }
213 //______________________________________________________________________________
214 int StiAux::AddIt(StiAux_t *add)
215 {
216  int n = fArr.GetSize();
217  fN++;
218  if (int(fN*sizeof(StiAux_t))>n) fArr.Set((n+100)*2);
219  memcpy(Get(fN),add,sizeof(StiAux_t));
220  return fN;
221 }
222 //______________________________________________________________________________
223 void StiAux::PrintIt(int id)
224 {
225 static const char* tit[] = {
226 "xnl[3]=","","", "xhl[3]=","","",
227 "ca=", "rho=",
228 "nYY=", "nZZ=",
229 "hYY=", "hZZ=",
230 "xng[3]=","","", "xhg[3]=","","",
231 "psi=", "dip=","chi2=",0};
232 
233  float *aux = (float*)Get(id);
234  printf("%d4 - ",id);
235  for (int i=0;tit[i];i++) { printf("%s%g ",tit[i],aux[i]);}
236  printf("\n");
237 
238 }
239 //______________________________________________________________________________
240 void StiDebug::Count(const char *key,double val)
241 {
242  if (mgDebug<2) return;
243  TH1 *&h = (TH1*&)myDebMap[key];
244  if (!h) { h = new TH1F(key,key,100,0.,0.);}
245  h->Fill(val);
246 }
247 //______________________________________________________________________________
248 void StiDebug::Count(const char *key,const char *val)
249 {
250  if (mgDebug<2) return;
251  TH1 *&h = (TH1*&)myDebMap[key];
252  if (!h) { h = new TH1F(key,key,0,0.,0.);}
253  h->Fill(val,1.);
254 }
255 //______________________________________________________________________________
256 void StiDebug::Count(const char *key,double val,double left,double rite)
257 {
258  if (mgDebug<2) return;
259  TH1 *&h = (TH1*&)myDebMap[key];
260  if (!h) { h = new TH1F(key,key,100,left,rite);}
261  h->Fill(val);
262 }
263 //______________________________________________________________________________
264 void StiDebug::Count(const char *key,double valx, double valy
265  ,double leftx,double ritex
266  ,double lefty,double ritey)
267 {
268  if (mgDebug<2) return;
269  TH1 *&h = (TH1*&)myDebMap[key];
270  if (!h) { h = new TH2F(key,key,100,leftx,ritex,100,lefty,ritey);}
271  h->Fill(valx,valy);
272 }
273 //______________________________________________________________________________
274 void StiDebug::Count(const char *key,double valx,double valy)
275 {
276  if (mgDebug<2) return;
277  TH1 *&h = (TH1*&)myDebMap[key];
278  if (!h) { h = new TH2F(key,key,100,0.,0.,100,0,0);}
279  h->Fill(valx,valy);
280 }
281 //______________________________________________________________________________
282 void StiDebug::Count(const char *key,const char *valx,double valy)
283 {
284  if (mgDebug<2) return;
285  TH2 *&h = (TH2*&)myDebMap[key];
286  if (!h) { h = new TH2F(key,key,100,0.,0.,100,0,0);}
287  h->Fill(valx,valy,1.);
288 }
289 //______________________________________________________________________________
290 void StiDebug::Sumary(int nHist)
291 {
292  if (!nHist) nHist=4;
293  printf("StiDebug::Sumary()\n");
294 
295  int nH = 0,n=0;
296  for (myTallyIter_t iter = mgTally.begin();iter != mgTally.end();++iter) {
297  int nn = (*iter).second; n++;
298  const char *kto = (*iter).first.c_str();
299  printf("%3d - Tally.%s = %d\n",n,kto,nn);
300  }
301 
302 
303 
304 
305  TH1 *H[10];
306  for (int numCha = 0; numCha<2; numCha++) {
307  nH = 0;n=0;
308  for (myDebIter_t iter = myDebMap.begin();iter != myDebMap.end();++iter) {
309  TH1 *h = (*iter).second; n++;
310  if (!h->IsA()->InheritsFrom("TH1")) continue;
311  if ( h->IsA()->InheritsFrom("TH2")) continue;
312 
313  int nEnt = h->GetEntries();
314  if (!nEnt) continue;
315  if ((h->GetXaxis()->GetLabels()==0) != (numCha == 0)) continue;
316  double mean = h->GetMean();
317  double rms = h->GetRMS();
318  printf("TH1 %2d - %12s:\t %5d %g(+-%g)\n",n,h->GetName(),nEnt,mean,rms);
319  if (nH==nHist) {Draw(nH,H);nH=0;}
320  if (numCha) h->LabelsOption(">V","X");
321  H[nH++] = h;
322  }
323  if (nH) Draw(nH,H);
324  }
325  for (myDebIter_t iter = myDebMap.begin();iter != myDebMap.end();++iter) {
326  TH1 *h = (*iter).second; n++;
327  if (!h->IsA()->InheritsFrom("TH2")) continue;
328  int nEnt = h->GetEntries();
329  if (!nEnt) continue;
330  double mean = h->GetMean();
331  double rms = h->GetRMS();
332  printf("TH2 %2d - %12s:\t %5d %g(+-%g)\n",n,h->GetName(),nEnt,mean,rms);
333  H[0]=h;Draw(1,H);
334  }
335 
336  if (!n) return;
337  while(!gSystem->ProcessEvents()){};
338 
339 }
340 
341 //______________________________________________________________________________
342 void StiDebug::Draw(int nH,TH1** H)
343 {
344 static int nCall=0; nCall++;
345 
346  TString ts("C_"); ts += H[0]->GetName();
347  TCanvas *&C = myCanMap[ts.Data()];
348  if (!C) C = new TCanvas(ts.Data(),ts.Data(),600,800);
349  C->Clear();C->Divide(1,nH);
350  for (int i=0;i<nH;i++) { C->cd(i+1); H[i]->Draw(); }
351 
352  C->Modified();C->Update();
353 }
354 
355 
StiKalmanTrack
Definition of Kalman Track.
Definition: StiKalmanTrack.h:82
StiHit
Definition: StiHit.h:51
StiHit::x_g
Float_t x_g() const
Return the global x, y, z values.
Definition: StiHit.h:73
StiAux_t
Definition: StiDebug.h:52
StiKalmanTrack.h
Definition of Kalman Track.
StiKalmanTrack::end
const StiKTNBidirectionalIterator & end() const
Definition: StiKalmanTrack.h:555
StiKalmanTrackNode::getCurvature
double getCurvature() const
Calculates and returns the tangent of the track pitch angle at this node.
Definition: StiKalmanTrackNode.h:126
StiKalmanTrack::begin
StiKTNBidirectionalIterator begin() const
Definition: StiKalmanTrack.h:540
Generated by   doxygen 1.8.5