//----------------------------------------------------------- // // Description: // 3D Straight Line fitter // // Author List: // Mathias Michel // //----------------------------------------------------------- // This Class' Header ------------------ #include "PndLmdLinFitTask.h" // C/C++ Headers ---------------------- #include // Collaborating Class Headers -------- #include "FairRootManager.h" #include "TClonesArray.h" #include "PndLinTrack.h" #include "TrackData/PndTrackCand.h" #include "PndSdsHit.h" #include "TrackData/PndTrackCandHit.h" #include "TFile.h" #include "TGeoTrack.h" #include "TGeoManager.h" #include "TLorentzVector.h" // Fit Classes ----------- #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include using namespace ROOT::Math; using namespace std; PndLmdLinFitTask::PndLmdLinFitTask() : FairTask("3D-Straight-Line-Fit") { fTCandBranchName = "LMDTrackCand"; fRecoBranchName = "LMDHitsStrip"; fTruePointBranch = "LMDPoint"; //True Points only for drawing! } PndLmdLinFitTask::~PndLmdLinFitTask() { cout<<"PndLmdLinFitTask::~PndLmdLinFitTask()"<GetObject(fTCandBranchName); if(fTCandArray==0) { Error("PndLmdLinFitTask::Init","trackcand-array not found!"); return kERROR; } fRecoArray=(TClonesArray*) ioman->GetObject(fRecoBranchName); if(fRecoArray==0) { Error("PndLmdLinFitTask::Init","reco-array not found!"); return kERROR; } fTrackArray = new TClonesArray("PndLinTrack"); ioman->Register("LMDTrack", "PndLmd", fTrackArray, kTRUE); std::cout << "-I- PndLmdLinFitTask: Initialisation successfull" << std::endl; return kSUCCESS; } void PndLmdLinFitTask::Exec(Option_t* opt) { // std::cout<<"PndLmdLinFitTask::Exec"<Delete(); Int_t ntcand=fTCandArray->GetEntriesFast(); // Detailed output if(fVerbose>1)std::cout<<" -I- PndLmdLinFitTask: contains "<2){ std::cout<< " Detailed Debug info on the candidates:"<At(itr); std::cout<< "TrackCand no. "<GetNHits()<<" hits."<GetNHits(); ihit++){ //fill Graph PndTrackCandHit theHit = trcnd->GetSortedHit(ihit); //get hit index = theHit.GetHitId(); detid = theHit.GetDetId(); std::cout << ihit << "\t" << index <20){ // std::cout<<"ntcand="<1)std::cout<<" -I- PndLmdLinFitTask: start Fitting "<At(track); const int numPts = trcnd->GetNHits(); //read how many points in this track /// Obtain first approximation ---------- TVector3 posSeed = trcnd->getPosSeed(); TVector3 dirSeed = trcnd->getDirSeed(); // if(dirSeed.Theta()<3e-2 || dirSeed.Theta()>5e-2 || fabs(dirSeed.Phi())>0.3){ // if(fVerbose>2) cout<<"Trk-cand doesn't pass throw limit (dirSeed.Theta() = "<2) std::cout << "Track: "<< track<< " Points: "<< numPts <GetSortedHit(ihit); //get hit Int_t index = theHit.GetHitId(); Int_t detId = theHit.GetDetId(); // if(fVerbose>2) std::cout << "Point: "<< ihit<< " index: "<< index <At(index); TVector3 addPos = addHit->GetPosition(); // addPos.Print(); fitme.SetPoint(ihit, addPos.X(), addPos.Y(), addPos.Z()); fitme.SetPointError(ihit, addHit->GetDx(), addHit->GetDy(), addHit->GetDz()); }//end of Hits in TCand Double_t parFit[6]; //fit-parameter TMatrixDSym *COVmatrix = new TMatrixDSym(6); Double_t accuracy = line3Dfit(numPts, &fitme, posSeed, dirSeed, parFit, COVmatrix); // if(accuracy>0 && accuracy<1e3){ PndLinTrack* trackfit = new PndLinTrack("Lumi", parFit[0], parFit[1], parFit[2], parFit[3], parFit[4], parFit[5], accuracy, firstHit, lastHit, track); trackfit->SetCovarianceMatrix(*COVmatrix); new((*fTrackArray)[rec_tkr]) PndLinTrack(*(trackfit)); //save Track delete trackfit;//TEST rec_tkr++; //} }// end of TCand's // Done-------------------------------------------------------------------------------------- if(fVerbose>2) std::cout<<"Fitting done"<SetObjectFit(gr); //min->SetFCN( *SumDistance2 ); min->SetFCN(*LocalFCN);//using local coordinate in FCN Double_t arglist[100]; arglist[0] = 1; min->ExecuteCommand("SET PRINT",arglist,1); double pStart[4] = {25.15,0.04,0.09,8.4e-5}; double pStartErr[4] = {4.5,0.004,4.5,0.004}; min->SetParameter(0,"x0",pStart[0],pStartErr[0],0,0); min->SetParameter(1,"Ax",pStart[1],pStartErr[1],0,0); min->SetParameter(2,"y0",pStart[2],pStartErr[2],0,0); min->SetParameter(3,"Ay",pStart[3],pStartErr[3],0,0); // //minimize step 1 arglist[0] = 1000; //number of functiona calls arglist[1] = 0.001; //tolerance min->ExecuteCommand("MIGRAD", arglist ,2); //if (minos) min->ExecuteCommand("MINOS",arglist,0); int nvpar,nparx; double amin,edm, errdef; min->GetStats(amin,edm,errdef,nvpar,nparx); if(fVerbose>1) min->PrintResults(1,amin); // get fit parameters and errors for (int i = 0; i <4; ++i){ fitpar[i] = min->GetParameter(i); fitparerr[i] = min->GetParError(i); } // return amin; Double_t chi2 = amin/(2.*Npoint-4); // cout<<"Chi^2 = "<2) cout<<"PndLmdLinFitTask::line3Dfit with SEED is used"<GetN(); Double_t ErrX1 = gr->GetErrorX(0); Double_t ErrY1 = gr->GetErrorY(0); Double_t ErrZ1 = gr->GetErrorZ(0); Double_t ErrX2 = gr->GetErrorX(1); Double_t ErrY2 = gr->GetErrorY(1); Double_t ErrZ2 = gr->GetErrorY(1); Double_t errRx = TMath::Hypot(ErrX1,ErrX2); Double_t errRy = TMath::Hypot(ErrY1,ErrY2); Double_t errRz = TMath::Hypot(ErrZ1,ErrZ2); //The default minimizer is Minuit TVirtualFitter::SetDefaultFitter("Minuit"); TVirtualFitter *min = TVirtualFitter::Fitter(0,5); min->SetObjectFit(gr); // min->SetFCN( *SumDistance2 ); min->SetFCN(*LocalFCN);//using local coordinate in FCN Double_t arglist[100]; arglist[0] = 1; min->ExecuteCommand("SET PRINT",arglist,1); if(fVerbose>2){ cout<<"posSeed:"<2){ for(int i=0;i<6;i++) cout<<"pStartErr["<SetParameter(0,"x0",pStart[0],pStartErr[0],0,0); min->SetParameter(1,"Ax",pStart[1],pStartErr[1],0,0); min->SetParameter(2,"y0",pStart[2],pStartErr[2],0,0); min->SetParameter(3,"Ay",pStart[3],pStartErr[3],0,0); min->SetParameter(4,"z0",pStart[4],0,0,0); min->SetParameter(5,"Az",pStart[5],0,0,0); // min->SetParameter(4,"z0",pStart[4],pStartErr[4],0,0); // min->SetParameter(5,"Az",pStart[5],pStartErr[5],0,0); // Now ready for minimization step arglist[0] = 1500; arglist[1] = 1.; min->ExecuteCommand("MIGRAD", arglist,2); ///Get results --------------------------------------------- int nvpar,nparx; double amin,edm, errdef; min->GetStats(amin,edm,errdef,nvpar,nparx); if(fVerbose>1) min->PrintResults(1,amin); Double_t fitparerr[6]; // get fit parameters for (int i = 0; i <6; ++i){ fitpar[i] = min->GetParameter(i); fitparerr[i] = min->GetParError(i); } // cout<<" before fitpar[5] = "<GetCovarianceMatrixElement(i,j); } } // (*covmatrix)(4,4) = 0.00013*gr->GetErrorZ(0)*gr->GetErrorZ(0); (*covmatrix)(4,4) = 0.0001*gr->GetErrorZ(0)*gr->GetErrorZ(0); // (*covmatrix)(4,4) = 5e-09; // (*covmatrix)(4,4) += (*covmatrix)(0,0)*pow(tan(2.326*TMath::Pi()/180.),2); double dp5_dp1 = fitpar[1]/fitpar[5]; double dp5_dp3 = fitpar[3]/fitpar[5]; double errdz2 = pow(dp5_dp1,2)*(*covmatrix)(1,1) + pow(dp5_dp3,2)*(*covmatrix)(3,3) + 2*fabs(dp5_dp1*dp5_dp3*(*covmatrix)(1,3)); // cout<<"pow(dp5_dp1,2)*(*covmatrix)(1,1) = "<