00001
00002
00003
00004
00005
00006
00007
00008
00010
00011 #include <cassert>
00012 #include <cmath>
00013 #include <set>
00014 #include <map>
00015
00016 #include "Algorithm/AlgConfig.h"
00017 #include "Candidate/CandContext.h"
00018 #include "CandFitTrack3/AlgFitTrack3.h"
00019 #include "CandFitTrack3/CandFitTrack3Handle.h"
00020 #include "CandFitTrack3/SwimObj3.h"
00021 #include "CandTrackSR/CandTrackSRHandle.h"
00022
00023 #include "CandTrackSR/TrackClusterSR.h"
00024 #include "Conventions/Mphysical.h"
00025 #include "Conventions/Munits.h"
00026 #include "Conventions/PlaneView.h"
00027 #include "DataUtil/GetMomFromRange.h"
00028 #include "MessageService/MsgService.h"
00029 #include "MinosObjectMap/MomNavigator.h"
00030 #include "Navigation/NavKey.h"
00031 #include "Navigation/NavSet.h"
00032 #include "RecoBase/CandSliceHandle.h"
00033 #include "RecoBase/CandStripHandle.h"
00034 #include "RecoBase/CandTrackHandle.h"
00035 #include "UgliGeometry/UgliGeomHandle.h"
00036 #include "Validity/VldContext.h"
00037 #include "Plex/PlexPlaneId.h"
00038
00039 #include "TFile.h"
00040 #include "TTree.h"
00041 #include "TF1.h"
00042 #include "TGraphErrors.h"
00043 #include "TMath.h"
00044
00045 ClassImp(AlgFitTrack3)
00046
00047
00048 CVSID("$Id: AlgFitTrack3.cxx,v 1.10 2007/11/11 08:37:40 rhatcher Exp $");
00049
00050
00051 AlgFitTrack3::AlgFitTrack3()
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 }
00089
00090
00091 AlgFitTrack3::~AlgFitTrack3()
00092 {
00093 MSG("FitTrack3", Msg::kDebug)
00094 << "AlgFitTrack3::~AlgFitTrack3()" << endl;
00095
00096 if(fUWeightMatrix) delete fUWeightMatrix;
00097 if(fVWeightMatrix) delete fVWeightMatrix;
00098 fTheta0i.clear();
00099 fThicknessi.clear();
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 }
00117
00118
00119 void AlgFitTrack3::CleanUp()
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();
00129 fThicknessi.clear();
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
00144
00145 }
00146
00147 void AlgFitTrack3::RunAlg(AlgConfig &ac, CandHandle &ch, CandContext &cx)
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
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
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
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
00236
00237
00238
00239 cfh.SetMomentumRange(GetMomFromRange(range));
00240
00241 } else {
00242 cfh.SetMomentumRange(0.);
00243 }
00244
00245 }
00246 else {
00247
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
00263
00264 SwimAndFillMaps();
00265 FillWeightMatrix();
00266 FillMatricesAandC();
00267 InvertMatrixAndGetParams();
00268
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
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
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
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
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
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
00372
00373
00374 cfh.SetMomentumRange(GetMomFromRange(range));
00375
00376 } else {
00377 cfh.SetMomentumRange(0.);
00378 }
00379 }
00380 else {
00381
00382
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
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
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
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)
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
00455
00456
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 }
00484
00485
00486
00487
00488 int AlgFitTrack3::InitializeTrkClstList(const CandTrackHandle *track0)
00489 {
00490 MSG("FitTrack3", Msg::kVerbose)
00491 << "AlgFitTrack3::InitializeTrkClstList()" << endl;
00492
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
00550 int lastShower=-1;
00551
00552 while (CandStripHandle *strip = stripItr()) {
00553
00554
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
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
00618
00619
00620
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
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);
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) {
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);
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);
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);
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
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
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;
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;
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;
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;
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
01026
01027
01028
01029
01030
01031
01032
01033
01034
01035
01036 return 1;
01037 }
01038
01039
01040
01041 void AlgFitTrack3::SwimAndFillMaps() {
01042 MSG("FitTrack3", Msg::kVerbose)
01043 << "AlgFitTrack3::SwimAndFillMaps()" << endl;
01044
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();
01062 fThicknessi.clear();
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
01075 SwimObj3 swimmer(startU,startV,startZ,startdUdZ,startdVdZ,
01076 double(charge*startMomentum),izfor,fMyVC);
01077
01078
01079
01080
01081
01082
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;
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 }
01149
01150
01151
01152 Double_t AlgFitTrack3::GetWiju(int row, int col) {
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 }
01167
01168
01169
01170
01171 Double_t AlgFitTrack3::GetWijv(int row, int col) {
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 }
01184
01185
01186
01187
01188
01189 void AlgFitTrack3::InvertMatrixAndGetParams() {
01190
01191 MSG("FitTrack3", Msg::kVerbose)
01192 << "AlgFitTrack3::InvertMatrixAndGetParams()" << endl;
01193
01194 fLastCharge=fCharge;
01195 fLastZ0=fZ0;
01196 fLastU0=fU0;
01197 fLastV0=fV0;
01198 fLastdUdZ0=fdUdZ0;
01199 fLastdVdZ0=fdVdZ0;
01200 fLastP0=fP0;
01201
01202
01203
01204
01205
01206 if(fMatrixA.Determinant()!=0) {
01207 TMatrixD invertedMatrixA(TMatrixD::kInverted,fMatrixA);
01208 TMatrixD solution(invertedMatrixA,TMatrixD::kMult,fMatrixC);
01209
01210
01211
01212
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 }
01296
01297
01298
01299 void AlgFitTrack3::FillMatricesAandC() {
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;
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;
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;
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;
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
01428 fChiSqU+=(newui-uim)*wiju*(newuj-ujm);
01429 fChiSqV+=(newvi-vim)*wijv*(newvj-vjm);
01430
01431
01432
01433 tempMatA[0][0]+=2.0*wiju;
01434 tempMatA[0][1]+=wiju*(zi+zj);
01435 tempMatA[0][2]+=wiju*(biu+bju);
01436
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
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
01447 tempMatA[3][2]+=wijv*(biv+bjv);
01448 tempMatA[3][3]+=2*wijv;
01449 tempMatA[3][4]+=wijv*(zi+zj);
01450
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
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
01464 for(int row=0;row<5;row++) {
01465
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 }
01473
01474
01475 void AlgFitTrack3::FillWeightMatrix() {
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
01503 mapUIndex[uptoUindex]=row;
01504 fReverseMapUIndex[row]=uptoUindex;
01505 uptoUindex++;
01506 }
01507 pos=fMeasuredV.find(planeRow);
01508 if(pos!=fMeasuredV.end()) {
01509
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
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
01609
01610
01611
01612 }
01613
01614
01615
01616 Double_t AlgFitTrack3::CalculateP(int k, int n)
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
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
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
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
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
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 }
01763
01764
01765
01766 Double_t AlgFitTrack3::CalculateD(int i, int k)
01767 {
01768 if(fDirection==1) {
01769 if(k<i) {
01770 return -1.0*CalculateD(k,i);
01771 }
01772 if(i<=fFirstPlane) return -1;
01773 if(k>fLastPlane) return -1;
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
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
01787 return 0;
01788 }
01789 else kDistance=pos->second;
01790 return kDistance-iDistance;
01791 }
01792 else if(fDirection==-1) {
01793 if(k>i) return -1.0*CalculateD(k,i);
01794 if(i>=fHighestPlane) return -1;
01795 if(k<fLowestPlane) return -1;
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
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
01809 return 0;
01810 }
01811 else kDistance=pos->second;
01812 return kDistance-iDistance;
01813
01814 }
01815 return 0;
01816 }
01817
01818
01819
01820
01821 void AlgFitTrack3::Trace(const char * ) const
01822 {
01823 }
