//----------------------------------------------------------- // File and Version Information: // $Id$ // // Description: // TPC Helitron/Plawa track fitting routine // // // Environment: // Processing of data from FOPI experiment S339 // // Author List: // Paul Buehler SMI/OEAW // //----------------------------------------------------------- // This Class' Header ------------------ #include "TpcHelTrackFittingTaskPhilipp.h" // C++ headers // Collaborating Class Headers -------- #include "FairRun.h" #include "FairRunAna.h" #include "FairRuntimeDb.h" #include "TVectorD.h" #include "TH2.h" #include "TCanvas.h" #include "TGraph.h" #include "TMultiGraph.h" #include "TLatex.h" #include "TLegend.h" #include "RKTrackRep.h" #include "GFTrack.h" #include "GFTrackCand.h" #include "RecoHits/GFAbsRecoHit.h" #include "PseudoSpacePoint.h" #include "GFKalman.h" #include "FOPIField.h" #include "FopiEvent.h" #include "CdcHit.h" #include "HeliHit.h" #include "HeliTrack.h" #include "PlawaTrack.h" #include "FopiForwardHit.h" #include "GFException.h" #define DEBUG 0 //----------------------------------------------------------- TpcHelTrackFittingTaskPhilipp::TpcHelTrackFittingTaskPhilipp() : fFopiEvBranchName("FopiEvent"), fTPCHelMatchParsBranchName("TpcHelMatchPars"), fTPCHelTrackBranchName("TpcHelPreFit"), fTPCHelTrackFitBranchName("TpcHelPostFit"), fTPCHelFitParBranchName("TpcHelFitPars"), fPersistence(kTRUE), fVerbose(kFALSE) { /*, frep(0) */ fRad = 3.1415926 / 180.; fTheMean = 1.352; fTheSigma = 1.68; fPhiMean = 8.698; fPhiSigma = 3.913; } //----------------------------------------------------------- TpcHelTrackFittingTaskPhilipp::~TpcHelTrackFittingTaskPhilipp() { } //----------------------------------------------------------- InitStatus TpcHelTrackFittingTaskPhilipp::Init() { //Get ROOT Manager ------------------ FairRootManager *ioman = FairRootManager::Instance(); if (ioman == 0) { Error("TpcHelTrackFittingTaskPhilipp::Init", "RootManager not instantiated!"); return kERROR; } // Get input ------------------------ fFopiEventArray = (TClonesArray *) ioman->GetObject(fFopiEvBranchName); if (fFopiEventArray == 0) { Error("TpcHelTrackFittingTaskPhilipp::Init", "FOPI event not found!"); return kERROR; } fMatchingPars = (TClonesArray *) ioman->GetObject(fTPCHelMatchParsBranchName); if (fMatchingPars == 0) { Error("TpcHelTrackFittingTaskPhilipp::Init", "Matching parameters array not found!!"); return kERROR; } fTPCHelTracks = (TClonesArray *) ioman->GetObject(fTPCHelTrackBranchName); if (fTPCHelTracks == 0) { Error("TpcHelTrackFittingTaskPhilipp::Init", "TPCHelTrack array not found!!"); return kERROR; } // prepare trees for output fTPCHelTrackFits = new TClonesArray("GFTrack"); ioman->Register(fTPCHelTrackFitBranchName, "Tpc", fTPCHelTrackFits, fPersistence); FitPars = new Double_t[18]; FitParsVec = new TVectorD(18); fFitPars = new TClonesArray("TVectorD"); ioman->Register(fTPCHelFitParBranchName, "Tpc", fFitPars, fPersistence); return kSUCCESS; } //----------------------------------------------------------- void TpcHelTrackFittingTaskPhilipp::SetParContainers() { // std::cout<<"TpcResidualTask::SetParContainers"<GetRuntimeDb(); if (!db) Fatal("SetParContainers", "No runtime database"); } //----------------------------------------------------------- bool fitsorter(TVector3 h1, TVector3 h2) { return(h1.Y() < h2.Y()); } void TpcHelTrackFittingTaskPhilipp::Exec(Option_t * opt) { if (fVerbose) fprintf(stderr, "TpcHelTrackFittingTaskPhilipp:Exec\n"); // Initialisations std::vector map; std::vector map2; TVector3 entry; TVector3 entry2; TVector3 pos0, mom0; TVector3 pos, mom; TVector3 poca, direction; TVector3 poca2, direction2; TMatrixT < double >cov0; TMatrixT < double >cov; TVectorD *MatchingParsVec = new TVectorD(23); GFDetPlane plane0 = GFDetPlane(TVector3(0., 0., 0.), TVector3(1., 0., 0.), TVector3(0., 1., 0.)); GFDetPlane plane1; GFAbsTrackRep *rep = NULL; GFAbsTrackRep *rep2 = NULL; // clean up result buffer if (fTPCHelTrackFits == 0) Fatal("TpcHelMatchingTask::Exec()", "No fTPCHelTrackFits"); fTPCHelTrackFits->Delete(); if (fFitPars == 0) Fatal("TpcHelMatchingTask::Exec()", "No fFitPars"); fFitPars->Delete(); // get FOPI event FopiEvent *fopiEv = (FopiEvent *) fFopiEventArray->At(0); if (fopiEv == 0) return; Int_t nrun = fopiEv->GetRunNum(); Int_t nspill = fopiEv->GetSpillNum(); Int_t nevent = fopiEv->GetEvtNum(); if (fVerbose) { fprintf(stderr, "Event: %i %i %i\n", nrun, nspill, nevent); } // how many tracks are there Int_t NumMatchPars = fTPCHelTracks->GetEntries(); if (fVerbose) fprintf(stderr, "Number of track combinations: %i\n", NumMatchPars); // prepare fitter GFKalman fitter; fitter.setNumIterations(10); // loop over matching Pars map.clear(); map2.clear(); Int_t ntracks = -1; Int_t counter = -1; Double_t themin = fTheMean - 3 * fTheSigma; Double_t themax = fTheMean + 5 * fTheSigma; Double_t phimin = 0;//fPhiMean - 3 * fPhiSigma; Double_t phimax = 60;//fPhiMean + 3 * fPhiSigma; // GFTrack **fitbuffer = new GFTrack*[NumMatchPars]; for (Int_t ii = 0; ii < NumMatchPars; ii++) { MatchingParsVec = (TVectorD *) fMatchingPars->At(ii); Double_t *MatchingPars = MatchingParsVec->GetMatrixArray(); Double_t deltaThe = MatchingPars[6] - MatchingPars[14]; Double_t deltaphi = MatchingPars[9] - MatchingPars[16]; if (deltaphi < -180) deltaphi += 360; else if (deltaphi > 180) deltaphi -= 360; deltaphi *= MatchingPars[17]; // fprintf(stderr,"TPCphi: %f, Helphi: %f ,p-Hel: %f\n", MatchingPars[9], MatchingPars[16], MatchingPars[17] ); // fprintf(stderr,"TPCthe: %f, Helthe: %f\n", MatchingPars[6], MatchingPars[14]); // fprintf(stderr,"deltaThe: %f ,deltaphi: %f\n",deltaThe, deltaphi); if (deltaThe > themin & deltaThe < themax & deltaphi > phimin & deltaphi < phimax) { // load track GFTrack *trk = (GFTrack *) fTPCHelTracks->At(ii); // trk->setCardinalRep(frep); /* Int_t numHits=trk->getNumHits(); for(Int_t jj=0;jjgetHit(jj); Hit->Print(); } */// read initial conditions // rep = trk->getCardinalRep(); // rep->getPosMom(plane0,pos0,mom0); rep = trk->getTrackRep(0); pos0 = rep->getPos(); mom0 = rep->getMom(); fitter.setInitialDirection(1); // fit track try { fitter.processTrack(trk); } catch(GFException & e) { e.what(); } fprintf(stderr, "Representation %i: %i\n", 0, rep->getStatusFlag()); fprintf(stderr, "Chi2: %f/%f\n", rep->getChiSqu(), rep->getRedChiSqu()); fprintf(stderr, "Charge: %f\n", rep->getCharge()); fprintf(stderr, "Initial/Fitted momentum: %f/%f\n", mom0.Mag(), rep->getMom().Mag()); rep2 = trk->getTrackRep(1); if (rep->getStatusFlag()==0 & rep->getRedChiSqu()!=0.0) { entry = TVector3(ii,rep->getRedChiSqu(),0); map.push_back(entry); } if (rep2->getStatusFlag()==0 & rep2->getRedChiSqu()!=0.0) { entry2 = TVector3(ii,rep2->getRedChiSqu(),0); map2.push_back(entry2); } } } // sort map according to chi2 sort (map.begin(),map.end(),fitsorter); Int_t cc = 0; while (map.size() > cc) { // TPC and Helitron track indices corresponding to combination with // minimal chi2 Double_t *MP = ((TVectorD*)fMatchingPars->At((map)[cc].X()))->GetMatrixArray(); Int_t iiref = MP[3]; Int_t jjref = MP[12]; // eliminate all combination containing either TPC track iiref or // Helitron track jjref Int_t pp = cc+1; while (map.size() > pp) { Double_t *MP2 = ((TVectorD*)fMatchingPars->At((map)[pp].X()))->GetMatrixArray(); Int_t iiloc = MP2[3]; Int_t jjloc = MP2[12]; if (iiloc==iiref || jjloc==jjref) { map.erase(map.begin()+pp); } else { pp++; } } cc++; } // sort map according to chi2 sort (map2.begin(),map2.end(),fitsorter); cc = 0; while (map2.size() > cc) { // TPC and Helitron track indices corresponding to combination with // minimal chi2 Double_t *MP = ((TVectorD*)fMatchingPars->At((map2)[cc].X()))->GetMatrixArray(); Int_t iiref = MP[3]; Int_t jjref = MP[12]; // eliminate all combination containing either TPC track iiref or // Helitron track jjref Int_t pp = cc+1; while (map2.size() > pp) { Double_t *MP2 = ((TVectorD*)fMatchingPars->At((map2)[pp].X()))->GetMatrixArray(); Int_t iiloc = MP2[3]; Int_t jjloc = MP2[12]; if (iiloc==iiref || jjloc==jjref) { map2.erase(map2.begin()+pp); } else { pp++; } } cc++; } for(Int_t ii=0; iiAt(map[ii].X()); Double_t *MatchingPars = MatchingParsVec->GetMatrixArray(); GFTrack *trk = (GFTrack *) fTPCHelTracks->At(map[ii].X()); fitter.setInitialDirection(1); // fit track try { fitter.processTrack(trk); } catch(GFException & e) { e.what(); } // find POCA (Point Of Closest Approach) of 0. TVector3 mainvertex(0.,0.,0.); try { rep->extrapolateToPoint(mainvertex,poca,direction); } catch (GFException& e) { e.what(); poca=TVector3(-999,-999,-999); } FitPars[0] = 0; //rep 0 FitPars[1] = MatchingPars[3]; //TPC track index FitPars[2] = MatchingPars[12]; //Hel track index FitPars[3] = MatchingPars[10]; //momentum TPC FitPars[4] = MatchingPars[11]; //charge TPC FitPars[5] = MatchingPars[17]; //momentum Hel FitPars[6] = MatchingPars[18]; //dedx Hel FitPars[7] = rep->getMom().Mag(); //momentum TpcHel match FitPars[8] = rep->getCharge(); //charge TpcHel match FitPars[9] = rep->getRedChiSqu(); //reduced CHi2 TpcHel match FitPars[10] = rep->getStatusFlag(); //status of Fit: 0 = fit, 1 = no fit FitPars[11] = poca.X(); //X of POCA FitPars[12] = poca.Y(); //Y of POCA FitPars[13] = poca.Z(); //Z of POCA FitPars[14] = MatchingPars[6]; //The TPC track FitPars[15] = MatchingPars[14]; //The Hel track FitPars[16] = MatchingPars[19]; //Matched PLAWA hit FitPars[17] = MatchingPars[22]; //Velocity PLAWA hit counter++; FitParsVec->SetElements(FitPars); new((*fFitPars)[counter]) TVectorD(*FitParsVec); // save fitted track on output array ntracks++; new((*fTPCHelTrackFits)[ntracks]) GFTrack(*trk); } for(Int_t ii=0; iiAt(map2[ii].X()); Double_t *MatchingPars = MatchingParsVec->GetMatrixArray(); GFTrack *trk = (GFTrack *) fTPCHelTracks->At(map2[ii].X()); fitter.setInitialDirection(1); // fit track try { fitter.processTrack(trk); } catch(GFException & e) { e.what(); } // find POCA (Point Of Closest Approach) of 0. TVector3 mainvertex(0.,0.,0.); try { rep2->extrapolateToPoint(mainvertex,poca2,direction2); } catch (GFException& e) { e.what(); poca2=TVector3(-999,-999,-999); } FitPars[0] = 1; //rep2 FitPars[1] = MatchingPars[3]; //TPC track index FitPars[2] = MatchingPars[12]; //Hel track index FitPars[3] = MatchingPars[10]; //momentum TPC FitPars[4] = MatchingPars[11]; //charge TPC FitPars[5] = MatchingPars[17]; //momentum Hel FitPars[6] = MatchingPars[18]; //dedx Hel FitPars[7] = rep2->getMom().Mag(); //momentum TpcHel match FitPars[8] = rep2->getCharge(); //charge TpcHel match FitPars[9] = rep2->getRedChiSqu(); //reduced CHi2 TpcHel match FitPars[10] = rep2->getStatusFlag(); //status of Fit: 0 = fit, 1 = no fit FitPars[11] = poca2.X(); //X of POCA FitPars[12] = poca2.Y(); //Y of POCA FitPars[13] = poca2.Z(); //Z of POCA FitPars[14] = MatchingPars[6]; //The TPC track FitPars[15] = MatchingPars[14]; //The Hel track FitPars[16] = MatchingPars[19]; //Matched PLAWA hit FitPars[17] = MatchingPars[22]; //Velocity PLAWA hit counter++; FitParsVec->SetElements(FitPars); new((*fFitPars)[counter]) TVectorD(*FitParsVec); // save fitted track on output array // ntracks++; // new((*fTPCHelTrackFits)[ntracks]) GFTrack(*trk); } // delete MatchingParsVec;MatchingParsVec=NULL; } ClassImp(TpcHelTrackFittingTaskPhilipp) //-----------------------------------------------------------