| 1 |
#if !defined(__CINT__) || defined(__MAKECINT__) |
| 2 |
|
| 3 |
|
| 4 |
#include <TString.h> |
| 5 |
#include <TH1F.h> |
| 6 |
#include <TH2F.h> |
| 7 |
#include <TMath.h> |
| 8 |
#include <TLine.h> |
| 9 |
#include <TPolyMarker.h> |
| 10 |
#include <TSelector.h> |
| 11 |
#include <TFile.h> |
| 12 |
|
| 13 |
#include <stdlib.h> |
| 14 |
#include <iostream> |
| 15 |
using namespace std; |
| 16 |
|
| 17 |
#include <PamLevel2.h> |
| 18 |
|
| 19 |
#endif |
| 20 |
//=============================================================================== |
| 21 |
// |
| 22 |
// |
| 23 |
// |
| 24 |
// |
| 25 |
//=============================================================================== |
| 26 |
class CaloAxis { |
| 27 |
|
| 28 |
public: |
| 29 |
vector<float> x;// independent coordinates (z) |
| 30 |
vector<float> y;// dependente coordinate (x/y) |
| 31 |
vector<float> w;// fit weight |
| 32 |
vector<float> q;// signal |
| 33 |
|
| 34 |
Float_t par[2]; |
| 35 |
Float_t covpar[2][2]; |
| 36 |
|
| 37 |
TPolyMarker *sel; |
| 38 |
TLine *axis; |
| 39 |
|
| 40 |
CaloAxis(){ |
| 41 |
for(Int_t i=0; i<2;i++){ |
| 42 |
par[i]=0; |
| 43 |
for(Int_t j=0; j<2;j++)covpar[i][j]=0; |
| 44 |
}; |
| 45 |
sel = 0; |
| 46 |
axis = 0; |
| 47 |
}; |
| 48 |
|
| 49 |
CaloAxis(Int_t nin, float* xin, float *yin){ |
| 50 |
for(Int_t i=0; i<nin; i++){ |
| 51 |
x.push_back(*xin++); |
| 52 |
y.push_back(*yin++); |
| 53 |
w.push_back(1.); |
| 54 |
q.push_back(0.); |
| 55 |
}; |
| 56 |
if(nin>2)Fit(); |
| 57 |
sel = 0; |
| 58 |
axis = 0; |
| 59 |
}; |
| 60 |
CaloAxis(Int_t nin, float* xin, float *yin, float* qin){ |
| 61 |
for(Int_t i=0; i<nin; i++){ |
| 62 |
x.push_back(*xin++); |
| 63 |
y.push_back(*yin++); |
| 64 |
q.push_back(*qin++); |
| 65 |
w.push_back(1.); |
| 66 |
}; |
| 67 |
if(nin>2)Fit(); |
| 68 |
sel = 0; |
| 69 |
axis = 0; |
| 70 |
}; |
| 71 |
CaloAxis(Int_t nin, float* xin, float *yin, float* qin, float* win){ |
| 72 |
for(Int_t i=0; i<nin; i++){ |
| 73 |
x.push_back(*xin++); |
| 74 |
y.push_back(*yin++); |
| 75 |
w.push_back(*win++); |
| 76 |
q.push_back(*qin++); |
| 77 |
}; |
| 78 |
if(nin>2)Fit(); |
| 79 |
sel = 0; |
| 80 |
axis = 0; |
| 81 |
}; |
| 82 |
|
| 83 |
void Add(float xin, float yin){ |
| 84 |
x.push_back(xin); |
| 85 |
y.push_back(yin); |
| 86 |
w.push_back(1.); |
| 87 |
q.push_back(0.); |
| 88 |
Int_t nin = x.size(); |
| 89 |
if(nin>2)Fit(); |
| 90 |
}; |
| 91 |
void Add(float xin, float yin, float qin){ |
| 92 |
x.push_back(xin); |
| 93 |
y.push_back(yin); |
| 94 |
q.push_back(qin); |
| 95 |
w.push_back(1.); |
| 96 |
Int_t nin = x.size(); |
| 97 |
|
| 98 |
if(nin>2)Fit(); |
| 99 |
}; |
| 100 |
void Add(float xin, float yin, float qin, float win){ |
| 101 |
x.push_back(xin); |
| 102 |
y.push_back(yin); |
| 103 |
w.push_back(win); |
| 104 |
q.push_back(qin); |
| 105 |
Int_t nin = x.size(); |
| 106 |
if(nin>2)Fit(); |
| 107 |
}; |
| 108 |
|
| 109 |
Int_t Fit(){ |
| 110 |
|
| 111 |
Float_t SSS = 0; |
| 112 |
Float_t SSX = 0; |
| 113 |
Float_t SXX = 0; |
| 114 |
|
| 115 |
Float_t SSY = 0; |
| 116 |
Float_t SXY = 0; |
| 117 |
|
| 118 |
Int_t np=0; |
| 119 |
for(Int_t i=0; i<(int)(x.size()); i++){ |
| 120 |
SSS += w[i]*w[i]; |
| 121 |
SSX += x[i]*w[i]*w[i]; |
| 122 |
SXX += x[i]*x[i]*w[i]*w[i]; |
| 123 |
SSY += y[i]*w[i]*w[i]; |
| 124 |
SXY += x[i]*y[i]*w[i]*w[i]; |
| 125 |
if(w[i])np++; |
| 126 |
} |
| 127 |
Float_t det = SSS*SXX-SSX*SSX; |
| 128 |
if(det==0)return 0; |
| 129 |
if(np<3)return 0; |
| 130 |
// cout << np << "points fitted -- ok "<<endl; |
| 131 |
par[0] = (SSY*SXX-SXY*SSX)/det; |
| 132 |
par[1] = (SXY*SSS-SSY*SSX)/det; |
| 133 |
return 1; |
| 134 |
}; |
| 135 |
|
| 136 |
void Reset(){ |
| 137 |
for(Int_t i=0; i<2;i++){ |
| 138 |
par[i]=0; |
| 139 |
for(Int_t j=0; j<2;j++)covpar[i][j]=0; |
| 140 |
}; |
| 141 |
x.clear(); |
| 142 |
y.clear(); |
| 143 |
w.clear(); |
| 144 |
q.clear(); |
| 145 |
if(sel)sel->Delete(); |
| 146 |
if(axis)axis->Delete(); |
| 147 |
|
| 148 |
} |
| 149 |
|
| 150 |
void Delete(){ |
| 151 |
Reset(); |
| 152 |
// delete x; |
| 153 |
// delete y; |
| 154 |
// delete w; |
| 155 |
// delete q; |
| 156 |
} |
| 157 |
|
| 158 |
void Print(){ |
| 159 |
for(Int_t i=0; i<(int)(x.size()); i++)if(w[i])cout << x[i] << " - "<<y[i]<<endl; |
| 160 |
} |
| 161 |
|
| 162 |
float GetPar(Int_t i){return par[i];}; |
| 163 |
|
| 164 |
int GetN(){return (int)(x.size());}; |
| 165 |
|
| 166 |
float* GetX(){ |
| 167 |
float *xout = new float[(int)(x.size())]; |
| 168 |
for(Int_t i=0; i<(int)(x.size()); i++)xout[i]=x[i]; |
| 169 |
return xout; |
| 170 |
} |
| 171 |
|
| 172 |
float* GetY(){ |
| 173 |
float *yout = new float[(int)(y.size())]; |
| 174 |
for(Int_t i=0; i<(int)(y.size()); i++) yout[i]=y[i]; |
| 175 |
return yout; |
| 176 |
} |
| 177 |
|
| 178 |
float* GetQ(){ |
| 179 |
float *qout = new float[(int)(q.size())]; |
| 180 |
for(int i=0; i<(int)(q.size()); i++) qout[i]=q[i]; |
| 181 |
return qout; |
| 182 |
} |
| 183 |
|
| 184 |
float GetChi2(){ |
| 185 |
float chi2=0; |
| 186 |
int nchi=0; |
| 187 |
for(Int_t i=0; i<(int)(y.size()); i++){ |
| 188 |
chi2 += w[i]*w[i]*(y[i]-par[0]-par[1]*x[i])*(y[i]-par[0]-par[1]*x[i]); |
| 189 |
nchi += (int)w[i]; |
| 190 |
}; |
| 191 |
if(nchi<3)return -999; |
| 192 |
chi2=chi2/(nchi-2); |
| 193 |
return chi2; |
| 194 |
} |
| 195 |
|
| 196 |
float GetQaxis(){ |
| 197 |
float qaxis=0; |
| 198 |
for(Int_t i=0; i<(int)(q.size()); i++) qaxis +=q[i]; |
| 199 |
// for(Int_t i=0; i<(int)(q.size()); i++) cout<<q[i]<<endl; |
| 200 |
return qaxis; |
| 201 |
} |
| 202 |
|
| 203 |
float GetQstrip(){ |
| 204 |
|
| 205 |
float qaxis=0; |
| 206 |
for(Int_t i=0; i<(int)(q.size()); i++) qaxis +=q[i]; |
| 207 |
if(q.size()==0) return 0.; |
| 208 |
return qaxis/q.size(); |
| 209 |
} |
| 210 |
|
| 211 |
float GetXmin(){ |
| 212 |
float zmin=9999; |
| 213 |
for(Int_t i=0; i<(int)(x.size()); i++) if(x[i]<zmin) zmin = x[i]; |
| 214 |
return zmin; |
| 215 |
} |
| 216 |
float GetXmax(){ |
| 217 |
float zmax=-9999; |
| 218 |
for(Int_t i=0; i<(int)(x.size()); i++) if(x[i]>zmax) zmax = x[i]; |
| 219 |
return zmax; |
| 220 |
|
| 221 |
} |
| 222 |
|
| 223 |
|
| 224 |
void Draw(){ |
| 225 |
|
| 226 |
if(GetN()<3)return; |
| 227 |
|
| 228 |
sel = new TPolyMarker(GetN(),GetY(),GetX()); |
| 229 |
sel->SetMarkerSize(0.3); |
| 230 |
sel->SetMarkerColor(5); |
| 231 |
|
| 232 |
axis = new TLine( par[0]+GetXmin()*par[1], GetXmin(), par[0]+GetXmax()*par[1], GetXmax() ); |
| 233 |
axis->SetLineWidth(1); |
| 234 |
axis->SetLineColor(5); |
| 235 |
|
| 236 |
sel->Draw(); |
| 237 |
axis->Draw(); |
| 238 |
} |
| 239 |
|
| 240 |
Int_t FitAxis( CaloLevel1* calo , Int_t view , Float_t toll ){ |
| 241 |
|
| 242 |
CaloStrip* cstrip = new CaloStrip(calo); |
| 243 |
|
| 244 |
// --------------------------- |
| 245 |
// first : all hits included |
| 246 |
// --------------------------- |
| 247 |
// CaloAxis *axis = new CaloAxis(); |
| 248 |
for(int ih=0;ih<calo->istrip;ih++){ |
| 249 |
int iv=-1; |
| 250 |
int ip=-1; |
| 251 |
int is=-1; |
| 252 |
calo->DecodeEstrip(ih,iv,ip,is); |
| 253 |
if(iv==view){ |
| 254 |
// cout<<ip<<" "<<is<<" "<< calo->DecodeEstrip(ih,iv,ip,is)<<" "<<endl; |
| 255 |
// float x = GetCoordXY(view,ip,is) ; |
| 256 |
// float z = GetCoordZ(view,ip) ; |
| 257 |
// float q = calo->GetEstrip(view,ip,is); |
| 258 |
cstrip->Set(view,ip,is); |
| 259 |
float x = 0; |
| 260 |
if( !view ) x = cstrip->GetX(); |
| 261 |
else x = cstrip->GetY(); |
| 262 |
float z = cstrip->GetZ(); |
| 263 |
float q = cstrip->GetE(); |
| 264 |
// if(view==0 && ip==21)q=0; |
| 265 |
if(q>0)this->Add(z,x,q); |
| 266 |
} |
| 267 |
|
| 268 |
} |
| 269 |
// cout << "(1) hits "<<this->GetN()<<endl; |
| 270 |
if(this->GetN()<3){ |
| 271 |
delete cstrip; |
| 272 |
return 0; |
| 273 |
} |
| 274 |
float P0; |
| 275 |
float P1; |
| 276 |
int dof; |
| 277 |
bool missing = false; |
| 278 |
int nmissing = 0; |
| 279 |
|
| 280 |
// --------------------------- |
| 281 |
// second : iteration |
| 282 |
// --------------------------- |
| 283 |
int niter=0; // number of iterations |
| 284 |
do{ |
| 285 |
niter++; |
| 286 |
float ttoll = toll; //tolerance (cm) |
| 287 |
if(niter==1)ttoll = 10.*toll; |
| 288 |
dof = this->GetN(); //previous fit |
| 289 |
P0 = this->par[0]; |
| 290 |
P1 = this->par[1]; |
| 291 |
this->Reset(); //fit reset |
| 292 |
// cout << "previous fit "<<dof<<" "<<P0<<" "<<P1<<endl; |
| 293 |
// cout << "(reset) hits "<<this->GetN()<<endl; |
| 294 |
missing = false; |
| 295 |
nmissing = 0; |
| 296 |
for(Int_t ip=0; ip<22; ip++){ |
| 297 |
// break the track if more than 2 planes are missing |
| 298 |
if(missing)nmissing++; |
| 299 |
if(nmissing==3)break; |
| 300 |
missing=true; |
| 301 |
// float z = GetCoordZ(view,ip) ; |
| 302 |
for(Int_t is=0; is<96; is++){ |
| 303 |
|
| 304 |
cstrip->Set(view,ip,is); |
| 305 |
float x = 0; |
| 306 |
if( !view ) x = cstrip->GetX(); |
| 307 |
else x = cstrip->GetY(); |
| 308 |
float z = cstrip->GetZ(); |
| 309 |
float q = cstrip->GetE(); |
| 310 |
|
| 311 |
if( q ){ |
| 312 |
// float x = GetCoordXY(view,ip,is) ; |
| 313 |
// float q = calo->GetEstrip(view,ip,is); |
| 314 |
// if(view==0 && ip==21)q=0; |
| 315 |
float d = fabs(x-P0-P1*z)/sqrt(P1*P1+1); |
| 316 |
// cout << d << " "<<ttoll<<endl; |
| 317 |
// if( q>0 && d < ttoll ){ |
| 318 |
if( d < ttoll ){ |
| 319 |
this->Add(z,x,q); |
| 320 |
// cout << "(-) hits "<<this->GetN()<<endl; |
| 321 |
missing=false; |
| 322 |
nmissing=0; |
| 323 |
} |
| 324 |
} |
| 325 |
} |
| 326 |
} |
| 327 |
// cout << "(2) hits "<<this->GetN()<<endl; |
| 328 |
if(niter==100)break; |
| 329 |
if(this->GetN()==0)break; |
| 330 |
}while(niter==1 || this->GetN() != dof); |
| 331 |
if(this->GetN()<3){ |
| 332 |
delete cstrip; |
| 333 |
return 0; |
| 334 |
} |
| 335 |
// --------------------------------------------- |
| 336 |
// third : selecting only closest strip-clusters |
| 337 |
// --------------------------------------------- |
| 338 |
P0 = this->par[0]; |
| 339 |
P1 = this->par[1]; |
| 340 |
this->Reset(); |
| 341 |
missing = false; |
| 342 |
nmissing = 0; |
| 343 |
for(Int_t ip=0; ip<22; ip++){ |
| 344 |
if(missing)nmissing++; |
| 345 |
if(nmissing==3)break; |
| 346 |
missing=true; |
| 347 |
int isclose = -1; |
| 348 |
float dmin = 999; |
| 349 |
// float z = GetCoordZ(view,ip) ; |
| 350 |
//search for closest hit strip |
| 351 |
for(Int_t is=0; is<96; is++){ |
| 352 |
cstrip->Set(view,ip,is); |
| 353 |
float x = 0; |
| 354 |
if( !view ) x = cstrip->GetX(); |
| 355 |
else x = cstrip->GetY(); |
| 356 |
float z = cstrip->GetZ(); |
| 357 |
float q = cstrip->GetE(); |
| 358 |
if( q ){ |
| 359 |
// float x = GetCoordXY(view,ip,is) ; |
| 360 |
float d = fabs(x-P0-P1*z)/sqrt(P1*P1+1); |
| 361 |
// cout << x <<" "<<z<<" "<<P0+P1*z<<" "<<d << endl; |
| 362 |
if(d<dmin && d<toll){ |
| 363 |
isclose =is; |
| 364 |
dmin =d; |
| 365 |
} |
| 366 |
} |
| 367 |
} |
| 368 |
if( isclose != -1 ){ |
| 369 |
Int_t isfrom = isclose-1; |
| 370 |
if(isfrom<0) isfrom=0; |
| 371 |
Int_t isto = isclose+2; |
| 372 |
if(isto>96) isto=96; |
| 373 |
for(Int_t is=isfrom; is< isto; is++){ |
| 374 |
cstrip->Set(view,ip,is); |
| 375 |
float x = 0; |
| 376 |
if( !view ) x = cstrip->GetX(); |
| 377 |
else x = cstrip->GetY(); |
| 378 |
float z = cstrip->GetZ(); |
| 379 |
float q = cstrip->GetE(); |
| 380 |
// float x = GetCoordXY(view,ip,is) ; |
| 381 |
// float q = calo->GetEstrip(view,ip,is); |
| 382 |
// if(view==0 && ip==21)q = 0; //mask last plane |
| 383 |
if(q){ |
| 384 |
this->Add(z,x,q); |
| 385 |
missing = false; |
| 386 |
nmissing = 0; |
| 387 |
} |
| 388 |
} |
| 389 |
} |
| 390 |
} |
| 391 |
|
| 392 |
// cout << " n. hits :"<<this->GetN()<<endl; |
| 393 |
// cout << " P0 P1 :"<<this->par[0]<< " "<<this->par[1]<<endl; |
| 394 |
// cout << " Chi2 :"<<this->GetChi2()<<endl; |
| 395 |
// cout << " Q axis :"<<this->GetQaxis()<<endl; |
| 396 |
// cout << " Q strip :"<<this->GetQstrip()<<endl; |
| 397 |
delete cstrip; |
| 398 |
return 1; |
| 399 |
} |
| 400 |
}; |
| 401 |
//=============================================================================== |
| 402 |
// |
| 403 |
// |
| 404 |
// |
| 405 |
// |
| 406 |
//=============================================================================== |
Parent Directory
| 
Revision Log