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