#include <AlgFitTrack3.h>
Inheritance diagram for AlgFitTrack3:
|
|
Definition at line 48 of file AlgFitTrack3.cxx. References MSG. 00049 :fIterations(0),fTrackNum(0),fUWeightMatrix(0),fVWeightMatrix(0),fMatrixA(5,5),fMatrixC(5,1) 00050 { 00051 MSG("FitTrack3", Msg::kDebug) 00052 << "Starting AlgFitTrack3::AlgFitTrack3()" << endl; 00053 00054 char* outPath = getenv ("OUTPATH") ; 00055 if (outPath == NULL) { 00056 sprintf(outPath,"."); 00057 } 00058 char* runnum = getenv ("RUNNUM"); 00059 if(runnum == NULL) { 00060 sprintf(runnum,"3"); 00061 } 00062 char filename[180]; 00063 sprintf(filename,"%s/debugFitTrack%s.root",outPath,runnum); 00064 fFile = new TFile(filename,"RECREATE"); 00065 fFile->cd(); 00066 00067 fSillyTree = new TTree("sillyTree","sillyTree"); 00068 fSillyTree->Branch("z",&fS_z,"z/D"); 00069 fSillyTree->Branch("measuredU",&fS_measuredU,"measuredU/D"); 00070 fSillyTree->Branch("measuredV",&fS_measuredV,"measuredV/D"); 00071 fSillyTree->Branch("swamU",&fS_swamU,"swamU/D"); 00072 fSillyTree->Branch("swamV",&fS_swamV,"swamV/D"); 00073 fSillyTree->Branch("isU",&fS_isU,"isU/I"); 00074 fSillyTree->Branch("trackNum",&fTrackNum,"trackNum/I"); 00075 fSillyTree->Branch("iteration",&fIterations,"iteration/I"); 00076 fReasonsForNotFitting = new TH1F("reasons","Why didn't I Fit?",20,-0.5,9.5); 00077 00078 fReasonsTree = new TTree("reasonsTree","reasonsTree"); 00079 fReasonsTree->Branch("totalU",&fR_totalU,"totalU/I"); 00080 fReasonsTree->Branch("totalV",&fR_totalV,"totalV/I"); 00081 fReasonsTree->Branch("showerU",&fR_showerU,"showerU/I"); 00082 fReasonsTree->Branch("showerV",&fR_showerV,"showerV/I"); 00083 fReasonsTree->Branch("validU",&fR_validU,"validU/I"); 00084 fReasonsTree->Branch("validV",&fR_validV,"validV/I"); 00085 }
|
|
|
Definition at line 88 of file AlgFitTrack3.cxx. References fActualZ, fDistanceFromStart, fErrorOnMeasure, fFile, fMeasuredU, fMeasuredV, fRangeThisPlane, fReverseMapUIndex, fReverseMapVIndex, fSwamdU, fSwamdV, fSwamU, fSwamV, fTheta0i, fThicknessi, fTrkClstList, and MSG. 00089 {
00090 MSG("FitTrack3", Msg::kDebug)
00091 << "AlgFitTrack3::~AlgFitTrack3()" << endl;
00092
00093 if(fUWeightMatrix) delete fUWeightMatrix;
00094 if(fVWeightMatrix) delete fVWeightMatrix;
00095 fTheta0i.clear(); // mean square scattering angle.
00096 fThicknessi.clear(); // thickness passed through at plane i.
00097 fSwamU.clear();
00098 fSwamV.clear();
00099 fSwamdU.clear();
00100 fSwamdV.clear();
00101 fMeasuredU.clear();
00102 fMeasuredV.clear();
00103 fActualZ.clear();
00104 fErrorOnMeasure.clear();
00105 fReverseMapUIndex.clear();
00106 fReverseMapVIndex.clear();
00107 fTrkClstList.Delete();
00108 fDistanceFromStart.clear();
00109 fRangeThisPlane.clear();
00110 fFile->Close();
00111
00112
00113 }
|
|
||||||||||||
|
Definition at line 1751 of file AlgFitTrack3.cxx. References fDirection, and fDistanceFromStart. Referenced by CalculateP(). 01752 {
01753 if(fDirection==1) { // forwards
01754 if(k<i) {
01755 return -1.0*CalculateD(k,i); //Will be negative distance
01756 }
01757 if(i<=fFirstPlane) return -1; // Can't do it that way
01758 if(k>fLastPlane) return -1; // Can't do it that way
01759 if(k==i) return 0;
01760 IntDoubleMap::iterator pos;
01761 Double_t iDistance=0;
01762 pos=fDistanceFromStart.find(i);
01763 if(pos==fDistanceFromStart.end()) {
01764 //Couldn't find that plane
01765 return 0;
01766 }
01767 else iDistance=pos->second;
01768 Double_t kDistance=0;
01769 pos=fDistanceFromStart.find(k);
01770 if(pos==fDistanceFromStart.end()) {
01771 //Couldn't find that plane
01772 return 0;
01773 }
01774 else kDistance=pos->second;
01775 return kDistance-iDistance;
01776 }
01777 else if(fDirection==-1) { //backwards
01778 if(k>i) return -1.0*CalculateD(k,i); //Will be negative distance
01779 if(i>=fHighestPlane) return -1; // Can't do it that way
01780 if(k<fLowestPlane) return -1; // Can't do it that way
01781 if(k==i) return 0;
01782 IntDoubleMap::iterator pos;
01783 Double_t iDistance=0;
01784 pos=fDistanceFromStart.find(i);
01785 if(pos==fDistanceFromStart.end()) {
01786 //Couldn't find that plane
01787 return 0;
01788 }
01789 else iDistance=pos->second;
01790 Double_t kDistance=0;
01791 pos=fDistanceFromStart.find(k);
01792 if(pos==fDistanceFromStart.end()) {
01793 //Couldn't find that plane
01794 return 0;
01795 }
01796 else kDistance=pos->second;
01797 return kDistance-iDistance;
01798
01799 }
01800 return 0;
01801 }
|
|
||||||||||||
|
Definition at line 1601 of file AlgFitTrack3.cxx. References abs(), CalculateD(), fDirection, fFirstPlane, fHighestPlane, fLowestPlane, fTheta0i, fThicknessi, and MSG. Referenced by FillWeightMatrix(). 01602 {
01603
01604 MSG("FitTrack3", Msg::kVerbose)
01605 << "AlgFitTrack3::CalculateP(" << k << "," << n << ")" << endl;
01606 if(k<1 || n<1) return -1;
01607 int highestIndex=abs(fLowestPlane-fHighestPlane);
01608 if(k>highestIndex || n>highestIndex) return -1;
01609
01610 int usek=k;
01611 int usen=n;
01612 if(k>n) {
01613 usek=n;
01614 usen=k;
01615 }
01616 if(fDirection==1) {
01617 if(k==n) {
01618 double val=0;
01619 // int kPlane=fFirstPlane+usek;
01620 int nPlane=fFirstPlane+usen;
01621 double lastTheta=0;
01622 double lastThickness=0;
01623 double thisTheta=0;
01624 double thisThickness=0;
01625 for(int i=1;i<=usen;i++) {
01626 int doingPlane=fFirstPlane+i;
01627
01628 IntDoubleMap::iterator thetaPos;
01629 IntDoubleMap::iterator thicknessPos;
01630 thetaPos = fTheta0i.find(doingPlane);
01631 if(thetaPos==fTheta0i.end()) {
01632 //Missing plane.
01633 if(lastTheta!=0) thisTheta=lastTheta;
01634 else {
01635 int uptoPlane=doingPlane;
01636 while(thetaPos==fTheta0i.end() &&
01637 uptoPlane<fHighestPlane) {
01638 uptoPlane++;
01639 thetaPos=fTheta0i.find(uptoPlane);
01640 }
01641 if(thetaPos==fTheta0i.end()) {
01642 thisTheta=0;
01643 }
01644 else thisTheta=thetaPos->second;
01645 }
01646 }
01647 else thisTheta=thetaPos->second;
01648
01649 thicknessPos = fThicknessi.find(doingPlane);
01650 if(thicknessPos==fThicknessi.end()) {
01651 //Missing plane.
01652 if(lastThickness!=0) thisThickness=lastThickness;
01653 else {
01654 int uptoPlane=doingPlane;
01655 while(thicknessPos==fThicknessi.end() &&
01656 uptoPlane<fHighestPlane) {
01657 uptoPlane++;
01658 thicknessPos=fThicknessi.find(uptoPlane);
01659 }
01660 if(thicknessPos==fThicknessi.end()) {
01661 thisThickness=0;
01662 }
01663 else thisThickness=thicknessPos->second;
01664 }
01665 }
01666 else thisThickness=thicknessPos->second;
01667
01668 val+=(thisTheta*thisTheta)*
01669 ((thisThickness*thisThickness/3.00)+
01670 CalculateD(doingPlane,nPlane)*
01671 (thisThickness+CalculateD(doingPlane,nPlane)));
01672 }
01673
01674 return val;
01675 }
01676 else {
01677 double val=0;
01678 int kPlane=fFirstPlane+usek;
01679 int nPlane=fFirstPlane+usen;
01680 double lastTheta=0;
01681 double lastThickness=0;
01682 double thisTheta=0;
01683 double thisThickness=0;
01684 for(int i=1;i<=usek;i++) {
01685 int doingPlane=fFirstPlane+i;
01686
01687 IntDoubleMap::iterator thetaPos;
01688 IntDoubleMap::iterator thicknessPos;
01689 thetaPos = fTheta0i.find(doingPlane);
01690 if(thetaPos==fTheta0i.end()) {
01691 //Missing plane.
01692 if(lastTheta!=0) thisTheta=lastTheta;
01693 else {
01694 int uptoPlane=doingPlane;
01695 while(thetaPos==fTheta0i.end() &&
01696 uptoPlane<fHighestPlane) {
01697 uptoPlane++;
01698 thetaPos=fTheta0i.find(uptoPlane);
01699 }
01700 if(thetaPos==fTheta0i.end()) {
01701 thisTheta=0;
01702 }
01703 else thisTheta=thetaPos->second;
01704 }
01705 }
01706 else thisTheta=thetaPos->second;
01707
01708 thicknessPos = fThicknessi.find(doingPlane);
01709 if(thicknessPos==fThicknessi.end()) {
01710 //Missing plane.
01711 if(lastThickness!=0) thisThickness=lastThickness;
01712 else {
01713 int uptoPlane=doingPlane;
01714 while(thicknessPos==fThicknessi.end() &&
01715 uptoPlane<fHighestPlane) {
01716 uptoPlane++;
01717 thicknessPos=fThicknessi.find(uptoPlane);
01718 }
01719 if(thicknessPos==fThicknessi.end()) {
01720 thisThickness=0;
01721 }
01722 else thisThickness=thicknessPos->second;
01723 }
01724 }
01725 else thisThickness=thicknessPos->second;
01726
01727 Double_t dik=CalculateD(doingPlane,kPlane);
01728 Double_t din=CalculateD(doingPlane,nPlane);
01729 MSG("FitTrack3", Msg::kVerbose)
01730 << "Summing (" << i << "," << k << "," << n << ")" << endl
01731 << "Theta0i " << thisTheta << endl
01732 << "xi (thickness) " << thisThickness << endl
01733 << "D(i,k) " << dik << endl
01734 << "D(i,n) " << din << endl;
01735
01736 val+=(thisTheta*thisTheta)*
01737 ((thisThickness*thisThickness/3.00)+
01738 (thisThickness/2.0)*(dik+din)+(dik*din));
01739 }
01740
01741 return val;
01742 }
01743
01744 }
01745
01746 return 0;
01747 }
|
|
|
Definition at line 116 of file AlgFitTrack3.cxx. References fActualZ, fDistanceFromStart, fErrorOnMeasure, fIterations, fMeasuredU, fMeasuredV, fRangeThisPlane, fReverseMapUIndex, fReverseMapVIndex, fSwamdU, fSwamdV, fSwamU, fSwamV, fTheta0i, fThicknessi, fTrkClstList, fUWeightMatrix, fVWeightMatrix, and MSG. Referenced by RunAlg(). 00117 {
00118 MSG("FitTrack3", Msg::kVerbose)
00119 << "Starting AlgFitTrack3::CleanUp()" << endl;
00120 fIterations=0;
00121 if(fUWeightMatrix) delete fUWeightMatrix;
00122 fUWeightMatrix=0;
00123 if(fVWeightMatrix) delete fVWeightMatrix;
00124 fVWeightMatrix=0;
00125 fTheta0i.clear(); // mean square scattering angle.
00126 fThicknessi.clear(); // thickness passed through at i.
00127 fSwamU.clear();
00128 fSwamV.clear();
00129 fSwamdU.clear();
00130 fSwamdV.clear();
00131 fMeasuredU.clear();
00132 fMeasuredV.clear();
00133 fActualZ.clear();
00134 fErrorOnMeasure.clear();
00135 fReverseMapUIndex.clear();
00136 fReverseMapVIndex.clear();
00137 fDistanceFromStart.clear();
00138 fRangeThisPlane.clear();
00139 fTrkClstList.Delete();
00140 // MSG("FitTrack3", Msg::kDebug)
00141 // << "Finished AlgFitTrack3::CleanUp()" << endl;
00142 }
|
|
|
Definition at line 1284 of file AlgFitTrack3.cxx. References abs(), fActualZ, fChiSqU, fChiSqV, fDirection, fdUdZ0, fdVdZ0, fHighestPlane, fLowestPlane, fMatrixA, fMatrixC, fMeasuredU, fMeasuredV, fP0, fSwamU, fSwamV, fU0, fV0, GetWiju(), GetWijv(), and MSG. Referenced by RunAlg(). 01284 {
01285
01286 MSG("FitTrack3", Msg::kVerbose)
01287 << "AlgFitTrack3::FillMatricesAandC()" << endl;
01288
01289 fChiSqU=0;
01290 fChiSqV=0;
01291 int sizeOfMatrix=abs(fLowestPlane-fHighestPlane);
01292
01293 Double_t tempMatA[5][5];
01294 Double_t tempMatC[5];
01295 for(int row=0;row<5;row++) {
01296 tempMatC[row]=0;
01297 for(int col=0;col<5;col++) {
01298 tempMatA[row][col]=0;
01299 }
01300 }
01301 for(int i=0;i<sizeOfMatrix;i++) {
01302 for(int j=0;j<sizeOfMatrix;j++) {
01303 int planeI=0;
01304 int planeJ=0;
01305 if(fDirection==1) {
01306 planeI=i+fLowestPlane+1;
01307 planeJ=j+fLowestPlane+1;
01308 }
01309 else {
01310 planeI=fHighestPlane-i-1;
01311 planeJ=fHighestPlane-j-1;
01312 }
01313
01314 Double_t wiju=GetWiju(i,j);
01315 Double_t wijv=GetWijv(i,j);
01316 if(wiju==0 && wijv==0) continue;
01317 MSG("FitTrack3", Msg::kVerbose)
01318 << "Wu(" << i << "," << j << ") = " << wiju << endl
01319 << "Wv(" << i << "," << j << ") = " << wijv << endl;
01320 Double_t zi=0;
01321 Double_t zj=0;
01322 Double_t biu=0;
01323 Double_t biv=0;
01324 Double_t bju=0;
01325 Double_t bjv=0;
01326 Double_t uim=0;
01327 Double_t ujm=0;
01328 Double_t vim=0;
01329 Double_t vjm=0;
01330 Double_t newui=0;
01331 Double_t newvi=0;
01332 Double_t newuj=0;
01333 Double_t newvj=0;
01334
01335
01336
01337 IntDoubleMap::iterator pos;
01338 pos=fActualZ.find(planeI);
01339 if(pos==fActualZ.end()) continue;
01340 else zi=(pos->second-fZ0);
01341 pos=fActualZ.find(planeJ);
01342 if(pos==fActualZ.end()) continue;
01343 else zj=(pos->second-fZ0);
01344
01345 if(wiju!=0) {
01346 pos=fSwamU.find(planeI);
01347 if(pos!=fSwamU.end()) {
01348 biu=(pos->second-fU0-fdUdZ0*(zi))*fP0;//*fCharge;
01349 newui=pos->second;
01350 }
01351 }
01352 if(wijv!=0) {
01353 pos=fSwamV.find(planeI);
01354 if(pos!=fSwamV.end()) {
01355 biv=(pos->second-fV0-fdVdZ0*(zi))*fP0;//*fCharge;
01356 newvi=pos->second;
01357 }
01358 }
01359
01360 if(wiju!=0) {
01361 pos=fSwamU.find(planeJ);
01362 if(pos!=fSwamU.end()) {
01363 bju=(pos->second-fU0-fdUdZ0*(zj))*fP0;//*fCharge;
01364 newuj=pos->second;
01365 }
01366 }
01367 if(wijv!=0) {
01368 pos=fSwamV.find(planeJ);
01369 if(pos!=fSwamV.end()) {
01370 bjv=(pos->second-fV0-fdVdZ0*(zj))*fP0;//*fCharge;
01371 newvj=pos->second;
01372 }
01373 }
01374
01375 if(wiju!=0) {
01376 pos=fMeasuredU.find(planeI);
01377 if(pos!=fMeasuredU.end()) uim=pos->second;
01378 pos=fMeasuredU.find(planeJ);
01379 if(pos!=fMeasuredU.end()) ujm=pos->second;
01380 }
01381 if(wijv!=0) {
01382 pos=fMeasuredV.find(planeI);
01383 if(pos!=fMeasuredV.end()) vim=pos->second;
01384 pos=fMeasuredV.find(planeJ);
01385 if(pos!=fMeasuredV.end()) vjm=pos->second;
01386 }
01387
01388 MSG("FitTrack3", Msg::kVerbose)
01389 << "Wu(" << i << "," << j << ") = " << wiju << endl
01390 << "Wv(" << i << "," << j << ") = " << wijv << endl
01391 << "uim " << uim << endl
01392 << "ujm " << ujm << endl
01393 << "vim " << vim << endl
01394 << "vjm " << vjm << endl
01395 << "ui " << newui << " u0 " << fU0
01396 << " dudz0 " << fdUdZ0 << " (zi-z0) " << zi
01397 << " p0 " << fP0 << endl
01398 << "uj " << newuj << " u0 " << fU0
01399 << " dudz0 " << fdUdZ0 << " (zj-z0) " << zj
01400 << " p0 " << fP0 << endl
01401 << "vi " << newvi << " v0 " << fV0
01402 << " dvdz0 " << fdVdZ0 << " (zi-z0) " << zi
01403 << " p0 " << fP0 << endl
01404 << "vj " << newvj << " v0 " << fV0
01405 << " dvdz0 " << fdVdZ0 << " (zj-z0) " << zj
01406 << " p0 " << fP0 << endl
01407 << "biu(" << i << ") = " << biu << endl
01408 << "bju(" << j << ") = " << bju << endl
01409 << "biv(" << i << ") = " << biv << endl
01410 << "bjv(" << j << ") = " << bjv << endl;
01411
01412 //Calculate Chisq.
01413 fChiSqU+=(newui-uim)*wiju*(newuj-ujm);
01414 fChiSqV+=(newvi-vim)*wijv*(newvj-vjm);
01415
01416 //Now do the silly sums.
01417 // Starting with row 0 of matrix A
01418 tempMatA[0][0]+=2.0*wiju;
01419 tempMatA[0][1]+=wiju*(zi+zj);
01420 tempMatA[0][2]+=wiju*(biu+bju);
01421 //Now row 1
01422 tempMatA[1][0]+=wiju*(zi+zj);
01423 tempMatA[1][1]+=2.0*wiju*(zi*zj);
01424 tempMatA[1][2]+=wiju*(biu*zj+bju*zi);
01425 //Now the long midddle row
01426 tempMatA[2][0]+=wiju*(biu+bju);
01427 tempMatA[2][1]+=wiju*(zi*bju+zj*biu);
01428 tempMatA[2][2]+=2.0*wiju*(biu*bju)+2.0*wijv*(biv*bjv);
01429 tempMatA[2][3]+=wijv*(biv+bjv);
01430 tempMatA[2][4]+=wijv*(biv*zj+bjv*zi);
01431 //Now row 3
01432 tempMatA[3][2]+=wijv*(biv+bjv);
01433 tempMatA[3][3]+=2*wijv;
01434 tempMatA[3][4]+=wijv*(zi+zj);
01435 //Now row 4
01436 tempMatA[4][2]+=wijv*(biv*zj+bjv*zi);
01437 tempMatA[4][3]+=wijv*(zi+zj);
01438 tempMatA[4][4]+=2*wijv*(zi*zj);
01439
01440 //And moving on to matrix C
01441 tempMatC[0]+=wiju*(uim+ujm);
01442 tempMatC[1]+=wiju*(uim*zj+ujm*zi);
01443 tempMatC[2]+=wiju*(uim*bju+ujm*biu)+wijv*(vim*bjv+vjm*biv);
01444 tempMatC[3]+=wijv*(vim+vjm);
01445 tempMatC[4]+=wijv*(vim*zj+vjm*zi);
01446 }
01447 }
01448 // cout << "Matrix C" << endl;
01449 for(int row=0;row<5;row++) {
01450 // cout << row << "\t" << tempMatC[row] << endl;
01451 fMatrixC[row][0]=tempMatC[row];
01452 for(int col=0;col<5;col++) {
01453 fMatrixA[row][col]=tempMatA[row][col];
01454 }
01455 }
01456
01457 }
|
|
|
Definition at line 1460 of file AlgFitTrack3.cxx. References abs(), CalculateP(), fDirection, fErrorOnMeasure, fHighestPlane, fLowestPlane, fMeasuredU, fMeasuredV, fReverseMapUIndex, fReverseMapVIndex, fSizeOfUWeightMatrix, fSizeOfVWeightMatrix, fUWeightMatrix, fVWeightMatrix, and MSG. Referenced by RunAlg(). 01460 {
01461
01462
01463 int sizeOfMatrix=abs(fLowestPlane-fHighestPlane);
01464 MSG("FitTrack3", Msg::kVerbose)
01465 << "AlgFitTrack3::FillWeightMatrix()\tsizeOfMatrix "
01466 << sizeOfMatrix << endl;
01467 TMatrixD invWeightMatrix(sizeOfMatrix,sizeOfMatrix);
01468
01469 int *mapUIndex = new int [fSizeOfUWeightMatrix];
01470 int *mapVIndex = new int [fSizeOfVWeightMatrix];
01471 fReverseMapUIndex.clear();
01472 fReverseMapVIndex.clear();
01473 int uptoUindex=0;
01474 int uptoVindex=0;
01475 IntDoubleMap::iterator pos;
01476 for(int row=0;row<sizeOfMatrix;row++) {
01477 int planeRow=0;
01478 if(fDirection==1) {
01479 planeRow=row+fLowestPlane+1;
01480 }
01481 else {
01482 planeRow=fHighestPlane-row-1;
01483 }
01484
01485 pos=fMeasuredU.find(planeRow);
01486 if(pos!=fMeasuredU.end()) {
01487 //Got a U plane;
01488 mapUIndex[uptoUindex]=row;
01489 fReverseMapUIndex[row]=uptoUindex;
01490 uptoUindex++;
01491 }
01492 pos=fMeasuredV.find(planeRow);
01493 if(pos!=fMeasuredV.end()) {
01494 //Got a U plane;
01495 mapVIndex[uptoVindex]=row;
01496 fReverseMapVIndex[row]=uptoVindex;
01497 uptoVindex++;
01498 }
01499 }
01500
01501 TMatrixD invUWeightMatrix(fSizeOfUWeightMatrix,fSizeOfUWeightMatrix);
01502 TMatrixD invVWeightMatrix(fSizeOfVWeightMatrix,fSizeOfVWeightMatrix);
01503 for(int row=0;row<fSizeOfUWeightMatrix;row++) {
01504
01505 int bigRow=mapUIndex[row];
01506 int planeRow=0;
01507 if(fDirection==1) {
01508 planeRow=bigRow+fLowestPlane+1;
01509 }
01510 else {
01511 planeRow=fHighestPlane-bigRow-1;
01512 }
01513 for(int col=0;col<fSizeOfUWeightMatrix;col++) {
01514 int bigCol=mapUIndex[col];
01515 int planeCol=0;
01516 if(fDirection==1) {
01517 planeCol=bigCol+fLowestPlane+1;
01518 }
01519 else {
01520 planeCol=fHighestPlane-bigCol-1;
01521 }
01522
01523 double epsilon=0;
01524 pos=fErrorOnMeasure.find(planeRow);
01525 if(pos==fErrorOnMeasure.end()) {
01526 invUWeightMatrix[row][col]=0;
01527 continue;
01528 }
01529 else epsilon = pos->second;
01530 if(bigRow!=bigCol) {
01531 epsilon=0;
01532 }
01533 double element = (epsilon*epsilon)+CalculateP(bigRow+1,bigCol+1);
01534 MSG("FitTrack3", Msg::kVerbose)
01535 << "W(" << bigRow+1 << "," << bigCol+1 << ") = "
01536 << element << "\tEpsilion = " << epsilon << endl;
01537 invUWeightMatrix[row][col]=element;
01538 }
01539 }
01540
01541 MSG("FitTrack3", Msg::kVerbose)
01542 << "Filled invUWeightMatrix" << endl;
01543 //Fill invVWeightMatrix
01544 for(int row=0;row<fSizeOfVWeightMatrix;row++) {
01545
01546 int bigRow=mapVIndex[row];
01547 int planeRow=0;
01548 if(fDirection==1) {
01549 planeRow=bigRow+fLowestPlane+1;
01550 }
01551 else {
01552 planeRow=fHighestPlane-bigRow-1;
01553 }
01554 for(int col=0;col<fSizeOfVWeightMatrix;col++) {
01555 int bigCol=mapVIndex[col];
01556 int planeCol=0;
01557 if(fDirection==1) {
01558 planeCol=bigCol+fLowestPlane+1;
01559 }
01560 else {
01561 planeCol=fHighestPlane-bigCol-1;
01562 }
01563
01564 double epsilon=0;
01565 pos=fErrorOnMeasure.find(planeRow);
01566 if(pos==fErrorOnMeasure.end()) {
01567 invVWeightMatrix[row][col]=0;
01568 continue;
01569 }
01570 else epsilon = pos->second;
01571 if(bigRow!=bigCol) {
01572 epsilon=0;
01573 }
01574 double element = (epsilon*epsilon)+CalculateP(bigRow+1,bigCol+1);
01575 MSG("FitTrack3", Msg::kVerbose)
01576 << "W(" << bigRow+1 << "," << bigCol+1 << ") = "
01577 << element << "\tEpsilion = " << epsilon << endl;
01578 invVWeightMatrix[row][col]=element;
01579 }
01580 }
01581
01582 MSG("FitTrack3", Msg::kVerbose)
01583 << "Filled invVWeightMatrix" << endl;
01584
01585 fUWeightMatrix = new TMatrixD(TMatrixD::kInverted,invUWeightMatrix);
01586
01587 MSG("FitTrack3", Msg::kVerbose)
01588 << "Inverted invUWeightMatrix" << endl;
01589 fVWeightMatrix = new TMatrixD(TMatrixD::kInverted,invVWeightMatrix);
01590
01591 MSG("FitTrack3", Msg::kVerbose)
01592 << "Inverted invUWeightMatrix" << endl;
01593 //fWeightMatrix = new TMatrixD(TMatrixD::kInverted,invWeightMatrix);
01594 // cout << "Silly weight matrices: " << endl;
01595 // fUWeightMatrix->Print();
01596 // fVWeightMatrix->Print();
01597 }
|
|
||||||||||||
|
Definition at line 1137 of file AlgFitTrack3.cxx. References fReverseMapUIndex, and MSG. Referenced by FillMatricesAandC(). 01137 {
01138 IntIntMap::iterator pos;
01139 pos=fReverseMapUIndex.find(row);
01140 if(pos==fReverseMapUIndex.end()) return 0;
01141 pos=fReverseMapUIndex.find(col);
01142 if(pos==fReverseMapUIndex.end()) return 0;
01143
01144 int newRow=fReverseMapUIndex[row];
01145 int newCol=fReverseMapUIndex[col];
01146 MSG("FitTrack3", Msg::kVerbose)
01147 << "GetWijv(" << row << "," << col << ")" << endl
01148 << "newRow " << newRow << " newCol " << newCol << endl;
01149
01150 return (*fUWeightMatrix)[newRow][newCol];
01151 }
|
|
||||||||||||
|
Definition at line 1156 of file AlgFitTrack3.cxx. References fReverseMapVIndex, and MSG. Referenced by FillMatricesAandC(). 01156 {
01157 IntIntMap::iterator pos;
01158 pos=fReverseMapVIndex.find(row);
01159 if(pos==fReverseMapVIndex.end()) return 0;
01160 pos=fReverseMapVIndex.find(col);
01161 if(pos==fReverseMapVIndex.end()) return 0;
01162 int newRow=fReverseMapVIndex[row];
01163 int newCol=fReverseMapVIndex[col];
01164 MSG("FitTrack3", Msg::kVerbose)
01165 << "GetWijv(" << row << "," << col << ")" << endl
01166 << "newRow " << newRow << " newCol " << newCol << endl;
01167 return (*fVWeightMatrix)[newRow][newCol];
01168 }
|
|
|
Definition at line 473 of file AlgFitTrack3.cxx. References abs(), TrackClusterSR::AddStrip(), fCharge, fDirection, fdUdZ0, fdVdZ0, fFirstPlane, fHighestPlane, fLastPlane, fLowestPlane, fMyVC, fP0, fParmMisalignmentError, fR_showerU, fR_showerV, fR_totalU, fR_totalV, fR_validU, fR_validV, fReasonsForNotFitting, fReasonsTree, fSizeOfUWeightMatrix, fSizeOfVWeightMatrix, fTrkClstList, fU0, fV0, fZ0, CandRecoHandle::GetBegPlane(), CandHandle::GetDaughterIterator(), VldContext::GetDetector(), CandRecoHandle::GetEndPlane(), CandRecoHandle::GetEndU(), CandRecoHandle::GetEndV(), CandRecoHandle::GetEndZ(), CandTrackHandle::GetMomentum(), TrackClusterSR::GetPlane(), CandStripHandle::GetPlane(), TrackClusterSR::GetPlaneView(), CandStripHandle::GetPlaneView(), VldContext::GetSimFlag(), CandTrackHandle::GetT(), CandRecoHandle::GetTimeSlope(), VldContext::GetTimeStamp(), TrackClusterSR::GetTPos(), TrackClusterSR::GetTPosError(), CandHandle::GetVldContext(), UgliGeomHandle::GetZExtent(), TrackClusterSR::GetZPos(), CandTrackHandle::IsInShower(), min(), MSG, and TrackClusterSR::SetTime3D(). Referenced by RunAlg(). 00474 {
00475 MSG("FitTrack3", Msg::kVerbose)
00476 << "AlgFitTrack3::InitializeTrkClstList()" << endl;
00477 // Int_t dplane = abs(track0->GetEndPlane()-track0->GetBegPlane());
00478 Int_t begplane = track0->GetBegPlane();
00479 Int_t endplane = track0->GetEndPlane();
00480 fFirstPlane=begplane;
00481 fLastPlane=endplane;
00482 const VldContext *myvc = track0->GetVldContext();
00483 fMyVC = new VldContext(myvc->GetDetector(),myvc->GetSimFlag(),
00484 myvc->GetTimeStamp());
00485 UgliGeomHandle ugh(*myvc);
00486 if(endplane<begplane) {
00487 Int_t temp=begplane;
00488 begplane=endplane;
00489 endplane=temp;
00490 }
00491 fLowestPlane=begplane;
00492 fHighestPlane=endplane;
00493 CandStripHandleItr stripItr(track0->GetDaughterIterator());
00494 CandStripHandleKeyFunc *stripKf = stripItr.CreateKeyFunc();
00495 stripKf->SetFun(CandStripHandle::KeyFromPlane);
00496 stripItr.GetSet()->AdoptSortKeyFunc(stripKf);
00497 stripKf = 0;
00498 Int_t izfor=1;
00499 if (track0->GetTimeSlope()>0.) {
00500 stripItr.SetDirection(1);
00501 izfor=1;
00502 }
00503 else {
00504 stripItr.SetDirection(-1);
00505 izfor=-1;
00506 }
00507 fDirection=izfor;
00508 Int_t oldplane=0;
00509 Int_t oldplane2=0;
00510 Int_t nplane=0;
00511 Int_t nplaneu=0;
00512 Int_t nplanev=0;
00513 Int_t isFirstPlaneU=0;
00514 stripItr.Reset();
00515 TrackClusterSR *oldtc=0;
00516
00517 fLowestPlane=1000;
00518 fHighestPlane=-1000;
00519
00520
00521 int numShowerPlanes=0;
00522 int numShowerStrips=0;
00523 map<int,int> showerPlanes;
00524 map<int,int>::iterator showerPos;
00525 TObjArray tempTrkClstList;
00526
00527 fR_totalU=0;
00528 fR_totalV=0;
00529 fR_showerU=0;
00530 fR_showerV=0;
00531 fR_validU=0;
00532 fR_validV=0;
00533 int lastTotal=-1;
00534 // int lastValid=-1;
00535 int lastShower=-1;
00536
00537 while (CandStripHandle *strip = stripItr()) {
00538 // cfh.SetInShower(strip,track0->IsInShower(strip));
00539 //Just counts the number of planes hit in each view.
00540 if(lastTotal!=strip->GetPlane()) {
00541 switch (strip->GetPlaneView()) {
00542 case PlaneView::kU:
00543 fR_totalU++;
00544 break;
00545 case PlaneView::kV:
00546 fR_totalV++;
00547 break;
00548 default:
00549 break;
00550 }
00551 }
00552 lastTotal=strip->GetPlane();
00553 if(!(track0->IsInShower(strip)>1)) {
00554 if (oldplane!=strip->GetPlane()) {
00555 if (strip->GetPlane()>=begplane) {
00556 if(strip->GetPlane()<fLowestPlane)
00557 fLowestPlane=strip->GetPlane();
00558 if(strip->GetPlane()>fHighestPlane)
00559 fHighestPlane=strip->GetPlane();
00560 nplane++;
00561 switch (strip->GetPlaneView()) {
00562 case PlaneView::kU:
00563 if(isFirstPlaneU==0) isFirstPlaneU=1;
00564 nplaneu++;
00565 fR_validU++;
00566 break;
00567 case PlaneView::kV:
00568 if(isFirstPlaneU==0) isFirstPlaneU=-1;
00569 nplanev++;
00570 fR_validV++;
00571 break;
00572 default:
00573 break;
00574 }
00575 }
00576 oldplane = strip->GetPlane();
00577 }
00578
00579
00580 if (strip->GetPlane()>=begplane) {
00581 // cfh.AddDaughterLink(*strip);
00582
00583 if (!oldtc) {
00584 oldtc = new TrackClusterSR(strip, fParmMisalignmentError);
00585 oldtc->SetTime3D(track0->GetT(oldtc->GetPlane()));
00586 oldplane2 = strip->GetPlane();
00587 }
00588 else if (strip->GetPlane()==oldplane2) {
00589 oldtc->AddStrip(strip);
00590 }
00591 else {
00592 TrackClusterSR *newtc = new TrackClusterSR(*oldtc);
00593 tempTrkClstList.Add(newtc);
00594 delete oldtc;
00595 oldplane2 = strip->GetPlane();
00596 oldtc = new TrackClusterSR(strip, fParmMisalignmentError);
00597 oldtc->SetTime3D(track0->GetT(oldtc->GetPlane()));
00598 }
00599 }
00600 }
00601 else {
00602 // MSG("FitTrack3", Msg::kDebug)
00603 // << "Got shower strip: plane "
00604 // << strip->GetPlane() << " strip " << strip->GetStrip()
00605 // << " InShower " << track0->IsInShower(strip) << endl;
00606
00607 if(lastShower!=strip->GetPlane()) {
00608 switch (strip->GetPlaneView()) {
00609 case PlaneView::kU:
00610 fR_showerU++;
00611 break;
00612 case PlaneView::kV:
00613 fR_showerV++;
00614 break;
00615 default:
00616 break;
00617 }
00618 }
00619 lastShower=strip->GetPlane();
00620 showerPos=showerPlanes.find(strip->GetPlane());
00621 if(showerPos!=showerPlanes.end()) {
00622 showerPlanes[strip->GetPlane()]++;
00623 numShowerStrips+=track0->IsInShower(strip);
00624 }
00625 else {
00626 numShowerPlanes++;
00627 numShowerStrips+=track0->IsInShower(strip);
00628 showerPlanes[strip->GetPlane()]=1;
00629 }
00630 }
00631
00632 }
00633
00634 if (oldtc) {
00635 TrackClusterSR *newtc = new TrackClusterSR(*oldtc);
00636 tempTrkClstList.Add(newtc);
00637 delete oldtc;
00638 }
00639
00640
00641 MSG("FitTrack3", Msg::kVerbose)
00642 << "Lowest Plane " << fLowestPlane
00643 << " Highest Plane " << fHighestPlane
00644 << " Direction " << izfor << endl;
00645
00646 //At this point each plane is a TrackClusterSR object in fTrkClstList
00647 MSG("FitTrack3", Msg::kVerbose)
00648 << "Filled fTrkClstList:\tnplane " << nplane
00649 << "\tnplaneu " << nplaneu << "\tnplanev" << nplanev << endl
00650 << "fTrkClstList.GetEntries() " << fTrkClstList.GetEntries() << endl;
00651
00652 if(nplaneu<8 || nplanev<8) {
00653 MSG("FitTrack3", Msg::kWarning)
00654 << "Too few non-shower planes to fit"
00655 << endl;
00656 fReasonsForNotFitting->Fill(1); //Too few planes total.
00657 fReasonsTree->Fill();
00658 fReasonsTree->Write("reasonsTree",TObject::kOverwrite);
00659 return -1;
00660 }
00661
00662
00663 Int_t startPlane[10];
00664 Int_t endPlane[10];
00665 Int_t numPlanesInSection[10];
00666 Int_t lastPlanenum=-1000;
00667 Int_t numSections=0;
00668 for (Int_t i=0; i<=tempTrkClstList.GetLast() ; i++) {
00669 TrackClusterSR *thisPlane =
00670 dynamic_cast<TrackClusterSR*>(tempTrkClstList.At(i));
00671 Int_t planenum=thisPlane->GetPlane();
00672 if(abs(planenum-lastPlanenum)>5) { //Missing 5 planes
00673 numSections++;
00674 startPlane[numSections-1]=planenum;
00675 endPlane[numSections-1]=planenum;
00676 numPlanesInSection[numSections-1]=1;
00677 }
00678 else {
00679 endPlane[numSections-1]=planenum;
00680 numPlanesInSection[numSections-1]++;
00681 }
00682
00683 lastPlanenum=planenum;
00684 }
00685 int biggestSection=0;
00686 int numInBiggestSection=0;
00687 for(int i=0; i<numSections; i++) {
00688 if(numPlanesInSection[i]>numInBiggestSection) {
00689 numInBiggestSection=numPlanesInSection[i];
00690 biggestSection=i;
00691 }
00692 }
00693
00694 int started=0;
00695 int ended=0;
00696 isFirstPlaneU=0;
00697 nplane=0;
00698 nplaneu=0;
00699 nplanev=0;
00700 for (Int_t i=0; i<=tempTrkClstList.GetLast() ; i++) {
00701 TrackClusterSR *thisPlane =
00702 dynamic_cast<TrackClusterSR*>(tempTrkClstList.At(i));
00703 Int_t planenum=thisPlane->GetPlane();
00704 if(planenum==startPlane[biggestSection]) {
00705 started=1;
00706 }
00707 if(started==1 && ended==0) {
00708 TrackClusterSR *newtc = new TrackClusterSR(*thisPlane);
00709 fTrkClstList.Add(newtc);
00710 nplane++;
00711 switch (thisPlane->GetPlaneView()) {
00712 case PlaneView::kU:
00713 if(isFirstPlaneU==0) isFirstPlaneU=1;
00714 nplaneu++;
00715 break;
00716 case PlaneView::kV:
00717 if(isFirstPlaneU==0) isFirstPlaneU=-1;
00718 nplanev++;
00719 break;
00720 default:
00721 break;
00722 }
00723 }
00724 if(planenum==endPlane[biggestSection]) ended=1;
00725 }
00726 tempTrkClstList.Delete();
00727 fLowestPlane=startPlane[biggestSection];
00728 fHighestPlane=endPlane[biggestSection];
00729 if(fLowestPlane>fHighestPlane) {
00730 int temp=fHighestPlane;
00731 fHighestPlane=fLowestPlane;
00732 fLowestPlane=temp;
00733 }
00734
00735
00736 fSizeOfVWeightMatrix=nplaneu;
00737 fSizeOfUWeightMatrix=nplanev;
00738 if(isFirstPlaneU==1) fSizeOfVWeightMatrix--;
00739 else fSizeOfUWeightMatrix--;
00740
00741 if(nplaneu<8 || nplanev<8) {
00742 MSG("FitTrack3", Msg::kWarning)
00743 << "Too few non-shower planes to fit in section"
00744 << endl;
00745 fReasonsForNotFitting->Fill(2); //Too few planes section.
00746 return -1;
00747 }
00748 numShowerStrips=0;
00749 numShowerPlanes=0;
00750 for(Int_t planenum=fLowestPlane;planenum<=fHighestPlane;planenum++) {
00751 showerPos=showerPlanes.find(planenum);
00752 if(showerPos!=showerPlanes.end()) {
00753 numShowerPlanes++;
00754 numShowerStrips+=showerPos->second;
00755 }
00756 }
00757
00758
00759 if(numShowerPlanes>7 && numShowerStrips>10) {
00760 MSG("FitTrack3", Msg::kWarning)
00761 << "Too many shower planes to fit: planes "
00762 << numShowerPlanes << " strips " << numShowerStrips
00763 << endl;
00764 fReasonsForNotFitting->Fill(3); //Too many shower planes.
00765 return -1;
00766 }
00767 Double_t *zUPos = new Double_t [nplanev];
00768 Double_t *zVPos = new Double_t [nplaneu];
00769 Double_t *uPos = new Double_t [nplanev];
00770 Double_t *vPos = new Double_t [nplaneu];
00771 Double_t *zUPosErr = new Double_t [nplanev];
00772 Double_t *zVPosErr = new Double_t [nplaneu];
00773 Double_t *uPosErr = new Double_t [nplanev];
00774 Double_t *vPosErr = new Double_t [nplaneu];
00775 Int_t uptoU=0;
00776 Int_t uptoV=0;
00777 for (Int_t i=0; i<=fTrkClstList.GetLast() ; i++) {
00778 TrackClusterSR *thisPlane =
00779 dynamic_cast<TrackClusterSR*>(fTrkClstList.At(i));
00780 Double_t tempZPos=thisPlane->GetZPos();
00781 Double_t tempTPos=thisPlane->GetTPos();
00782 Double_t tempZPosErr=0.005;
00783 Double_t tempTPosErr=thisPlane->GetTPosError();
00784
00785 switch (thisPlane->GetPlaneView()) {
00786 case PlaneView::kU:
00787 zVPos[uptoV]=tempZPos;
00788 vPos[uptoV]=tempTPos;
00789 zVPosErr[uptoV]=tempZPosErr;
00790 vPosErr[uptoV]=tempTPosErr;
00791 uptoV++;
00792 break;
00793 case PlaneView::kV:
00794 zUPos[uptoU]=tempZPos;
00795 uPos[uptoU]=tempTPos;
00796 zUPosErr[uptoU]=tempZPosErr;
00797 uPosErr[uptoU]=tempTPosErr;
00798 uptoU++;
00799 break;
00800 default:
00801 break;
00802 }
00803 }
00804 Double_t uFigureOfStraightness=0;
00805 Double_t vFigureOfStraightness=0;
00806 {
00807 TGraphErrors grU(nplanev,zUPos,uPos,zUPosErr,uPosErr);
00808 TGraphErrors grV(nplaneu,zVPos,vPos,zVPosErr,vPosErr);
00809 TF1 fitty("fitty","pol1");
00810 grU.Fit("fitty","Q");
00811 Double_t uSlope=fitty.GetParameter(0);
00812 Double_t uSlopeErr=fitty.GetParError(0);
00813 Double_t uIntercept=fitty.GetParameter(1);
00814 Double_t uInterceptErr=fitty.GetParError(1);
00815 Double_t uChiSq=fitty.GetChisquare();
00816 Double_t uNDF=fitty.GetNDF();
00817 uFigureOfStraightness=uChiSq/uNDF;
00818 grV.Fit("fitty","Q");
00819 Double_t vSlope=fitty.GetParameter(0);
00820 Double_t vSlopeErr=fitty.GetParError(0);
00821 Double_t vIntercept=fitty.GetParameter(1);
00822 Double_t vInterceptErr=fitty.GetParError(1);
00823 Double_t vChiSq=fitty.GetChisquare();
00824 Double_t vNDF=fitty.GetNDF();
00825 vFigureOfStraightness=vChiSq/vNDF;
00826 MSG("FitTrack3", Msg::kVerbose)
00827 << endl
00828 << "uSlope: " << uSlope << endl
00829 << "uSlopeErr: " << uSlopeErr << endl
00830 << "uIntercept: " << uIntercept << endl
00831 << "uInterceptErr: " << uInterceptErr << endl
00832 << "uChiSq: " << uChiSq << endl
00833 << "uNDF: " << uNDF << endl
00834 << "uFigureOfStraightness: " << uFigureOfStraightness << endl
00835 << "vSlope: " << vSlope << endl
00836 << "vSlopeErr: " << vSlopeErr << endl
00837 << "vIntercept: " << vIntercept << endl
00838 << "vInterceptErr: " << vInterceptErr << endl
00839 << "vChiSq: " << vChiSq << endl
00840 << "vNDF: " << vNDF << endl
00841 << "vFigureOfStraightness: " << vFigureOfStraightness << endl;
00842 }
00843 delete [] zUPos;
00844 delete [] zVPos;
00845 delete [] uPos;
00846 delete [] vPos;
00847 delete [] zUPosErr;
00848 delete [] zVPosErr;
00849 delete [] uPosErr;
00850 delete [] vPosErr;
00851
00852 if(uFigureOfStraightness<0.4 || vFigureOfStraightness<0.4) {
00853 MSG("FitTrack3", Msg::kWarning)
00854 << "Too straight to bother fitting: uChiSq/uNDF "
00855 << uFigureOfStraightness << " vChiSq/vNDF "
00856 << vFigureOfStraightness << endl;
00857 fReasonsForNotFitting->Fill(4); //Too straight
00858 return -1;
00859 }
00860
00861 TrackClusterSR *firstPlane =
00862 dynamic_cast<TrackClusterSR*>(fTrkClstList.At(0));
00863 Double_t startU=0;
00864 Double_t startV=0;
00865 Double_t startZ=firstPlane->GetZPos();
00866 Double_t startdUdZ=0;
00867 Double_t startdVdZ=0;
00868 Double_t startMomentum=track0->GetMomentum()*1.1;
00869
00870
00871 Float_t zextent[2];
00872 ugh.GetZExtent(zextent[0],zextent[1]);
00873 Double_t endz=track0->GetEndZ();
00874 Double_t endu=track0->GetEndU();
00875 Double_t endv=track0->GetEndV();
00876 Double_t endx= .70710678*(endu-endv);
00877 Double_t endy= .70710678*(endu+endv);
00878
00879 Double_t d2endz=min(endz-zextent[0],zextent[1]-endz);
00880 Double_t d2endv=0;
00881 Double_t d2endu=0;
00882 Double_t d2endx=0;
00883 Double_t d2endy=0;
00884 switch (fMyVC->GetDetector()) {
00885 case DetectorType::kNear:
00886
00887 break;
00888 case DetectorType::kFar:
00889 // assume 8 m octagon
00890 d2endu = min(4.-endu,4.+endu);
00891 d2endv = min(4.-endv,4.+endv);
00892 d2endx = min(4.-endx,4.+endx);
00893 d2endy = min(4.-endy,4.+endy);
00894 break;
00895 case DetectorType::kCalDet:
00896 // assume 1 m rectangle
00897 break;
00898 default:
00899 MSG("EventSR",Msg::kWarning) << "Detector type " << fMyVC->GetDetector()
00900 << " not supported.\n";
00901 break;
00902 }
00903 Double_t d2endmin=min(min(min(d2endu,d2endv),min(d2endx,d2endy)),d2endz);
00904 MSG("FitTrack3", Msg::kDebug)
00905 << endl
00906 << "d2endmin: " << d2endmin << endl
00907 << "d2endz: " << d2endz << endl
00908 << "d2endu: " << d2endu << endl
00909 << "d2endv: " << d2endv << endl
00910 << "d2endx: " << d2endx << endl
00911 << "d2endy: " << d2endy << endl;
00912
00913 if(d2endmin<0.1) startMomentum+=4; //Not contained
00914
00915 fZ0=startZ;
00916 fU0=startU;
00917 fV0=startV;
00918 fdUdZ0=startdUdZ;
00919 fdVdZ0=startdVdZ;
00920 fP0=startMomentum;
00921 switch (firstPlane->GetPlaneView()) {
00922 case PlaneView::kU:
00923 startV=firstPlane->GetTPos();
00924 break;
00925 case PlaneView::kV:
00926 startU=firstPlane->GetTPos();
00927 break;
00928 default:
00929 break;
00930 }
00931
00932 MSG("FitTrack3", Msg::kVerbose)
00933 << "Start Momentum " << fP0
00934 << " Start U " << startU
00935 << " Start V " << startV << endl;
00936 int doneSame=0;
00937 int doneOther=0;
00938 double firstOtherZ=0;
00939 double firstOtherTPos=0;
00940 int gotFirstOther=0;
00941 for (Int_t i=1; i<=fTrkClstList.GetLast() ; i++) {
00942 TrackClusterSR *thisPlane =
00943 dynamic_cast<TrackClusterSR*>(fTrkClstList.At(i));
00944 if(thisPlane->GetPlaneView()==firstPlane->GetPlaneView()
00945 &&!doneSame) {
00946 Double_t newZ=thisPlane->GetZPos();
00947 Double_t newTPos=thisPlane->GetTPos();
00948 Double_t startTPos=firstPlane->GetTPos();
00949 if(newZ==startZ) continue; //should never happen
00950 Double_t dTdZ=(newTPos-startTPos)/(newZ-startZ);
00951 switch (firstPlane->GetPlaneView()) {
00952 case PlaneView::kU:
00953 startdVdZ=dTdZ;
00954 break;
00955 case PlaneView::kV:
00956 startdUdZ=dTdZ;
00957 break;
00958 default:
00959 break;
00960 }
00961 doneSame=1;
00962 }
00963 else if(!doneOther) {
00964 if(!gotFirstOther) {
00965 firstOtherZ=thisPlane->GetZPos();
00966 firstOtherTPos=thisPlane->GetTPos();
00967 gotFirstOther=1;
00968 }
00969 else {
00970 Double_t newZ=thisPlane->GetZPos();
00971 Double_t newTPos=thisPlane->GetTPos();
00972 if(newZ==firstOtherZ) continue; //should never happen
00973 Double_t dTdZ=(newTPos-firstOtherTPos)/(newZ-firstOtherZ);
00974 Double_t startT=firstOtherTPos+(startZ-firstOtherZ)*dTdZ;
00975 switch (firstPlane->GetPlaneView()) {
00976 case PlaneView::kU:
00977 startdUdZ=dTdZ;
00978 startU=startT;
00979 break;
00980 case PlaneView::kV:
00981 startdVdZ=dTdZ;
00982 startV=startT;
00983 break;
00984 default:
00985 break;
00986 }
00987 doneOther=1;
00988 }
00989
00990
00991
00992 }
00993 if(doneSame && doneOther) break;
00994
00995 }
00996 Int_t charge=-1; //Probably is.
00997 fCharge=charge;
00998 fZ0=startZ;
00999 fU0=startU;
01000 fV0=startV;
01001 fdUdZ0=startdUdZ;
01002 fdVdZ0=startdVdZ;
01003 fP0=startMomentum;
01004 MSG("FitTrack3", Msg::kVerbose)
01005 << "Start Momentum " << fP0 << endl
01006 << " Start U " << fU0 << endl
01007 << " Start V " << fV0 << endl
01008 << " Start dUdZ " << fdUdZ0 << endl
01009 << " Start dVdZ " << fdVdZ0 << endl;
01010 // Have now initialised the following variables.
01011 //
01012 // Double_t startU;
01013 // Double_t startV;
01014 // Double_t startZ;
01015 // Double_t startdUdZ;
01016 // Double_t startdVdZ;
01017 // Int_t startPlane;
01018 // Double_t startMomentum;
01019 // Int_t izfor;
01020
01021 return 1;
01022 }
|
|
|
Definition at line 1174 of file AlgFitTrack3.cxx. References fCharge, fDdUdZ0, fDdVdZ0, fDP0, fDU0, fdUdZ0, fDV0, fdVdZ0, fIsNaturalP0, fLastCharge, fLastdUdZ0, fLastdVdZ0, fLastP0, fLastU0, fLastV0, fLastZ0, fMatrixA, fMatrixC, fMomFromRange, fP0, fU0, fV0, fZ0, and MSG. Referenced by RunAlg(). 01174 {
01175
01176 MSG("FitTrack3", Msg::kVerbose)
01177 << "AlgFitTrack3::InvertMatrixAndGetParams()" << endl;
01178
01179 fLastCharge=fCharge;
01180 fLastZ0=fZ0; //Doesn't change
01181 fLastU0=fU0;
01182 fLastV0=fV0;
01183 fLastdUdZ0=fdUdZ0;
01184 fLastdVdZ0=fdVdZ0;
01185 fLastP0=fP0;
01186
01187 // cout << endl << endl << "Matrix A =" << endl;
01188 // fMatrixA.Print();
01189 // cout << endl << endl << "Matrix C =" << endl;
01190 // fMatrixC.Print();
01191 if(fMatrixA.Determinant()!=0) {
01192 TMatrixD invertedMatrixA(TMatrixD::kInverted,fMatrixA);
01193 TMatrixD solution(invertedMatrixA,TMatrixD::kMult,fMatrixC);
01194 // cout << endl << endl << "Inverted Matrix A =" << endl;
01195 // invertedMatrixA.Print();
01196 // cout << endl << endl << "Matrix B =" << endl;
01197 // solution.Print();
01198 fU0=solution[0][0];
01199 fdUdZ0=solution[1][0];
01200 if(solution[2][0]!=0) fP0=1.0/solution[2][0];
01201 fV0=solution[3][0];
01202 fdVdZ0=solution[4][0];
01203 if(fP0<0) {
01204 fP0=-fP0;
01205 fCharge=-fCharge;
01206 }
01207
01208 if(invertedMatrixA[0][0]>0) {
01209 fDU0=TMath::Sqrt(invertedMatrixA[0][0]);
01210 }
01211 else {
01212 fDU0=0;
01213 }
01214 if(invertedMatrixA[1][1]>0) {
01215 fDdUdZ0=TMath::Sqrt(invertedMatrixA[1][1]);
01216 }
01217 else {
01218 fDdUdZ0=0;
01219 }
01220 Double_t fDOneOverP0;
01221 if(invertedMatrixA[2][2]>0) {
01222 fDOneOverP0=TMath::Sqrt(invertedMatrixA[2][2]);
01223 }
01224 else {
01225 fDOneOverP0=0;
01226 }
01227 if(solution[2][0]!=0) fDP0=fDOneOverP0/solution[2][0];
01228 if(fDP0<0) fDP0=-fDP0;
01229 if(invertedMatrixA[3][3]>0) {
01230 fDV0=TMath::Sqrt(invertedMatrixA[3][3]);
01231 }
01232 else {
01233 fDV0=0;
01234 }
01235 if(invertedMatrixA[4][4]>0) {
01236 fDdVdZ0=TMath::Sqrt(invertedMatrixA[4][4]);
01237 }
01238 else {
01239 fDdVdZ0=0;
01240 }
01241 fIsNaturalP0=1;
01242 if(fP0<0.8*fMomFromRange) {
01243 fIsNaturalP0=0;
01244 MSG("FitTrack3",Msg::kError)
01245 << "Momentum has gone too low. Now: " << fP0 << " from range: "
01246 << fMomFromRange << " will try: " << fMomFromRange*0.9 << endl;
01247 fP0=fMomFromRange*0.9;
01248 }
01249
01250 if(fDdVdZ0<=0 || fDV0<=0 || fDOneOverP0<=0 || fDdUdZ0<=0 || fDU0<=0) {
01251 MSG("FitTrack3",Msg::kError)
01252 << "Negative errors. It's all gone horribly wrong"
01253 << endl;
01254 fIsNaturalP0=0;
01255 }
01256
01257 MSG("FitTrack3", Msg::kDebug)
01258 << endl
01259 << "Charge: " << fCharge << " was " << fLastCharge << endl
01260 << "Z0: " << fZ0 << " was " << fLastZ0 << endl
01261 << "U0: " << fU0 << " +/- " << fDU0 << " was " << fLastU0 << endl
01262 << "V0: " << fV0 << " +/- " << fDV0 << " was " << fLastV0 << endl
01263 << "dUdZ0: " << fdUdZ0 << " +/- " << fDdUdZ0
01264 << " was " << fLastdUdZ0 << endl
01265 << "dVdZ0: " << fdVdZ0 << " +/- " << fDdVdZ0
01266 << " was " << fLastdVdZ0 << endl
01267 << "P0: " << fP0 << " +/- " << fDP0 << " was " << fLastP0 << endl;
01268
01269
01270 }
01271 else {
01272 MSG("FitTrack3",Msg::kError)
01273 << "Can't invert MatrixA"
01274 << endl;
01275 fP0+=1;
01276 fIsNaturalP0=0;
01277 }
01278
01279
01280 }
|
|
||||||||||||||||
|
Implements AlgBase. Definition at line 144 of file AlgFitTrack3.cxx. References CandHandle::AddDaughterLink(), CleanUp(), fCharge, fChiSqU, fChiSqV, fDdUdZ0, fDdVdZ0, fDirection, fDP0, fDU0, fdUdZ0, fDV0, fdVdZ0, fFile, FillMatricesAandC(), FillWeightMatrix(), fIsNaturalP0, fLastP0, fMomFromRange, fP0, fParmMaxIterate, fParmMisalignmentError, fRangeThisPlane, fReasonsForNotFitting, fSwamdU, fSwamdV, fSwamU, fSwamV, fTrackNum, fTrkClstList, fU0, fV0, CandRecoHandle::GetCandSlice(), CandFitTrackHandle::GetChi2(), CandContext::GetDataIn(), CandHandle::GetDaughterIterator(), CandRecoHandle::GetDirCosU(), CandRecoHandle::GetDirCosV(), CandRecoHandle::GetDirCosZ(), CandTrackHandle::GetdS(), CandFitTrackHandle::GetEMCharge(), CandRecoHandle::GetEndDirCosU(), CandRecoHandle::GetEndDirCosV(), CandRecoHandle::GetEndDirCosZ(), CandRecoHandle::GetEndPlane(), CandRecoHandle::GetEndT(), CandRecoHandle::GetEndU(), CandRecoHandle::GetEndV(), CandRecoHandle::GetEndZ(), Registry::GetInt(), CandTrackHandle::GetMomentum(), CandFitTrackHandle::GetMomentumCurve(), CandFitTrackHandle::GetMomentumRange(), CandFitTrackHandle::GetPass(), TrackClusterSR::GetPlane(), CandTrackHandle::GetRange(), TrackClusterSR::GetStripList(), CandRecoHandle::GetTimeOffset(), CandRecoHandle::GetTimeSlope(), CandTrackHandle::GetU(), CandTrackHandle::GetV(), CandRecoHandle::GetVtxPlane(), CandRecoHandle::GetVtxT(), CandRecoHandle::GetVtxU(), CandRecoHandle::GetVtxV(), CandRecoHandle::GetVtxZ(), InitializeTrkClstList(), InvertMatrixAndGetParams(), CandTrackHandle::IsInShower(), CandTrackHandle::IsTPosValid(), MSG, CandRecoHandle::SetCandSlice(), CandFitTrackHandle::SetChi2(), CandRecoHandle::SetDirCosU(), CandRecoHandle::SetDirCosV(), CandRecoHandle::SetDirCosZ(), CandTrackHandle::SetdS(), CandFitTrack3Handle::SetdUdZ(), CandFitTrack3Handle::SetdUdZ0(), CandFitTrack3Handle::SetdUdZ0Err(), CandFitTrack3Handle::SetdUdZ0Initial(), CandFitTrack3Handle::SetdVdZ(), CandFitTrack3Handle::SetdVdZ0(), CandFitTrack3Handle::SetdVdZ0Err(), CandFitTrack3Handle::SetdVdZ0Initial(), CandFitTrackHandle::SetEMCharge(), CandRecoHandle::SetEndDirCosU(), CandRecoHandle::SetEndDirCosV(), CandRecoHandle::SetEndDirCosZ(), CandRecoHandle::SetEndPlane(), CandRecoHandle::SetEndT(), CandRecoHandle::SetEndU(), CandRecoHandle::SetEndV(), CandRecoHandle::SetEndZ(), CandFitTrack3Handle::SetInitialQP(), CandTrackHandle::SetInShower(), CandFitTrackHandle::SetMomentumCurve(), CandFitTrackHandle::SetMomentumRange(), CandFitTrack3Handle::SetNIterate(), CandFitTrack3Handle::SetP0(), CandFitTrack3Handle::SetP0Err(), CandFitTrack3Handle::SetP0Initial(), CandFitTrackHandle::SetPass(), CandTrackHandle::SetRange(), CandRecoHandle::SetTimeOffset(), CandRecoHandle::SetTimeSlope(), CandTrackHandle::SetU(), CandFitTrack3Handle::SetU0(), CandFitTrack3Handle::SetU0Err(), CandFitTrack3Handle::SetU0Initial(), CandTrackHandle::SetV(), CandFitTrack3Handle::SetV0(), CandFitTrack3Handle::SetV0Err(), CandFitTrack3Handle::SetV0Initial(), CandRecoHandle::SetVtxPlane(), CandRecoHandle::SetVtxT(), CandRecoHandle::SetVtxU(), CandRecoHandle::SetVtxV(), CandRecoHandle::SetVtxZ(), SwimAndFillMaps(), and tc. 00145 {
00146 MSG("FitTrack3", Msg::kVerbose) << "Starting AlgFitTrack3::RunAlg()" << endl;
00147
00148 fParmMisalignmentError = ac.GetInt("MisalignmentError")*Munits::mm;
00149
00150 fParmMaxIterate = ac.GetInt("MaxIterate");
00151
00152 assert(ch.InheritsFrom("CandFitTrack3Handle"));
00153 CandFitTrack3Handle &cfh = dynamic_cast<CandFitTrack3Handle &>(ch);
00154
00155 cfh.SetPass(0);
00156
00157 assert(cx.GetDataIn());
00158
00159 assert(cx.GetDataIn()->InheritsFrom("TObjArray"));
00160
00161
00162 const CandTrackHandle *track0 = 0;
00163 const TObjArray *cxin = dynamic_cast<const TObjArray *>(cx.GetDataIn());
00164 for (Int_t i=0; i<=cxin->GetLast(); i++) {
00165 TObject *tobj = cxin->At(i);
00166 if (tobj->InheritsFrom("CandTrackHandle")) {
00167 track0 = dynamic_cast<CandTrackHandle*>(tobj);
00168 }
00169 }
00170
00171 assert(track0);
00172
00173 fTrackNum++;
00174 //Set some standard CandTrack/CandFitTrack thingies.
00175 cfh.SetCandSlice(track0->GetCandSlice());
00176 cfh.SetDirCosU(track0->GetDirCosU());
00177 cfh.SetDirCosV(track0->GetDirCosV());
00178 cfh.SetDirCosZ(track0->GetDirCosZ());
00179 cfh.SetVtxU(track0->GetVtxU());
00180 cfh.SetVtxV(track0->GetVtxV());
00181 cfh.SetVtxZ(track0->GetVtxZ());
00182 cfh.SetVtxT(track0->GetVtxT());
00183 cfh.SetVtxPlane(track0->GetVtxPlane());
00184 cfh.SetEndDirCosU(track0->GetEndDirCosU());
00185 cfh.SetEndDirCosV(track0->GetEndDirCosV());
00186 cfh.SetEndDirCosZ(track0->GetEndDirCosZ());
00187 cfh.SetEndU(track0->GetEndU());
00188 cfh.SetEndV(track0->GetEndV());
00189 cfh.SetEndZ(track0->GetEndZ());
00190 cfh.SetEndT(track0->GetEndT());
00191 cfh.SetEndPlane(track0->GetEndPlane());
00192 cfh.SetTimeSlope(track0->GetTimeSlope());
00193 cfh.SetTimeOffset(track0->GetTimeOffset());
00194 cfh.SetMomentumRange(track0->GetMomentum());
00195 fMomFromRange=track0->GetMomentum();
00196
00197
00198 if(InitializeTrkClstList(track0)==-1) {
00199 //Can't go on
00200
00201 CandStripHandleItr stripItr(track0->GetDaughterIterator());
00202 CandStripHandleKeyFunc *stripKf = stripItr.CreateKeyFunc();
00203 stripKf->SetFun(CandStripHandle::KeyFromPlane);
00204 stripItr.GetSet()->AdoptSortKeyFunc(stripKf);
00205 stripKf = 0;
00206 Int_t izfor=1;
00207 if (track0->GetTimeSlope()>0.) {
00208 stripItr.SetDirection(1);
00209 izfor=1;
00210 }
00211 else {
00212 stripItr.SetDirection(-1);
00213 izfor=-1;
00214 }
00215 stripItr.Reset();
00216 while (CandStripHandle *strip = stripItr()) {
00217 cfh.AddDaughterLink(*strip);
00218 cfh.SetInShower(strip,track0->IsInShower(strip));
00219 }
00220 for (Int_t iplane = cfh.GetEndPlane(); iplane*fDirection>=cfh.GetVtxPlane()*fDirection; iplane-=fDirection) {
00221 if (track0->IsTPosValid(iplane)) {
00222 // u and v coordinates have been calculated for this plane
00223 cfh.SetRange(iplane,track0->GetRange(iplane));
00224
00225 cfh.SetU(iplane,track0->GetU(iplane));
00226 cfh.SetV(iplane,track0->GetV(iplane));
00227 cfh.SetdS(iplane,track0->GetdS(iplane));
00228 }
00229 }
00230 Double_t range = cfh.GetRange();
00231 if (range>0.) {
00232 cfh.SetMomentumRange(0.105658*exp(1.0363*log(range)-4.383));
00233 // taken from PDG R/M vs p/M plot for Fe
00234 } else {
00235 cfh.SetMomentumRange(0.);
00236 }
00237
00238 }
00239 else {
00240 //All the maps have been filled with the initial values.
00241 cfh.SetU0Initial(fU0);
00242 cfh.SetV0Initial(fV0);
00243 cfh.SetdUdZ0Initial(fdUdZ0);
00244 cfh.SetdVdZ0Initial(fdVdZ0);
00245 cfh.SetP0Initial(fP0);
00246 if(fP0>0) {
00247 cfh.SetInitialQP(fCharge/fP0);
00248 }
00249
00250
00251 SwimAndFillMaps();
00252 FillWeightMatrix();
00253 FillMatricesAandC();
00254 InvertMatrixAndGetParams();
00255 //First iteration don't know if we have the correct sign.
00256 //But the clever algorithm switches sign all by itself.
00257 SwimAndFillMaps();
00258 FillWeightMatrix();
00259 FillMatricesAandC();
00260 InvertMatrixAndGetParams();
00261 //Set Parameters in cfh
00262 cfh.SetU0(fU0);
00263 cfh.SetV0(fV0);
00264 cfh.SetdUdZ0(fdUdZ0);
00265 cfh.SetdVdZ0(fdVdZ0);
00266 cfh.SetP0(fP0);
00267
00268 cfh.SetU0Err(fDU0);
00269 cfh.SetV0Err(fDV0);
00270 cfh.SetdUdZ0Err(fDdUdZ0);
00271 cfh.SetdVdZ0Err(fDdVdZ0);
00272 cfh.SetP0Err(fDP0);
00273 int numIterations=1;
00274 int numIterations2=0;
00275 while((fabs(fLastP0-fP0)>0.2 || fIsNaturalP0==0 )&& numIterations<=fParmMaxIterate) {
00276 SwimAndFillMaps();
00277 FillWeightMatrix();
00278 FillMatricesAandC();
00279 InvertMatrixAndGetParams();
00280 //Set Parameters in cfh
00281 cfh.SetU0(fU0);
00282 cfh.SetV0(fV0);
00283 cfh.SetdUdZ0(fdUdZ0);
00284 cfh.SetdVdZ0(fdVdZ0);
00285 cfh.SetP0(fP0);
00286
00287 cfh.SetU0Err(fDU0);
00288 cfh.SetV0Err(fDV0);
00289 cfh.SetdUdZ0Err(fDdUdZ0);
00290 cfh.SetdVdZ0Err(fDdVdZ0);
00291 cfh.SetP0Err(fDP0);
00292 numIterations++;
00293 }
00294 if(fabs(fLastP0-fP0)>0.2 && fIsNaturalP0==1) {
00295 fP0=(fLastP0+fP0)/2.0;
00296
00297 while((fabs(fLastP0-fP0)>0.2 || fIsNaturalP0==0 )&& numIterations2<=fParmMaxIterate) {
00298 SwimAndFillMaps();
00299 FillWeightMatrix();
00300 FillMatricesAandC();
00301 InvertMatrixAndGetParams();
00302 //Set Parameters in cfh
00303 cfh.SetU0(fU0);
00304 cfh.SetV0(fV0);
00305 cfh.SetdUdZ0(fdUdZ0);
00306 cfh.SetdVdZ0(fdVdZ0);
00307 cfh.SetP0(fP0);
00308
00309 cfh.SetU0Err(fDU0);
00310 cfh.SetV0Err(fDV0);
00311 cfh.SetdUdZ0Err(fDdUdZ0);
00312 cfh.SetdVdZ0Err(fDdVdZ0);
00313 cfh.SetP0Err(fDP0);
00314 numIterations2++;
00315 }
00316 }
00317 cfh.SetNIterate(numIterations2+numIterations);
00318
00319 //Quick check to see if our answer makes any sense what so ever
00320 Double_t totalRange=0;
00321 for (int i=0;i<=fTrkClstList.GetLast(); i++) {
00322 TrackClusterSR *tc = dynamic_cast<TrackClusterSR*>(fTrkClstList.At(i));
00323 Int_t iplane=tc->GetPlane();
00324
00325 IntDoubleMap::iterator pos;
00326 pos=fRangeThisPlane.find(iplane);
00327 if(pos!=fRangeThisPlane.end()) {
00328 totalRange+=pos->second;
00329 }
00330 }
00331 if(totalRange < track0->GetRange() / 3.0) {
00332 //Something messed up.
00333 CandStripHandleItr stripItr(track0->GetDaughterIterator());
00334 CandStripHandleKeyFunc *stripKf = stripItr.CreateKeyFunc();
00335 stripKf->SetFun(CandStripHandle::KeyFromPlane);
00336 stripItr.GetSet()->AdoptSortKeyFunc(stripKf);
00337 stripKf = 0;
00338 Int_t izfor=1;
00339 if (track0->GetTimeSlope()>0.) {
00340 stripItr.SetDirection(1);
00341 izfor=1;
00342 }
00343 else {
00344 stripItr.SetDirection(-1);
00345 izfor=-1;
00346 }
00347 stripItr.Reset();
00348 while (CandStripHandle *strip = stripItr()) {
00349 cfh.AddDaughterLink(*strip);
00350 cfh.SetInShower(strip,track0->IsInShower(strip));
00351 }
00352 for (Int_t iplane = cfh.GetEndPlane(); iplane*fDirection>=cfh.GetVtxPlane()*fDirection; iplane-=fDirection) {
00353 if (track0->IsTPosValid(iplane)) {
00354 // u and v coordinates have been calculated for this plane
00355 cfh.SetRange(iplane,track0->GetRange(iplane));
00356
00357 cfh.SetU(iplane,track0->GetU(iplane));
00358 cfh.SetV(iplane,track0->GetV(iplane));
00359 cfh.SetdS(iplane,track0->GetdS(iplane));
00360 }
00361 }
00362 Double_t range = cfh.GetRange();
00363 if (range>0.) {
00364 cfh.SetMomentumRange(0.105658*exp(1.0363*log(range)-4.383));
00365 // taken from PDG R/M vs p/M plot for Fe
00366 } else {
00367 cfh.SetMomentumRange(0.);
00368 }
00369 }
00370 else {
00371
00372 //Whoo! have hopefully fitted the little blighter
00373 for (int i=0;i<=fTrkClstList.GetLast(); i++) {
00374 TrackClusterSR *tc = dynamic_cast<TrackClusterSR*>(fTrkClstList.At(i));
00375 for (int j=0; j<=tc->GetStripList()->GetLast(); j++) {
00376 CandStripHandle *strip = dynamic_cast<CandStripHandle*>
00377 (tc->GetStripList()->At(j));
00378 cfh.AddDaughterLink(*strip);
00379 cfh.SetInShower(strip,track0->IsInShower(strip));
00380 Int_t thisPlane=tc->GetPlane();
00381 if((fDirection==1 && thisPlane==fLowestPlane) ||
00382 (fDirection==-1 && thisPlane==fHighestPlane)) {
00383 cfh.SetU(thisPlane,fU0);
00384 cfh.SetV(thisPlane,fV0);
00385 }
00386 IntDoubleMap::iterator pos;
00387 pos=fSwamU.find(thisPlane);
00388 if(pos!=fSwamU.end()) cfh.SetU(thisPlane,pos->second);
00389 pos=fSwamV.find(thisPlane);
00390 if(pos!=fSwamV.end()) cfh.SetV(thisPlane,pos->second);
00391 pos=fSwamdU.find(thisPlane);
00392 if(pos!=fSwamdU.end()) cfh.SetdUdZ(thisPlane,pos->second);
00393 pos=fSwamdV.find(thisPlane);
00394 if(pos!=fSwamdV.end()) cfh.SetdVdZ(thisPlane,pos->second);
00395
00396
00397 }
00398 }
00399
00400 // set ds, range
00401 Double_t totalRange=0;
00402 for (Int_t iplane = cfh.GetEndPlane(); iplane*fDirection>=cfh.GetVtxPlane()*fDirection; iplane-=fDirection) {
00403 if (track0->IsTPosValid(iplane)) {
00404 // u and v coordinates have been calculated for this plane
00405 cfh.SetdS(iplane,track0->GetdS(iplane));
00406 cfh.SetRange(iplane,totalRange);
00407 IntDoubleMap::iterator pos;
00408 pos=fRangeThisPlane.find(iplane);
00409 if(pos!=fRangeThisPlane.end()) {
00410 totalRange+=pos->second;
00411 }
00412 }
00413 }
00414 MSG("FitTrack3", Msg::kVerbose) << "Total Range: " << totalRange
00415 << " CandTrackSR range: "
00416 << track0->GetRange() << endl;
00417 int definetlyFail=0;
00418 if(totalRange < track0->GetRange() / 3.0) {
00419 definetlyFail=1;
00420 for (Int_t iplane = cfh.GetEndPlane();
00421 iplane*fDirection>=cfh.GetVtxPlane()*fDirection;
00422 iplane-=fDirection) {
00423 if (track0->IsTPosValid(iplane)) {
00424 // u and v coordinates have been calculated for this plane
00425 cfh.SetRange(iplane,track0->GetRange(iplane));
00426 }
00427 }
00428 }
00429
00430
00431
00432
00433
00434 MSG("FitTrack3", Msg::kVerbose) << "Finished after "
00435 << numIterations << " iterations." << endl;
00436 if(fabs(fLastP0-fP0)<0.2 && fIsNaturalP0==1 && fP0>0.1 && !definetlyFail) //Sometimes get silly results
00437 cfh.SetPass(1);
00438 MSG("FitTrack3", Msg::kVerbose) << "Did it pass?? "
00439 << cfh.GetPass() << endl;
00440
00441
00442 Double_t range = cfh.GetRange();
00443 if (range>0.) {
00444 cfh.SetMomentumRange(0.105658*exp(1.0363*log(range)-4.383));
00445 // taken from PDG R/M vs p/M plot for Fe
00446 } else {
00447 cfh.SetMomentumRange(0.);
00448 }
00449
00450
00451 MSG("FitTrack3", Msg::kVerbose) << "Set Momentum Range to: "
00452 << cfh.GetMomentumRange() << endl;
00453 cfh.SetMomentumCurve(fP0);
00454 MSG("FitTrack3", Msg::kVerbose) << "Set Momentum Curve to: "
00455 << cfh.GetMomentumCurve() << endl;
00456 cfh.SetEMCharge(fCharge);
00457 MSG("FitTrack3", Msg::kVerbose) << "Set Charge to: "
00458 << cfh.GetEMCharge() << endl;
00459 cfh.SetChi2(fChiSqU+fChiSqV);
00460 MSG("FitTrack3", Msg::kVerbose) << "Set Chisq to: "
00461 << cfh.GetChi2() << endl;
00462
00463 }
00464 }
00465 CleanUp();
00466 fFile->cd();
00467 fReasonsForNotFitting->Write("reasons",TObject::kOverwrite);
00468 }
|
|
|
Definition at line 1026 of file AlgFitTrack3.cxx. References PlexPlaneId::AsString(), fActualZ, fDistanceFromStart, SwimObj3::fdudz, SwimObj3::fdvdz, fErrorOnMeasure, fFile, fIterations, fMeasuredU, fMeasuredV, fMyVC, fRangeThisPlane, fS_isU, fS_measuredU, fS_measuredV, fS_swamU, fS_swamV, fS_z, fSillyTree, fSwamdU, fSwamdV, fSwamU, fSwamV, fTheta0i, fThicknessi, fTrkClstList, SwimObj3::fu, SwimObj3::fv, SwimObj3::GetPath(), PlexPlaneId::GetPlane(), TrackClusterSR::GetPlane(), UgliGeomHandle::GetPlaneIdFromZ(), TrackClusterSR::GetPlaneView(), SwimObj3::GetRange(), SwimObj3::GetTheta(), TrackClusterSR::GetTPos(), TrackClusterSR::GetTPosError(), TrackClusterSR::GetZPos(), PlexPlaneId::IsSteel(), MSG, and SwimObj3::SwimTo(). Referenced by RunAlg(). 01026 {
01027 MSG("FitTrack3", Msg::kVerbose)
01028 << "AlgFitTrack3::SwimAndFillMaps()" << endl;
01029 // Have now initialised the following variables.
01030 //
01031 fIterations++;
01032
01033
01034 Double_t startU=fU0;
01035 Double_t startV=fV0;
01036 Double_t startZ=fZ0;
01037 Double_t startdUdZ=fdUdZ0;
01038 Double_t startdVdZ=fdVdZ0;
01039 Double_t startMomentum=fP0;
01040 Int_t charge=fCharge;
01041 Int_t izfor=fDirection;
01042
01043 UgliGeomHandle ugh(*fMyVC);
01044
01045
01046 fTheta0i.clear(); // mean square scattering angle.
01047 fThicknessi.clear(); // thickness .
01048 fSwamU.clear();
01049 fSwamV.clear();
01050 fSwamdU.clear();
01051 fSwamdV.clear();
01052 fMeasuredU.clear();
01053 fMeasuredV.clear();
01054 fActualZ.clear();
01055 fErrorOnMeasure.clear();
01056 fDistanceFromStart.clear();
01057 fRangeThisPlane.clear();
01058
01059 //Make ourselves a swimmer object.
01060 SwimObj3 swimmer(startU,startV,startZ,startdUdZ,startdVdZ,
01061 double(charge*startMomentum),izfor,fMyVC);
01062
01063 // Double_t lastU=startU;
01064 // Double_t lastV=startV;
01065 // Double_t lastZ=startZ;
01066
01067 // int numDone;
01068 Double_t totalDistanceFromStart=0;
01069 for(Int_t i=1; i<=fTrkClstList.GetLast(); i++) {
01070 TrackClusterSR *thisPlane =
01071 dynamic_cast<TrackClusterSR*>(fTrkClstList.At(i));
01072 Int_t newPlane=thisPlane->GetPlane();
01073 Double_t newZ=thisPlane->GetZPos();
01074 MSG("FitTrack3", Msg::kVerbose)
01075 << "Plane: " << newPlane
01076 << "\tZ Pos: " << newZ << endl;
01077
01078 fActualZ[newPlane]=newZ;
01079 swimmer.SwimTo(newZ);
01080 Double_t pathFromLast=swimmer.GetPath();
01081
01082 totalDistanceFromStart+=pathFromLast;
01083 fDistanceFromStart[newPlane]=pathFromLast;
01084 Double_t thisTheta=swimmer.GetTheta();
01085 fTheta0i[newPlane]=thisTheta;
01086 MSG("FitTrack3", Msg::kVerbose)
01087 << "Plane: " << newPlane << endl
01088 << "Z Pos: " << newZ << endl
01089 << "Path From Last Plane " << pathFromLast << endl
01090 << "Theta0i " << thisTheta << endl;
01091
01092 PlexPlaneId planeId= ugh.GetPlaneIdFromZ(newZ);
01093 MSG("FitTrack3", Msg::kVerbose)
01094 << "Plane: " << newPlane << endl
01095 << "From Plane Id: " << planeId.GetPlane() << endl
01096 << "IsSteel: " << planeId.IsSteel() << endl
01097 << planeId.AsString() << endl;
01098 fRangeThisPlane[newPlane]=swimmer.GetRange();
01099
01100 Double_t thisPlaneThickness=pathFromLast; //Seems I need this number and not xi
01101
01102 fThicknessi[newPlane]=thisPlaneThickness;
01103
01104 fErrorOnMeasure[newPlane]=thisPlane->GetTPosError();
01105 fS_z=newZ;
01106 fSwamdV[newPlane]=swimmer.fdvdz;
01107 fSwamdU[newPlane]=swimmer.fdudz;
01108 fSwamV[newPlane]=swimmer.fv;
01109 fSwamU[newPlane]=swimmer.fu;
01110 switch (thisPlane->GetPlaneView()) {
01111 case PlaneView::kU:
01112 fMeasuredV[newPlane]=thisPlane->GetTPos();
01113 fS_measuredV=thisPlane->GetTPos();
01114 fS_swamV=swimmer.fv;
01115 fS_isU=0;
01116 break;
01117 case PlaneView::kV:
01118 fMeasuredU[newPlane]=thisPlane->GetTPos();
01119 fS_measuredU=thisPlane->GetTPos();
01120 fS_swamU=swimmer.fu;
01121 fS_isU=1;
01122 break;
01123 default:
01124 break;
01125 }
01126 fSillyTree->Fill();
01127 }
01128
01129 fFile->cd();
01130 fSillyTree->Write("sillyTree",TObject::kOverwrite);
01131
01132
01133 }
|
|
|
Reimplemented from AlgBase. Definition at line 1806 of file AlgFitTrack3.cxx. 01807 {
01808 }
|
|
|
Definition at line 63 of file AlgFitTrack3.h. Referenced by CleanUp(), FillMatricesAandC(), SwimAndFillMaps(), and ~AlgFitTrack3(). |
|
|
Definition at line 74 of file AlgFitTrack3.h. Referenced by InitializeTrkClstList(), InvertMatrixAndGetParams(), and RunAlg(). |
|
|
Definition at line 99 of file AlgFitTrack3.h. Referenced by FillMatricesAandC(), and RunAlg(). |
|
|
Definition at line 100 of file AlgFitTrack3.h. Referenced by FillMatricesAandC(), and RunAlg(). |
|
|
Definition at line 93 of file AlgFitTrack3.h. Referenced by InvertMatrixAndGetParams(), and RunAlg(). |
|
|
Definition at line 94 of file AlgFitTrack3.h. Referenced by InvertMatrixAndGetParams(), and RunAlg(). |
|
|
Definition at line 70 of file AlgFitTrack3.h. Referenced by CalculateD(), CalculateP(), FillMatricesAandC(), FillWeightMatrix(), InitializeTrkClstList(), and RunAlg(). |
|
|
Definition at line 53 of file AlgFitTrack3.h. Referenced by CalculateD(), CleanUp(), SwimAndFillMaps(), and ~AlgFitTrack3(). |
|
|
Definition at line 90 of file AlgFitTrack3.h. Referenced by InvertMatrixAndGetParams(), and RunAlg(). |
|
|
Definition at line 91 of file AlgFitTrack3.h. Referenced by InvertMatrixAndGetParams(), and RunAlg(). |
|
|
Definition at line 78 of file AlgFitTrack3.h. Referenced by FillMatricesAandC(), InitializeTrkClstList(), InvertMatrixAndGetParams(), and RunAlg(). |
|
|
Definition at line 92 of file AlgFitTrack3.h. Referenced by InvertMatrixAndGetParams(), and RunAlg(). |
|
|
Definition at line 79 of file AlgFitTrack3.h. Referenced by FillMatricesAandC(), InitializeTrkClstList(), InvertMatrixAndGetParams(), and RunAlg(). |
|
|
Definition at line 60 of file AlgFitTrack3.h. Referenced by CleanUp(), FillWeightMatrix(), SwimAndFillMaps(), and ~AlgFitTrack3(). |
|
|
Definition at line 131 of file AlgFitTrack3.h. Referenced by RunAlg(), SwimAndFillMaps(), and ~AlgFitTrack3(). |
|
|
Definition at line 66 of file AlgFitTrack3.h. Referenced by CalculateP(), and InitializeTrkClstList(). |
|
|
Definition at line 69 of file AlgFitTrack3.h. Referenced by CalculateP(), FillMatricesAandC(), FillWeightMatrix(), and InitializeTrkClstList(). |
|
|
Definition at line 72 of file AlgFitTrack3.h. Referenced by InvertMatrixAndGetParams(), and RunAlg(). |
|
|
Definition at line 102 of file AlgFitTrack3.h. Referenced by CleanUp(), and SwimAndFillMaps(). |
|
|
Definition at line 82 of file AlgFitTrack3.h. Referenced by InvertMatrixAndGetParams(). |
|
|
Definition at line 86 of file AlgFitTrack3.h. Referenced by InvertMatrixAndGetParams(). |
|
|
Definition at line 87 of file AlgFitTrack3.h. Referenced by InvertMatrixAndGetParams(). |
|
|
Definition at line 88 of file AlgFitTrack3.h. Referenced by InvertMatrixAndGetParams(), and RunAlg(). |
|
|
Definition at line 67 of file AlgFitTrack3.h. Referenced by InitializeTrkClstList(). |
|
|
Definition at line 84 of file AlgFitTrack3.h. Referenced by InvertMatrixAndGetParams(). |
|
|
Definition at line 85 of file AlgFitTrack3.h. Referenced by InvertMatrixAndGetParams(). |
|
|
Definition at line 83 of file AlgFitTrack3.h. Referenced by InvertMatrixAndGetParams(). |
|
|
Definition at line 68 of file AlgFitTrack3.h. Referenced by CalculateP(), FillMatricesAandC(), FillWeightMatrix(), and InitializeTrkClstList(). |
|
|
Definition at line 125 of file AlgFitTrack3.h. Referenced by FillMatricesAandC(), and InvertMatrixAndGetParams(). |
|
|
Definition at line 126 of file AlgFitTrack3.h. Referenced by FillMatricesAandC(), and InvertMatrixAndGetParams(). |
|
|
Definition at line 61 of file AlgFitTrack3.h. Referenced by CleanUp(), FillMatricesAandC(), FillWeightMatrix(), SwimAndFillMaps(), and ~AlgFitTrack3(). |
|
|
Definition at line 62 of file AlgFitTrack3.h. Referenced by CleanUp(), FillMatricesAandC(), FillWeightMatrix(), SwimAndFillMaps(), and ~AlgFitTrack3(). |
|
|
Definition at line 73 of file AlgFitTrack3.h. Referenced by InvertMatrixAndGetParams(), and RunAlg(). |
|
|
Definition at line 129 of file AlgFitTrack3.h. Referenced by InitializeTrkClstList(), and SwimAndFillMaps(). |
|
|
Definition at line 80 of file AlgFitTrack3.h. Referenced by FillMatricesAandC(), InitializeTrkClstList(), InvertMatrixAndGetParams(), and RunAlg(). |
|
|
Definition at line 48 of file AlgFitTrack3.h. Referenced by RunAlg(). |
|
|
Definition at line 49 of file AlgFitTrack3.h. Referenced by InitializeTrkClstList(), and RunAlg(). |
|
|
Definition at line 137 of file AlgFitTrack3.h. Referenced by InitializeTrkClstList(). |
|
|
Definition at line 138 of file AlgFitTrack3.h. Referenced by InitializeTrkClstList(). |
|
|
Definition at line 135 of file AlgFitTrack3.h. Referenced by InitializeTrkClstList(). |
|
|
Definition at line 136 of file AlgFitTrack3.h. Referenced by InitializeTrkClstList(). |
|
|
Definition at line 139 of file AlgFitTrack3.h. Referenced by InitializeTrkClstList(). |
|
|
Definition at line 140 of file AlgFitTrack3.h. Referenced by InitializeTrkClstList(). |
|
|
Definition at line 52 of file AlgFitTrack3.h. Referenced by CleanUp(), RunAlg(), SwimAndFillMaps(), and ~AlgFitTrack3(). |
|
|
Definition at line 133 of file AlgFitTrack3.h. Referenced by InitializeTrkClstList(), and RunAlg(). |
|
|
Definition at line 134 of file AlgFitTrack3.h. Referenced by InitializeTrkClstList(). |
|
|
Definition at line 105 of file AlgFitTrack3.h. Referenced by CleanUp(), FillWeightMatrix(), GetWiju(), and ~AlgFitTrack3(). |
|
|
Definition at line 106 of file AlgFitTrack3.h. Referenced by CleanUp(), FillWeightMatrix(), GetWijv(), and ~AlgFitTrack3(). |
|
|
Definition at line 148 of file AlgFitTrack3.h. Referenced by SwimAndFillMaps(). |
|
|
Definition at line 144 of file AlgFitTrack3.h. Referenced by SwimAndFillMaps(). |
|
|
Definition at line 145 of file AlgFitTrack3.h. Referenced by SwimAndFillMaps(). |
|
|
Definition at line 146 of file AlgFitTrack3.h. Referenced by SwimAndFillMaps(). |
|
|
Definition at line 147 of file AlgFitTrack3.h. Referenced by SwimAndFillMaps(). |
|
|
Definition at line 143 of file AlgFitTrack3.h. Referenced by SwimAndFillMaps(). |
|
|
Definition at line 132 of file AlgFitTrack3.h. Referenced by SwimAndFillMaps(). |
|
|
Definition at line 103 of file AlgFitTrack3.h. Referenced by FillWeightMatrix(), and InitializeTrkClstList(). |
|
|
Definition at line 104 of file AlgFitTrack3.h. Referenced by FillWeightMatrix(), and InitializeTrkClstList(). |
|
|
Definition at line 58 of file AlgFitTrack3.h. Referenced by CleanUp(), RunAlg(), SwimAndFillMaps(), and ~AlgFitTrack3(). |
|
|
Definition at line 59 of file AlgFitTrack3.h. Referenced by CleanUp(), RunAlg(), SwimAndFillMaps(), and ~AlgFitTrack3(). |
|
|
Definition at line 56 of file AlgFitTrack3.h. Referenced by CleanUp(), FillMatricesAandC(), RunAlg(), SwimAndFillMaps(), and ~AlgFitTrack3(). |
|
|
Definition at line 57 of file AlgFitTrack3.h. Referenced by CleanUp(), FillMatricesAandC(), RunAlg(), SwimAndFillMaps(), and ~AlgFitTrack3(). |
|
|
Definition at line 54 of file AlgFitTrack3.h. Referenced by CalculateP(), CleanUp(), SwimAndFillMaps(), and ~AlgFitTrack3(). |
|
|
Definition at line 55 of file AlgFitTrack3.h. Referenced by CalculateP(), CleanUp(), SwimAndFillMaps(), and ~AlgFitTrack3(). |
|
|
Definition at line 120 of file AlgFitTrack3.h. Referenced by RunAlg(). |
|
|
Definition at line 128 of file AlgFitTrack3.h. Referenced by CleanUp(), InitializeTrkClstList(), RunAlg(), SwimAndFillMaps(), and ~AlgFitTrack3(). |
|
|
Definition at line 76 of file AlgFitTrack3.h. Referenced by FillMatricesAandC(), InitializeTrkClstList(), InvertMatrixAndGetParams(), and RunAlg(). |
|
|
Definition at line 123 of file AlgFitTrack3.h. Referenced by CleanUp(), and FillWeightMatrix(). |
|
|
Definition at line 77 of file AlgFitTrack3.h. Referenced by FillMatricesAandC(), InitializeTrkClstList(), InvertMatrixAndGetParams(), and RunAlg(). |
|
|
Definition at line 124 of file AlgFitTrack3.h. Referenced by CleanUp(), and FillWeightMatrix(). |
|
|
Definition at line 75 of file AlgFitTrack3.h. Referenced by InitializeTrkClstList(), and InvertMatrixAndGetParams(). |
1.3.9.1