00001 #include <stdlib.h>
00002 #include <stdio.h>
00003 #include <assert.h>
00004 #include "StiDebug.h"
00005 #include "TMath.h"
00006 #include "TROOT.h"
00007 #include "TCanvas.h"
00008 #include "TGraph.h"
00009 #include "TSystem.h"
00010 #include "TObjArray.h"
00011 #include "TH1.h"
00012 #include "Sti/StiKalmanTrack.h"
00013 #include "Sti/StiKalmanTrackNode.h"
00014 static int myReady=0;
00015 TObjArray *StiDebug::mgHist=0;
00016 TObjArray *StiDebug::mgTally=0;
00017
00018
00019 void StiDebug::Init()
00020 {
00021 if (gROOT->IsBatch()) return;
00022 }
00023
00024 int StiDebug::iFlag(const char *flagName,int dflt)
00025 {
00026 const char *val = gSystem->Getenv(flagName);
00027 if (!val) return dflt;
00028 printf("\nStiDebug::iFlag: %s = %s\n\n",flagName,val);
00029
00030 return atoi(val);
00031 }
00032
00033
00034 double StiDebug::dFlag(const char *flagName, double dflt)
00035 {
00036 const char *val = gSystem->Getenv(flagName);
00037 if (!val) return dflt;
00038 printf("\nStiDebug::dFlag: %s = %s\n\n",flagName,val);
00039
00040 return atof(val);
00041 }
00042
00043 int StiDebug::tally(const char *name,int val)
00044 {
00045 if (!mgTally) mgTally = new TObjArray();
00046 TObject *to = mgTally->FindObject(name);
00047 if (!to) mgTally->Add((to=new TNamed(name,"")));
00048 int tally = to->GetUniqueID();
00049 tally +=val;
00050 to->SetUniqueID(tally);
00051 return tally;
00052 }
00053
00054 void StiDebug::hist(const char *name,double val)
00055 {
00056 if (!mgHist) return;
00057 TH1 *hist = (TH1*)mgHist->FindObject(name);
00058 if (!hist) return;
00059 hist->Fill(val);
00060 }
00061
00062 void StiDebug::Finish()
00063 {
00064 printf("\n");
00065 for (int kase=0;kase<2;kase++) {
00066 switch(kase) {
00067 case 0:{
00068 if (!mgTally) break;
00069 int nT = mgTally->GetLast()+1;
00070 if (!nT) break;
00071 for (int iT=0;iT<nT;iT++) {
00072 TObject *to = mgTally->At(iT);
00073 if (!to) continue;
00074 printf("StiDebug::Tally %s=%d\n",to->GetName(),to->GetUniqueID());
00075 }
00076 break;}
00077
00078 case 1:{
00079 if (gROOT->IsBatch()) break;
00080 if (!mgHist) break;
00081 int nH = mgHist->GetLast()+1;
00082 if (!nH) break;
00083 int n = 0;
00084 TCanvas *tc=0;
00085 for (int iH=0;iH<nH;iH++) {
00086 TH1 *th = (TH1*)mgHist->At(iH);
00087 if (!th) continue;
00088 if((n%3)==0) {
00089 tc = new TCanvas("C1","",600,800);
00090 tc->Divide(1,3);
00091 }
00092 tc->cd((n%3)+1); th->Draw();
00093 tc->Modified() ; tc->Update();
00094 n++;
00095 }}
00096 }
00097 }
00098
00099 }
00100
00101 void StiDebug::Break(int kase)
00102 {
00103 static int myBreak=-2005;
00104 if (kase!=myBreak) return;
00105 printf("*** Break(%d) ***\n",kase);
00106 }
00107
00108 void StiDebug::FpeOn()
00109 {
00110 gSystem->SetFPEMask(kInvalid | kDivByZero | kOverflow );
00111 printf("*** Float Point Exception is ON ***\n");
00112
00113 }
00114
00115 void StiDebug::show(StiKalmanTrack *kt)
00116 {
00117
00118
00119
00120
00121
00122
00123 static TCanvas *myCanvas=0;
00124 static TGraph *graph[3][3] = {{0,0,0},{0,0,0},{0,0,0}};
00125
00126 float X[3][3][100],Y[3][3][100];
00127 int N[3][3];
00128
00129 if (!myCanvas) {myCanvas = new TCanvas("C1","",600,800);
00130 myCanvas->Divide(1,3);}
00131 if(kt) {
00132 double curv = 0,xPrev,yPrev; int nCurv = 0;
00133 for (int i=0;i<9;i++) {delete graph[0][i];graph[0][i]=0;}
00134 StiKTNBidirectionalIterator it;
00135 for (int ig=0;ig<3;ig++) {
00136 int n=0;
00137 double s = 0;
00138 xPrev = 3e33;
00139 for (it=kt->begin();it!=kt->end();++it) {
00140 StiKalmanTrackNode *node = &(*it);
00141 if (!node->isValid()) continue;
00142
00143 double xNode = node->x_g();
00144 double yNode = node->y_g();
00145 if (xPrev<3e33) {
00146 double ds = sqrt(pow(xNode-xPrev,2)+pow(yNode-yPrev,2));
00147 double si = 0.5*ds*curv; if (si>0.99) si=0.99;
00148 if (si>0.01) ds = 2*asin(si)/curv;
00149 s += ds;
00150 }
00151 xPrev = xNode;
00152 yPrev = yNode;
00153
00154 StiHit *hit = node->getHit();
00155 if (hit && node->getChi2()>1000) hit=0;
00156 if (ig==2 && !hit) continue;
00157 if (ig==0) { curv += node->getCurvature(); nCurv++; continue;}
00158 if (ig==1) {
00159 X[0][ig][n] = node->x_g();
00160 Y[0][ig][n] = node->y_g();
00161 Y[2][ig][n] = node->z_g();
00162 } else {
00163 X[0][ig][n] = hit->x_g();
00164 Y[0][ig][n] = hit->y_g();
00165 Y[2][ig][n] = hit->z_g();
00166 }
00167
00168 if (n) {
00169 float xh = X[0][ig][n]-X[0][ig][0];
00170 float yh = Y[0][ig][n]-Y[0][ig][0];
00171 float rh = xh*xh+yh*yh+1E-10;
00172 X[1][ig][n-1] = xh/rh;
00173 Y[1][ig][n-1] = yh/rh;
00174 }
00175 X[2][ig][n]=s;
00176 n++;
00177 }
00178 if (ig==0) { curv=fabs(curv)/nCurv; continue;}
00179 N[0][ig] = n;
00180 N[1][ig] = n-1;
00181 N[2][ig] = n;
00182 }
00183
00184 for (int ip=0;ip<3;ip++) {
00185 double xMin=999,xMax=-999,yMin=999,yMax=-999;
00186 for (int ig=1;ig<3;ig++) {
00187 for (int j=0;j<N[ip][ig];j++) {
00188 double x = X[ip][ig][j];
00189 if (xMin> x) xMin = x;
00190 if (xMax< x) xMax = x;
00191 double y = Y[ip][ig][j];
00192 if (yMin> y) yMin = y;
00193 if (yMax< y) yMax = y;
00194 }
00195 }
00196 X[ip][0][0] = xMin; Y[ip][0][0] = yMin;
00197 X[ip][0][1] = xMin; Y[ip][0][1] = yMax;
00198 X[ip][0][2] = xMax; Y[ip][0][2] = yMin;
00199 X[ip][0][3] = xMax; Y[ip][0][3] = yMax;
00200 N[ip][0] = 4;
00201 }
00202 static const char *opt[]={"AP","Same CP","Same *"};
00203 for (int ip=0;ip<3;ip++) {
00204 for (int ig =0;ig<3;ig++) {
00205 graph[ip][ig] = new TGraph(N[ip][ig] , X[ip][ig], Y[ip][ig]);
00206 if(ig==2) graph[ip][ig]->SetMarkerColor(kRed);
00207 myCanvas->cd(ip+1); graph[ip][ig]->Draw(opt[ig]);
00208 }
00209 }
00210
00211 }
00212
00213
00214 if (!myCanvas) return;
00215 myCanvas->Modified();
00216 myCanvas->Update();
00217 myReady = 2005;
00218 while(!gSystem->ProcessEvents()){};
00219 }
00220
00221 ClassImp(StiAux)
00222
00223 StiAux::StiAux():TDataSet("StiAux")
00224 {
00225 fN=0;
00226 }
00227
00228 StiAux_t* StiAux::Get(int id) const
00229 {
00230 if (id>fN) return 0;
00231 return ((StiAux_t*)(fArr.GetArray())) + id-1;
00232 }
00233
00234 int StiAux::AddIt(StiAux_t *add)
00235 {
00236 int n = fArr.GetSize();
00237 fN++;
00238 if (int(fN*sizeof(StiAux_t))>n) fArr.Set((n+100)*2);
00239 memcpy(Get(fN),add,sizeof(StiAux_t));
00240 return fN;
00241 }
00242
00243 void StiAux::PrintIt(int id)
00244 {
00245 static const char* tit[] = {
00246 "xnl[3]=","","", "xhl[3]=","","",
00247 "ca=", "rho=",
00248 "nYY=", "nZZ=",
00249 "hYY=", "hZZ=",
00250 "xng[3]=","","", "xhg[3]=","","",
00251 "psi=", "dip=","chi2=",0};
00252
00253 float *aux = (float*)Get(id);
00254 printf("%d4 - ",id);
00255 for (int i=0;tit[i];i++) { printf("%s%g ",tit[i],aux[i]);}
00256 printf("\n");
00257
00258 }
00259
00260
00261
