12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457145814591460146114621463146414651466146714681469147014711472147314741475147614771478147914801481148214831484148514861487148814891490149114921493149414951496149714981499150015011502150315041505150615071508150915101511151215131514151515161517151815191520152115221523152415251526152715281529153015311532153315341535153615371538153915401541154215431544154515461547154815491550155115521553155415551556155715581559156015611562156315641565156615671568156915701571157215731574157515761577157815791580158115821583158415851586158715881589159015911592159315941595159615971598159916001601160216031604160516061607160816091610161116121613161416151616161716181619162016211622162316241625162616271628162916301631163216331634163516361637163816391640164116421643164416451646164716481649165016511652165316541655165616571658165916601661166216631664166516661667166816691670167116721673167416751676167716781679168016811682168316841685168616871688168916901691169216931694169516961697169816991700170117021703170417051706170717081709171017111712171317141715171617171718171917201721172217231724172517261727172817291730173117321733173417351736173717381739174017411742174317441745174617471748174917501751175217531754175517561757175817591760176117621763176417651766176717681769177017711772177317741775177617771778177917801781178217831784178517861787178817891790179117921793179417951796179717981799180018011802180318041805180618071808180918101811181218131814181518161817181818191820182118221823182418251826182718281829183018311832183318341835183618371838183918401841184218431844184518461847184818491850185118521853185418551856185718581859186018611862186318641865186618671868186918701871187218731874187518761877187818791880188118821883188418851886188718881889189018911892189318941895189618971898189919001901190219031904190519061907190819091910191119121913191419151916191719181919192019211922192319241925192619271928192919301931193219331934193519361937193819391940194119421943194419451946194719481949195019511952195319541955195619571958195919601961196219631964196519661967196819691970197119721973197419751976197719781979198019811982198319841985198619871988198919901991199219931994199519961997199819992000200120022003200420052006200720082009201020112012201320142015201620172018201920202021202220232024202520262027202820292030203120322033203420352036203720382039204020412042204320442045204620472048204920502051205220532054205520562057205820592060206120622063206420652066206720682069207020712072207320742075207620772078207920802081208220832084208520862087208820892090209120922093209420952096209720982099210021012102210321042105210621072108210921102111211221132114211521162117211821192120212121222123212421252126212721282129213021312132213321342135213621372138213921402141214221432144214521462147214821492150215121522153215421552156215721582159216021612162216321642165216621672168216921702171217221732174217521762177217821792180218121822183218421852186218721882189219021912192219321942195219621972198219922002201220222032204220522062207220822092210221122122213221422152216221722182219222022212222222322242225222622272228222922302231223222332234223522362237223822392240224122422243224422452246224722482249225022512252225322542255225622572258225922602261226222632264226522662267226822692270227122722273227422752276227722782279228022812282228322842285228622872288228922902291229222932294229522962297229822992300230123022303230423052306230723082309231023112312231323142315231623172318231923202321232223232324232523262327232823292330233123322333233423352336233723382339234023412342234323442345234623472348234923502351235223532354235523562357235823592360236123622363236423652366236723682369237023712372237323742375237623772378237923802381238223832384238523862387238823892390239123922393239423952396239723982399240024012402240324042405240624072408240924102411241224132414241524162417241824192420242124222423242424252426242724282429243024312432243324342435243624372438243924402441244224432444244524462447244824492450245124522453245424552456245724582459246024612462246324642465246624672468246924702471247224732474247524762477247824792480248124822483248424852486248724882489249024912492249324942495249624972498249925002501250225032504250525062507250825092510251125122513251425152516251725182519252025212522252325242525252625272528252925302531253225332534253525362537253825392540254125422543254425452546254725482549255025512552255325542555255625572558255925602561256225632564256525662567256825692570257125722573257425752576257725782579258025812582258325842585258625872588258925902591259225932594259525962597259825992600260126022603260426052606260726082609261026112612261326142615261626172618261926202621262226232624262526262627262826292630263126322633263426352636263726382639264026412642264326442645264626472648264926502651265226532654265526562657265826592660266126622663266426652666266726682669267026712672267326742675267626772678267926802681268226832684268526862687268826892690269126922693269426952696269726982699270027012702270327042705270627072708270927102711271227132714271527162717271827192720 |
- /// \class MpdTpcKalmanFilter
- ///
- /// Kalman filter track reconstructor in the MPD central detector
- /// \author Alexander Zinchenko (LHEP, JINR, Dubna)
- #include "MpdTpcKalmanFilter.h"
- #include "MpdCodeTimer.h"
- #include "MpdTpcHit.h"
- #include "MpdTpcSectorGeo.h"
- #include "MpdTpcKalmanTrack.h"
- #include "MpdKalmanFilter.h"
- #include "MpdKalmanTrack.h"
- #include "MpdKalmanGeoScheme.h"
- #include "MpdKalmanHit.h"
- #include "MpdTpcDedxTask.h"
- #include "MpdVertexZfinder.h"
- #include "MpdTpcFoundHit.h"
- #include "MpdTpcSectorGeo.h"
- #include "TpcCluster.h"
- #include "TpcGeoPar.h"
- #include "TpcPoint.h"
- #include "TpcLheTrack.h"
- #include "MpdMCTrack.h"
- #include "FairEventHeader.h"
- #include "FairField.h"
- #include "FairMCEventHeader.h"
- #include "FairRootManager.h"
- #include "FairRunAna.h"
- #include "FairRuntimeDb.h"
- #include "FairTask.h"
- //#include "MpdTofGeoPar.h" //AZ
- #include <TClonesArray.h>
- #include <TGeoManager.h>
- #include <TGeoMatrix.h>
- #include <TGeoPgon.h>
- #include <TGeoTrd1.h>
- #include <TGeoTube.h>
- #include <TRandom.h>
- #include <TROOT.h>
- #ifdef XROOTD
- #include "TXNetFile.h"
- #endif
- #ifdef _OPENMP
- #include "omp.h"
- #include "sys/time.h"
- #endif
- const Double_t MpdTpcKalmanFilter::fgkChi2Cut = 20; //50; //20; //100;
- FILE *lunTpc = 0x0; //fopen("dl.dat","w");
- //__________________________________________________________________________
- MpdTpcKalmanFilter::MpdTpcKalmanFilter()
- : FairTask("TPC Kalman filter"),
- fNofEvents(0),
- fNTracks(0),
- fNPass(1),
- fHits(0x0),
- fKHits(new TClonesArray("MpdKalmanHit")),
- fTracks(new TClonesArray("MpdTpcKalmanTrack")),
- fTrackCand(new TClonesArray("MpdTpcKalmanTrack")),
- fLHEtracks(0x0),
- fMCtracks(0x0),
- fTpcPoints(0x0),
- fLayPointers(0x0),
- fhLays(0x0),
- fVertZ(0.0),
- fZflag(0),
- fModular(0),
- //fPadPlane(TpcPadPlane::Instance())
- //fSecGeo(MpdTpcSectorGeo::Instance())
- //fCache(new std::map<Double_t,tuple<TMatrixD,TMatrixD,TMatrixD,TMatrixD> >)
- fCache(new std::map<Double_t,matrix4>)
- {
- /// Default constructor
- fUseMCHit = kTRUE;
- }
- //__________________________________________________________________________
- MpdTpcKalmanFilter::MpdTpcKalmanFilter(const char *name,
- const char *title)
- //: FairTask(name),
- : FairTask("TPC Kalman filter"),
- fNofEvents(0),
- fNTracks(0),
- fNPass(1),
- fHits(0x0),
- fKHits(new TClonesArray("MpdKalmanHit")),
- fTracks(new TClonesArray("MpdTpcKalmanTrack")),
- fTrackCand(new TClonesArray("MpdTpcKalmanTrack")),
- fLHEtracks(0x0),
- fMCtracks(0x0),
- fTpcPoints(0x0),
- fLayPointers(0x0),
- fhLays(0x0),
- fVertZ(0.0),
- fZflag(0),
- fModular(0),
- //fPadPlane(TpcPadPlane::Instance())
- //fSecGeo(MpdTpcSectorGeo::Instance())
- //fCache(new std::map<Double_t,tuple<TMatrixD,TMatrixD,TMatrixD,TMatrixD> >)
- fCache(new std::map<Double_t,matrix4>)
- {
- /// Constructor
- FairTask *dedx = new MpdTpcDedxTask();
- Add(dedx);
- fUseMCHit = kTRUE;
- }
- //__________________________________________________________________________
- MpdTpcKalmanFilter::~MpdTpcKalmanFilter()
- {
- /// Destructor
- //new delete fHits;
- delete fKHits;
- delete fTracks;
- delete fTrackCand;
- delete [] fLayPointers;
- delete fhLays;
- }
- //__________________________________________________________________________
- InitStatus MpdTpcKalmanFilter::Init() {
- cout << "InitStatus MpdTpcKalmanFilter::Init\n\n";
- //fVerbose = 10;
- FairRootManager *manager = FairRootManager::Instance();
- fSecGeo = MpdTpcSectorGeo::Instance();
- fhLays = new TH1F("hLays","TPC layers",150,0,150);
- //fHits = (TClonesArray*) manager->GetObject("TpcCluster");
- if (fUseMCHit) fHits = (TClonesArray*) manager->GetObject("TpcHit");
- //else fHits = (TClonesArray*) manager->GetObject("MpdTpcHit");
- else fHits = (TClonesArray*) manager->GetObject("TpcRecPoint");
- if ( ! fHits ) {
- cout << "-I- "<< GetName() << "::Init: No MpdTpcHit array!" << endl;
- return kERROR;
- }
- //fKHits->SetOwner(kTRUE);
- fLHEtracks = (TClonesArray *)manager->GetObject("LheGeantTrack");
- fMCtracks = (TClonesArray *)manager->GetObject("MCTrack");
- fTpcPoints = (TClonesArray *)manager->GetObject("TpcPoint");
- #ifdef GEANT
- //fListMCtracks = (TClonesArray *)fManager->GetObject("MCTrack");
- #endif
- fNTracks = fNofEvents = 0;
- fNPass = 2; //2;
- Register();
- // Get event number offset
- /* TFile* f = new TFile(FairRootManager::Instance()->GetInFile()->GetName(),"READ"); // EL
- #ifdef XROOTD
- if (!f->IsOpen()) f = new TXNetFile(FairRootManager::Instance()->GetInFile()->GetName(),"READ");
- #endif
- cout << f->GetName() << " " << manager->GetObject("MCEventHeader.") << endl;
- TTree *cbmsim = (TTree*) f->FindObjectAny("cbmsim");
- FairMCEventHeader *head = 0x0; //new FairMCEventHeader();
- cbmsim->SetBranchAddress("MCEventHeader.",&head);
- cbmsim->GetEntry(0);
- Int_t event = head->GetEventID() - 1;
- f->Close();
- delete f;
- fNofEvents += event;
- */
- /*TH1F *hZ = */new TH1F ("hZ","Z of vertex",150,-75,75);
- /*TH1F *hPt = */new TH1F ("hPt","Pt of cand.",100,0,50);
- /*TH1F *hChi2 = */new TH1F ("hChi2","Chi2 of track",100,0,50);
- /*TH1F *hNhits = */new TH1F ("hNhits","Points per track",100,0,100);
- return kSUCCESS;
- }
- // ----- Private method SetParContainers -------------------------------
- void MpdTpcKalmanFilter::SetParContainers() {
- //return;
- // Get run and runtime database
- FairRunAna* run = FairRunAna::Instance();
- if ( ! run ) Fatal("SetParContainers", "No analysis run");
- FairRuntimeDb* db = run->GetRuntimeDb();
- if ( ! db ) Fatal("SetParContainers", "No runtime database");
- // Get TPC geometry parameter container
- db->getContainer("TpcGeoPar");
- db->getContainer("MpdTofGeoPar");
- }
- //__________________________________________________________________________
- InitStatus MpdTpcKalmanFilter::ReInit()
- {
- return kSUCCESS;
- }
- //__________________________________________________________________________
- void MpdTpcKalmanFilter::Reset()
- {
- ///
- cout << " MpdTpcKalmanFilter::Reset " << endl;
- ProcInfo_t proc;
- gSystem->GetProcInfo(&proc);
- cout << " User CPU time: " << proc.fCpuUser << " Memory: resident " << proc.fMemResident << ", virtual " << proc.fMemVirtual << endl;
- fhLays->Reset();
- fKHits->Delete();
- fTracks->Delete();
- fTrackCand->Delete();
- delete [] fLayPointers;
- fLayPointers = 0x0;
- for (Int_t i = 0; i < fgkLays; ++i)
- for (Int_t j = 0; j < fgkSecs; ++j) fLaySecBegPointers[i][j] = fLaySecEndPointers[i][j] = -1;
- fNTracks = 0;
- }
- //__________________________________________________________________________
- void MpdTpcKalmanFilter::Register()
- {
- ///
- FairTask *itsHF = (FairTask*) FairRun::Instance()->GetTask("Ideal STS hit Producer");
- TClonesArray *itsP = (TClonesArray*) FairRootManager::Instance()->GetObject("StsPoint");
- if (itsHF && itsP) FairRootManager::Instance()->
- Register("TpcKalmanTrack","MpdKalmanTrack", fTracks, kFALSE);
- else FairRootManager::Instance()->
- Register("TpcKalmanTrack","MpdKalmanTrack", fTracks, kTRUE);
- FairRootManager::Instance()->Register("TpcKalmanHit","KfHit", fKHits, kFALSE);
- }
- //__________________________________________________________________________
- void MpdTpcKalmanFilter::Finish()
- {
- ///
- fKHits->Delete();
- fTracks->Delete();
- fTrackCand->Delete();
- //((TH1F*)gROOT->FindObjectAny("hChi2"))->Write();
- //((TH1F*)gROOT->FindObjectAny("hNhits"))->Write();
- //((TH1F*)gROOT->FindObjectAny("hZ"))->Write();
- }
- //__________________________________________________________________________
- void MpdTpcKalmanFilter::Exec(Option_t * option)
- {
- ///
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
- static Int_t first = 1;
- if (first && fHits->GetEntriesFast()) {
- first = 0;
- // Check geometry
- if (fHits->First()->GetUniqueID()) fModular = 1;
- cout << "\n !!!!! ***** ***** ***** !!!!! " << endl;
- if (fModular) {
- FillGeoScheme();
- cout << " ***** TPC with modular design of readout chambers ***** " << endl;
- } else cout << " ***** TPC with cylindrical design ***** " << endl;
- cout << " !!!!! ***** ***** ***** !!!!! " << endl;
- }
- cout << "\n\n *** Event # " << ++fNofEvents << endl;
- cout << " ===== MpdTpcKalmanFilter =====\n";
- Reset();
- // Create Kalman hits
- if (fHits->GetEntriesFast() == 0) return;
- TString name = fHits->UncheckedAt(0)->GetName();
- cout << " Name: " << name << endl;
- if (name.Contains("TpcCluster")) Cluster2KalmanHits();
- else if (fModular == 0) MakeKalmanHits();
- else MakeKalmanHitsModul();
- // Evaluate vertex Z
- MpdVertexZfinder *vertexZ = (MpdVertexZfinder*) FairRun::Instance()->GetTask("MpdVertexZfinder");
- if (vertexZ) {
- vertexZ->SetHits(fKHits,fhLays);
- fVertZ = vertexZ->FindZ(fLayPointers, fZflag);
- }
- for (Int_t i = 0; i < fNPass; ++i) {
- fTrackCand->Delete();
- GetTrackSeeds(i);
- if (!fUseMCHit && i == 0) GetTrackSeedsEndCaps();
- cout << " Total number of hits for tracking: " << fHits->GetEntriesFast() << endl;
- cout << " Total number of track seeds: " << fTrackCand->GetEntriesFast() << endl;
- DoTracking(i);
- StoreTracks();
- cout << " Total number of found tracks: " << fTracks->GetEntriesFast() << endl;
- if (i != fNPass - 1) ExcludeHits(); // exclude used hits
- }
- RemoveShorts(); // remove short tracks (Nhits < 4)
- cout << " Total number of found long tracks: " << fTracks->GetEntriesFast() << endl;
- if (0) {
- // To improve dL-resolution
- Int_t ntr = fTracks->GetEntriesFast();
- for (Int_t i = 0; i < ntr; ++i) {
- MpdKalmanTrack *track = (MpdKalmanTrack*) fTracks->UncheckedAt(i);
- //track->GetCovariance()->Print();
- //track->GetParamNew()->Print();
- MpdKalmanFilter::Instance()->Refit(track, -1);
- track->SetLength(0);
- MpdKalmanFilter::Instance()->Refit(track, 1, kTRUE);
- //track->GetParamNew()->Print();
- //track->GetCovariance()->Print();
- }
- }
- AddHits(); // add hit objects to tracks
- // Apply track merging procedure (for clusters) - 11.10.2015
- if (!fUseMCHit) {
- for (Int_t ipass = 0; ipass < 3; ++ipass) {
- MergeTracks(ipass);
- }
- }
- GoToBeamLine(); // propagate tracks to the beam line
- GoOut(); // go outward for matching with outer detectors
- #ifdef GEANT
- //CheckTracks();
- #ifdef TRACK
- //PrintTracks(0);
- #endif
- #endif
-
- #ifdef MASSPROD
- //fTpcPoints->Clear("C");
- // fTstPoints->Clear("C");
- #endif
- /*
- FairMCEventHeader *mchead = (FairMCEventHeader*) FairRootManager::Instance()->GetObject("MCEventHeader.");
- Int_t event = mchead->GetEventID();
- //FairEventHeader *head = (FairEventHeader*) FairRootManager::Instance()->GetObject("EventHeader.");
- FairEventHeader *head = FairRunAna::Instance()->GetEvtHeader();
- head->SetEventId(event);
- */
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
- }
- //__________________________________________________________________________
- void MpdTpcKalmanFilter::FillGeoScheme()
- {
- /// Fill Kalman filter geometry manager info
- MpdKalmanGeoScheme *geo = MpdKalmanFilter::Instance()->GetGeo();
- TVector2 p2;
- Int_t nRows[2] = {fSecGeo->NofRowsReg(0), fSecGeo->NofRowsReg(1)};
- Int_t nSec = fSecGeo->NofSectors();
- Double_t dPhi = fSecGeo->Dphi();
- // Get size along Z
- TString volName = "tpc01sv";
- TGeoVolume *sv = gGeoManager->GetVolume(volName);
- TGeoShape *shape = sv->GetShape();
- Double_t dZ = ((TGeoTube*)shape)->GetDZ();
- if (TString(shape->ClassName()) == "TGeoTube") {
- // Old sensitive volume
- TGeoVolume *endPlateG10 = gGeoManager->GetVolume("tpc01ppS");
- dZ -= ((TGeoTube*)endPlateG10->GetShape())->GetDZ() * 2;
- TGeoVolume *endPlateAl = gGeoManager->GetVolume("tpc01bpS");
- dZ -= ((TGeoTube*)endPlateAl->GetShape())->GetDZ() * 2;
- } else if (TString(shape->ClassName()) == "TGeoPgon") {
- // New sensitive volume
- Double_t phi1 = ((TGeoPgon*)shape)->Phi1();
- Double_t loc[3] = {100, 0, 10}, glob[3] = {0};
- gGeoManager->FindNode(loc[0],loc[1],loc[2]);
- if (TString(gGeoManager->GetPath()).Contains("row")) gGeoManager->CdUp(); // layered geometry
- if (TString(gGeoManager->GetPath()).Contains("sec")) gGeoManager->CdUp();
- gGeoManager->LocalToMaster(loc,glob);
- phi1 += -TMath::ATan2 (glob[1],glob[0]) * TMath::RadToDeg(); // due to rotation
- if (TMath::Abs(phi1-dPhi*TMath::RadToDeg()/2) > 0.001)
- Fatal("MpdTpcKalmanFilter::FillGeoScheme()"," !!! Inconsistent sens. volume parameters !!! ");
- dZ = ((TGeoPgon*)shape)->GetDZ() * 2;
- //FairGeoNode* membrane = (FairGeoNode*) (passNodes->FindObject("tpc01mb"));
- //if (!membrane) membrane = (FairGeoNode*) (sensNodes->FindObject("tpc01mb")); // it was a test version !
- TGeoVolume *membrane = gGeoManager->GetVolume("tpc01mb");
- dZ += ((TGeoTube*)membrane->GetShape())->GetDZ(); // membrane
- } else Fatal("MpdTpcKalmanFilter::FillGeoScheme()"," !!! Unknown sensitive volume shape !!! ");
- // Use TGeoManager to describe readout chamber geometry
- //TGeoMedium *medium = gGeoManager->GetMedium("TPCmixture");
- //TGeoVolume *tpcVolOrig = gGeoManager->GetVolume("tpc01sv");
- TGeoMedium *medium = gGeoManager->GetMedium("air");
- TGeoVolume *tpcVolOrig = gGeoManager->GetVolume("tpcChamber1");
- TGeoVolume *tpcVol = new TGeoVolume("tpc01sv_new", tpcVolOrig->GetShape(), medium);
- //TGeoVolume* tpcsec[2];
- //Double_t rads[2];
- Double_t dx1 = fSecGeo->GetMinY() * TMath::Tan(dPhi/2);
- Double_t dx2 = fSecGeo->GetMaxY() * TMath::Tan(dPhi/2);
- Double_t dy = dZ, dz = (fSecGeo->GetMaxY() - fSecGeo->GetMinY()) / 2.;
- Double_t rad = (fSecGeo->GetMaxY() + fSecGeo->GetMinY()) / 2.;
- TGeoVolume *tpcsec = gGeoManager->MakeTrd1("tpcsec", medium, dx1, dx2, dy, dz), *vol = NULL;
- //TGeoVolume * Division(const char* name, const char* mother, Int_t iaxis, Int_t ndiv, Double_t start, Double_t step, Int_t numed = 0, Option_t* option = "");
- //gGeoManager->Division("tpcrow", "tpcsec", 3, nRows, 0., 0., 0, "N");
- // The code below is based on TGeoTrd1::Divide()
- Int_t ndiv = nRows[0] + nRows[1];
- Double_t start = -dz, dx1n, dx2n, zmax, step = fSecGeo->PadHeight(), zmin = start - step; // end = dz,
- for (Int_t id = 0; id < ndiv; ++id) {
- zmin += step;
- if (id == nRows[0]) step = fSecGeo->PadHeight(1);
- zmax = zmin + step;
- dx1n = 0.5 * (dx1 * (dz - zmin) + dx2 * (dz + zmin)) / dz;
- dx2n = 0.5 * (dx1 * (dz - zmax) + dx2 * (dz + zmax)) / dz;
- shape = new TGeoTrd1(dx1n, dx2n, dy, step/2.);
- vol = new TGeoVolume("tpcrow", shape, tpcsec->GetMedium());
- TGeoHMatrix *matr = new TGeoHMatrix();
- Double_t transl[3] = {0.0, 0.0, zmin+step/2};
- matr->SetTranslation(transl);
- tpcsec->AddNode(vol, id+1, matr);
- }
- //
- for (Int_t isec = 0; isec < nSec; ++isec) {
- Double_t phi = dPhi * isec;
- TGeoTranslation t(rad*TMath::Cos(phi), rad*TMath::Sin(phi), 0.);
- TGeoRotation r;
- Double_t phiDeg = TMath::RadToDeg() * phi;
- //r.SetAngles(90,phiDeg+90,0,0,90,phiDeg);
- r.SetAngles(90,phiDeg-90,0,0,90,phiDeg);
- TGeoCombiTrans c(t,r);
- TGeoHMatrix *matr = new TGeoHMatrix(c);
- //gGeoManager->RegisterMatrix(matr);
- // AddNode(const TGeoVolume* vol, Int_t copy_no, TGeoMatrix* mat = 0, Option_t* option = "")
- tpcVol->AddNode(tpcsec, isec+1, matr);
- }
- //gGeoManager->ClearPhysicalNodes();
- //gGeoManager->ClearThreadData();
- //gGeoManager->ClearThreadsMap();
- gGeoManager->ReplaceVolume(tpcVolOrig, tpcVol);
- //gGeoManager->RefreshPhysicalNodes();
- // Fill GeoScheme
- TVector3 p3norm, p3pos;
- Double_t padHeight = fSecGeo->PadHeight();
- Double_t rStart = fSecGeo->GetMinY() - padHeight / 2;
- Int_t nRows2 = fSecGeo->NofRows();
- for (Int_t isec = 0; isec < nSec; ++isec) {
- Double_t phi = dPhi * isec;
- p3norm.SetXYZ(TMath::Cos(phi), TMath::Sin(phi), 0.);
- Double_t x0 = rStart * p3norm.X(), y0 = rStart * p3norm.Y();
- padHeight = fSecGeo->PadHeight();
- for (Int_t lay = 0; lay < nRows2; ++lay) {
- Int_t detId = lay * 1000000 + isec;
- geo->SetNormal(detId, p3norm); // normal vector
- if (lay >= nRows[0]) padHeight = fSecGeo->PadHeight(1); // outer ROC region
- x0 += padHeight * p3norm.X();
- y0 += padHeight * p3norm.Y();
- TGeoNode *node = gGeoManager->FindNode(x0, y0, 0.);
- geo->SetPath(detId, gGeoManager->GetPath()); // TGeo path
- geo->SetDetId(gGeoManager->GetPath(), detId); // path to detId
- const Double_t *transl = gGeoManager->GetCurrentMatrix()->GetTranslation();
- p3pos.SetXYZ(transl[0], transl[1], 0.);
- geo->SetGlobalPos(detId, p3pos); // padrow position
- TGeoTrd1 *shape1 = (TGeoTrd1*) node->GetVolume()->GetShape();
- if (TString(shape1->Class()->GetName()) != "TGeoTrd1") { cout <<
- " !!! MpdTpcKalmanFilter::FillGeoScheme(): Wrong shape " << endl; exit(0); }
- p2.Set((shape1->GetDx1()+shape1->GetDx2())/2, shape1->GetDy());
- geo->SetSize(detId, p2); // sector size
- }
- //exit(0);
- }
- }
- //__________________________________________________________________________
- void MpdTpcKalmanFilter::RemoveShorts()
- {
- /// Remove short tracks (Nhits < 4)
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
- Int_t nReco = fTracks->GetEntriesFast();
- for (Int_t i = 0; i < nReco; ++i) {
- MpdKalmanTrack *track = (MpdKalmanTrack*) fTracks->UncheckedAt(i);
- if (track->GetNofHits() < 4) fTracks->RemoveAt(i);
- }
- fTracks->Compress();
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
- }
- //__________________________________________________________________________
- void MpdTpcKalmanFilter::GoToBeamLine()
- {
- /// Propagate tracks to the beam line (mult. scat. and energy loss corrections)
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
- const Int_t nR = 7;
- //const Double_t rad[nR] = {27.035, 27.57, 28.105, 30.64, 33.15, 33.665, 34.178};
- //const Double_t dx[nR] = {0.0031, 0.00085, 0.0031, 0.0003, 0.001, 0.00085, 0.001};
- const Double_t rad[nR] = {27.00, 27.01, 27.31, 27.32, 33.82, 33.83, 34.13};
- const Double_t dx[nR] = {7.8e-4, 1.194e-2, 7.8e-4, 3.3e-4, 7.8e-4, 1.194e-2, 9.6e-4};
- Bool_t ok = 0;
- Int_t nReco = fTracks->GetEntriesFast();
- for (Int_t i = 0; i < nReco; ++i) {
- MpdTpcKalmanTrack *track = (MpdTpcKalmanTrack*) fTracks->UncheckedAt(i);
- if (track->GetNodeNew() == "") continue;
- MpdKalmanHit hit;
- hit.SetDetectorID(-999);
- hit.SetType(MpdKalmanHit::kFixedR);
- for (Int_t j = nR-1; j >= 0; --j) {
- hit.SetPos(rad[j]);
- //MpdKalmanFilter::Instance()->GetGeo()->SetGlobalPos(&hit,TVector3(rad[j],0.,0.),kTRUE); // put hit in GeoScheme
- ok = MpdKalmanFilter::Instance()->PropagateToHit(track,&hit,kTRUE);
- if (!ok) {
- fTracks->RemoveAt(i);
- break;
- }
- //continue;
- // Add multiple scattering
- TMatrixDSym *cov = track->Weight2Cov();
- Double_t th = track->GetParamNew(3);
- Double_t cosTh = TMath::Cos(th);
- Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, dx[j]);
- (*cov)(2,2) += (angle2 / cosTh / cosTh);
- (*cov)(3,3) += angle2;
- //cov->Print();
- Int_t iok = 0;
- MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
- track->SetWeight(*cov);
- }
- if (!ok) continue;
- TMatrixDSym *cov = track->Weight2Cov();
- // Update track
- track->SetParam(*track->GetParamNew());
- track->SetPos(track->GetPosNew());
- // Correct for energy loss
- //*
- Double_t theta = TMath::PiOver2() - track->GetParam(3);
- Double_t pt = 1. / track->GetParam(4);
- Double_t ptCor = CorrectForLoss(TMath::Abs(pt), theta, track->GetTrackID()); // ionization loss correction
- track->SetParam(4,TMath::Sign(1./ptCor,pt));
- track->SetParamNew(*track->GetParam());
- Double_t sigmaPt = CorrectForLossFluct(TMath::Abs(ptCor), theta, 0.1396, 1); // ionization loss correction fluct.
- (*cov)(4,4) += sigmaPt*sigmaPt;
- //*/
-
- //*
- hit.SetPos(0.);
- hit.SetMeas(0,track->GetParam(2));
- //hit.SetRphi(track->GetParam(2)); // track Phi - check if it is necessary !!!!!!!
- //MpdKalmanFilter::Instance()->GetGeo()->SetGlobalPos(&hit,TVector3(0.,0.,0.),kTRUE); // put hit in GeoScheme
- //*/
- /*
- MpdKalmanHitZ hit;
- hit.SetZ(0.);
- */
- //cout << i << " " << track->GetTrackID() << " " << track->GetLength() << " " << ((MpdKalmanHitR*)track->GetHits()->First())->GetLength() << endl;
- Double_t pos = track->GetPosNew();
- //Double_t pos = ((MpdKalmanHitR*)track->GetHits()->Last())->GetR(); // position at last hit
- //MpdKalmanFilter::Instance()->PropagateParamR(track, &hit, kTRUE);
- //track->GetCovariance()->Print();
- //(track->GetWeightAtHit()->Invert()).Print();
- MpdKalmanFilter::Instance()->PropagateToHit(track, &hit, kTRUE);
- //track->Weight2Cov();
- //track->GetCovariance()->Print();
- track->SetPos(pos); // save position after passing inner shell
- track->SetParam(*track->GetParamNew()); // !!! track params at PCA
- //cout << i << " " << track->GetTrackID() << " " << track->GetLength() << " " << ((MpdKalmanHitR*)track->GetHits()->First())->GetLength() << endl;
- track->SetLengAtHit(track->GetLength()-track->GetLengAtHit()); // track length from PCA to the nearest hit
- }
- fTracks->Compress();
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
- }
- //__________________________________________________________________________
- void MpdTpcKalmanFilter::StoreTracks()
- {
- /// Transfer tracks from fTrackCand to fTracks
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
- Int_t nFound = fTracks->GetEntriesFast();
- for (Int_t i = 0; i < fNTracks; ++i) {
- MpdTpcKalmanTrack *track = (MpdTpcKalmanTrack*) fTrackCand->UncheckedAt(i);
- track->Weight2Cov();
- new ((*fTracks)[nFound++]) MpdTpcKalmanTrack(*track);
- //fTrackCand->RemoveAt(i);
- }
- fNTracks = 0;
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
- }
- //__________________________________________________________________________
- void MpdTpcKalmanFilter::AddHits(Int_t indx0)
- {
- /// Add hit objects to tracks and compute number of wrongly assigned hits
- /// (hits with ID different from ID of the track)
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
- Int_t ib = 0, ie = fTracks->GetEntriesFast() - 1;
- if (indx0 >= 0) ib = ie = indx0; // process one track
- for (Int_t i = ib; i <= ie; ++i) {
- MpdTpcKalmanTrack *track = (MpdTpcKalmanTrack*) fTracks->UncheckedAt(i);
- Int_t nHits = track->GetNofHits();
- track->SetNofHits(nHits);
- TClonesArray &trHits = *track->GetTrHits();
- //Int_t idOld = track->GetTrackID(); // 24.12.2018
- SetTrackID(track); // set track ID as ID of majority of its hits
- Int_t nWrong = 0, motherID = track->GetTrackID();
- // Get track mother ID
- MpdMCTrack *mctrack = (MpdMCTrack*) fMCtracks->UncheckedAt(motherID);
- while (mctrack->GetMotherId() >= 0) {
- motherID = mctrack->GetMotherId();
- mctrack = (MpdMCTrack*) fMCtracks->UncheckedAt(mctrack->GetMotherId());
- }
- TObjArray *hits = track->GetHits();
- for (Int_t j = 0; j < nHits; ++j) {
- MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(j);
- MpdKalmanHit *newHit = new (trHits[j]) MpdKalmanHit(*hit);
- MpdTpcHit *p = GetTpcHit(hit);
- Int_t motherID1 = p->GetTrackID();
- //AZ 19-07-15 if (!fUseMCHit) motherID1 = p->GetRefIndex();
- // Get point mother ID
- mctrack = (MpdMCTrack*) fMCtracks->UncheckedAt(motherID1);
- while (mctrack->GetMotherId() >= 0) {
- motherID1 = mctrack->GetMotherId();
- mctrack = (MpdMCTrack*) fMCtracks->UncheckedAt(mctrack->GetMotherId());
- }
- if (motherID1 != motherID) nWrong++;
- hits->AddAt(newHit, j); // replace content - 13.10.2015
- }
- track->SetNofWrong(nWrong);
- track->SetLastLay(((MpdKalmanHit*)track->GetHits()->First())->GetLayer());
- }
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
- }
- //__________________________________________________________________________
- void MpdTpcKalmanFilter::SetTrackID(MpdTpcKalmanTrack* track)
- {
- /// Set track ID as ID of majority of its hits
- const Int_t idMax = 9999;
- vector<Int_t> ids(idMax);
- Int_t nHits = track->GetNofHits(), locMax = 0, size = idMax;
- TObjArray *hits = track->GetHits();
- for (Int_t i = 0; i < nHits; ++i) {
- MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(i);
- Int_t id = GetHitID(hit);
- if (id >= size) { size = id + 1; ids.resize(size); }
- ++ids[id];
- if (ids[id] > ids[locMax]) locMax = id;
- }
- if (ids[locMax] > 1) track->SetTrackID(locMax);
- }
- //__________________________________________________________________________
- void MpdTpcKalmanFilter::GetTrackSeeds(Int_t iPass)
- {
- /// Build track seeds from central detector hits
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
- MpdKalmanGeoScheme *geo = MpdKalmanFilter::Instance()->GetGeo();
- Int_t layMax0 = ((MpdKalmanHit*)fKHits->First())->GetLayer();
- TH1F *hPt = (TH1F*) gROOT->FindObjectAny("hPt");
- TH1F *hZ = (TH1F*) gROOT->FindObjectAny("hZ");
- // Estimate Z-position of vertex
- // Loop over layers
- Int_t layBeg = layMax0, layEnd = layBeg - 2, iDir = -1, dLays = 6; //10;
- //Int_t layMax = layMax0, dLays = 6; //10;
- Double_t zToler = 10;
- if (iPass > 0) {
- //25.12.2018 layBeg = ((MpdKalmanHit*)fKHits->Last())->GetLayer();
- layBeg = 6;
- //25.12.2018 layEnd = layBeg + 2;
- layEnd = layBeg + 4;
- //25.12.2018 iDir = 1;
- iDir = 2;
- //25.12.2018 dLays = 1;
- dLays = 2;
- zToler = 30;
- }
- if (fZflag < 0) zToler = 50; // bad quality of vertex estimate
- MpdTpcKalmanTrack *track = 0x0;
- MpdTpcHit *h1 = 0x0;//, *h2 = 0x0;
- TVector3 pos1, pos2, posLoc;
- Double_t rad1, phi1, rad2, phi2, params[3];
- for (Int_t lay = layBeg; lay != layEnd; lay+=iDir) {
- Int_t nHits1 = (Int_t) fhLays->GetBinContent(lay+1,0), isec = -1;
- Int_t nHits2 = (Int_t) fhLays->GetBinContent(lay+dLays*TMath::Sign(1,iDir)+1,0);
- //cout << "Hits: " << nHits1 << " " << nHits2 << endl;
- // Loop over hits in first layer
- for (Int_t ihit1 = 0; ihit1 < nHits1; ++ihit1) {
- MpdKalmanHit *hit1 = (MpdKalmanHit*) fKHits->UncheckedAt(fLayPointers[lay]+ihit1);
- if (hit1->GetIndex() >= 0) h1 = (MpdTpcHit*) fHits->UncheckedAt(hit1->GetIndex());
- if (hit1->GetFlag() & MpdKalmanHit::kUsed) continue;
- if (fModular) {
- //posLoc.SetXYZ(-hit1->GetMeas(0), hit1->GetDist(), hit1->GetMeas(1));
- posLoc.SetXYZ(hit1->GetMeas(0), hit1->GetDist(), hit1->GetMeas(1));
- isec = hit1->GetDetectorID() % 1000000;
- //pos1 = fPadPlane->fromSectorReferenceFrame(posLoc, isec);
- fSecGeo->Local2Global(isec, posLoc, pos1);
- rad1 = pos1.Pt();
- phi1 = pos1.Phi();
- //cout << posLoc.X() << " " << posLoc.Y() << " " << posLoc.Z() << endl;
- //cout << pos1.X() << " " << pos1.Y() << " " << pos1.Z() << endl;
- //cout << rad1 << " " << phi1 << endl;
- } else {
- rad1 = hit1->GetPos();
- phi1 = hit1->GetMeas(0) / rad1;
- }
- // Loop over hits in second (closer to the beam for the first pass) layer
- for (Int_t ihit2 = 0; ihit2 < nHits2; ++ihit2) {
- MpdKalmanHit *hit2 = (MpdKalmanHit*) fKHits->UncheckedAt(fLayPointers[lay+dLays*TMath::Sign(1,iDir)]+ihit2);
- //if (hit2->GetIndex() >= 0) h2 = (MpdTpcHit*) fHits->UncheckedAt(hit2->GetIndex());
- //if (h2->GetTrackID() != h1->GetTrackID()) continue; //!!! for test - exact ID matching
- if (hit2->GetFlag() & MpdKalmanHit::kUsed) continue;
- if (fModular) {
- //posLoc.SetXYZ(-hit2->GetMeas(0), hit2->GetDist(), hit2->GetMeas(1));
- posLoc.SetXYZ(hit2->GetMeas(0), hit2->GetDist(), hit2->GetMeas(1));
- Int_t isec1 = hit2->GetDetectorID() % 1000000;
- //if (TMath::Abs(isec-isec1) > 1 && TMath::Abs(isec-isec1) < fPadPlane->nSectors()-1) continue;
- //pos2 = fPadPlane->fromSectorReferenceFrame(posLoc, isec1);
- if (TMath::Abs(isec-isec1) > 1 && TMath::Abs(isec-isec1) < fSecGeo->NofSectors()-1) continue;
- fSecGeo->Local2Global(isec1, posLoc, pos2);
- rad2 = pos2.Pt();
- if ((rad2 - rad1) * TMath::Sign(1,iDir) < 0.5) continue;
- phi2 = pos2.Phi();
- } else {
- rad2 = hit2->GetPos();
- phi2 = hit2->GetMeas(0) / rad2;
- }
- Double_t dPhi = MpdKalmanFilter::Instance()->Proxim(phi1,phi2) - phi1;
- if (TMath::Abs(dPhi) > TMath::PiOver4()) continue;
- Double_t pt = EvalPt(hit1, hit2, params);
- //Double_t zvert = hit2->GetMeas(1) - (hit1->GetMeas(1)-hit2->GetMeas(1)) / (rad1-rad2) * rad2;
- Double_t zvert = hit2->GetMeas(1) - (hit1->GetMeas(1) - hit2->GetMeas(1)) / params[0] * params[1];
- if (hZ && iPass) hZ->Fill(zvert);
- if (TMath::Abs(zvert-fVertZ) > zToler) continue;
- //printf("%4d %4d %4d %2d %2d %5.3f %5.3f %7.3f %7.3f %4d %4d %7.3f %6.3f\n", fNTracks, ihit1, ihit2,
- // hit1->GetLayer(), hit2->GetLayer(), rad1, rad2, hit1->GetMeas(1), hit2->GetMeas(1),
- // h1->GetTrackID(), h2->GetTrackID(), zvert, pt);
- if (hPt) hPt->Fill(TMath::Abs(pt));
- //if (TMath::Abs(pt) > 5) cout << " ***!!! " << pt << endl;
- //if (TMath::Abs(pt) < 0.05) continue; // skip low-Pt tracks
- if (TMath::Abs(pt) < 0.02) continue; // skip low-Pt tracks
- if (TMath::Abs(pt) > 10.) continue; // skip high-Pt tracks
- TVector3 vertex(0,0,0);
- if (fModular) {
- if (iPass == 0) track = new ((*fTrackCand)[fNTracks++])
- MpdTpcKalmanTrack(hit1, hit2, vertex, pos1, pos2, pt, params);
- else track = new ((*fTrackCand)[fNTracks++])
- MpdTpcKalmanTrack(hit2, hit1, vertex, pos2, pos1, pt, params);
- TString s1 = geo->Path (hit1->GetDetectorID());
- track->SetNodeNew(s1);
- //cout << hit1->GetDetectorID() << " " << track->GetNodeNew() << endl;
- } else {
- if (iPass == 0) track = new ((*fTrackCand)[fNTracks++]) MpdTpcKalmanTrack(hit1, hit2,
- vertex, pt, params);
- else track = new ((*fTrackCand)[fNTracks++]) MpdTpcKalmanTrack(hit2, hit1, vertex,
- pt, params);
- }
- if (lay <= layMax0 - 1) track->SetDirection(MpdKalmanTrack::kOutward);
- if (iPass > 0) {
- track->GetHits()->Clear();
- track->GetHits()->Add(hit1);
- //track->GetHits()->Add(hit2);
- //if (fModular) track->SetNodeNew(geo->Path(hit2->GetDetectorID())); // 29-apr-12
- // 7-may-12
- track->SetParam(1,hit1->GetMeas(1)); // hit closest to the beam line
- track->SetPos(TMath::Min(rad1,rad2));
- //31.12.2018 track->SetParam(0,track->GetParam(0)/track->GetPosNew()*track->GetPos());
- //31.12.2018
- if (rad1 < rad2) track->SetParam(0,MpdKalmanFilter::Instance()->Proxim(track->GetParam(0)/rad2,pos1.Phi())*rad1);
- else track->SetParam(0,MpdKalmanFilter::Instance()->Proxim(track->GetParam(0)/rad1,pos2.Phi())*rad2);
- //
- track->SetPosNew(track->GetPos());
- track->SetParamNew(*track->GetParam());
- }
- //*
- if (fModular) {
- /*TString s1, s2;
- s1 += hit1->GetDetectorID() % 1000000; // sector No.
- s2 += hit2->GetDetectorID() % 1000000; // sector No.
- if (iPass > 0) { track->SetNode(s1); track->SetNodeNew(s2); }
- else { track->SetNode(s2); track->SetNodeNew(s1); }
- */
- /*
- TString s1 = geo->Path(hit1->GetDetectorID());
- TString s2 = geo->Path(hit2->GetDetectorID());
- if (iPass > 0) track->SetNodeNew(s2); // trick
- else track->SetNodeNew(s1);
- */
- }
- //*/
- if (h1) track->SetTrackID(h1->GetTrackID());
- //((MpdKalmanTrack*)fTrackCand->UncheckedAt(fNTracks-1))->GetParam()->Print();
- //cout << fNTracks-1 << endl;
- //((TpcLheKalmanTrack*)fTrackCand->UncheckedAt(fNTracks-1))->Weight2Cov()->Print();
- }
- }
- }
- //hLays->Write();
- //hZ->Write();
- //hPt->Write();
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
- }
- //__________________________________________________________________________
- void MpdTpcKalmanFilter::GetTrackSeedsEndCaps()
- {
- /// Build track seeds from end-cap regions (inward tracks from hits near endcaps)
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
- MpdKalmanGeoScheme *geo = MpdKalmanFilter::Instance()->GetGeo();
- Int_t layBeg = fSecGeo->NofRows() - 3, layEnd = 4, iDir = -1, dLays = 5;
- Double_t zToler = 10, zMax = 160;
- if (fZflag < 0) zToler = 50; // bad quality of vertex estimate
- MpdTpcKalmanTrack *track = 0x0;
- MpdTpcHit *h1 = 0x0;//, *h2 = 0x0;
- TVector3 pos1, pos2, posLoc;
- Double_t rad1, phi1, rad2, phi2, params[3];
- // Loop over layers
- for (Int_t lay = layBeg; lay != layEnd; lay+=iDir) {
- Int_t nHits1 = (Int_t) fhLays->GetBinContent(lay+1,0), isec = -1;
- Int_t nHits2 = (Int_t) fhLays->GetBinContent(lay+dLays*iDir+1,0);
- if (nHits1 <= 0 || nHits2 <= 0) continue;
- // Loop over hits in first layer
- for (Int_t ihit1 = 0; ihit1 < nHits1; ++ihit1) {
- MpdKalmanHit *hit1 = (MpdKalmanHit*) fKHits->UncheckedAt(fLayPointers[lay]+ihit1);
- if (hit1->GetIndex() >= 0) h1 = (MpdTpcHit*) fHits->UncheckedAt(hit1->GetIndex());
- if (hit1->GetFlag() & MpdKalmanHit::kUsed) continue;
- if (TMath::Abs(hit1->GetMeas(1)) > zMax) continue; // too close to end-cap
- if (fModular) {
- //posLoc.SetXYZ(-hit1->GetMeas(0), hit1->GetDist(), hit1->GetMeas(1));
- posLoc.SetXYZ(hit1->GetMeas(0), hit1->GetDist(), hit1->GetMeas(1));
- isec = hit1->GetDetectorID() % 1000000;
- fSecGeo->Local2Global(isec, posLoc, pos1);
- rad1 = pos1.Pt();
- phi1 = pos1.Phi();
- Double_t tanThe = rad1 / TMath::Abs(pos1.Z()-fVertZ);
- Int_t row = fSecGeo->PadRow(hit1->GetDetectorID());
- Double_t padHeight = fSecGeo->PadHeight();
- if (row >= fSecGeo->NofRowsReg(0)) padHeight = fSecGeo->PadHeight(1);
- if (TMath::Abs(pos1.Z()) < zMax - padHeight / tanThe) continue; // too far from end-cap
- //cout << posLoc.X() << " " << posLoc.Y() << " " << posLoc.Z() << endl;
- //cout << pos1.X() << " " << pos1.Y() << " " << pos1.Z() << endl;
- //cout << rad1 << " " << phi1 << endl;
- } else {
- rad1 = hit1->GetPos();
- phi1 = hit1->GetMeas(0) / rad1;
- }
- // Loop over hits in second (closer to the beam for the first pass) layer
- for (Int_t ihit2 = 0; ihit2 < nHits2; ++ihit2) {
- MpdKalmanHit *hit2 = (MpdKalmanHit*) fKHits->UncheckedAt(fLayPointers[lay+dLays*iDir]+ihit2);
- //if (hit2->GetIndex() >= 0) h2 = (MpdTpcHit*) fHits->UncheckedAt(hit2->GetIndex());
- //if (h2->GetTrackID() != h1->GetTrackID()) continue; //!!! for test - exact ID matching
- if (hit2->GetFlag() & MpdKalmanHit::kUsed) continue;
- if (fModular) {
- //posLoc.SetXYZ(-hit2->GetMeas(0), hit2->GetDist(), hit2->GetMeas(1));
- posLoc.SetXYZ(hit2->GetMeas(0), hit2->GetDist(), hit2->GetMeas(1));
- Int_t isec1 = hit2->GetDetectorID() % 1000000;
- //if (TMath::Abs(isec-isec1) > 1 && TMath::Abs(isec-isec1) < fPadPlane->nSectors()-1) continue;
- if (TMath::Abs(isec-isec1) > 1 && TMath::Abs(isec-isec1) < fSecGeo->NofSectors()-1) continue;
- fSecGeo->Local2Global(isec1, posLoc, pos2);
- rad2 = pos2.Pt();
- if ((rad2 - rad1) * iDir < 0.5) continue;
- phi2 = pos2.Phi();
- } else {
- rad2 = hit2->GetPos();
- phi2 = hit2->GetMeas(0) / rad2;
- }
- Double_t dPhi = MpdKalmanFilter::Instance()->Proxim(phi1,phi2) - phi1;
- if (TMath::Abs(dPhi) > TMath::PiOver4()) continue;
- Double_t pt = EvalPt(hit1, hit2, params);
- // Double_t zvert = hit2->GetMeas(1) - (hit1->GetMeas(1)-hit2->GetMeas(1)) / (rad1-rad2) * rad2;
- Double_t zvert = hit2->GetMeas(1) - (hit1->GetMeas(1) - hit2->GetMeas(1)) / params[0] * params[1];
- if (TMath::Abs(zvert-fVertZ) > zToler) continue;
- //cout << fNTracks << " " << ihit1 << " " << ihit2 << " " << hit1->GetLayer() << " " << hit2->GetLayer() << " " << rad1 << " " << rad2 << " " << hit1->GetMeas(1) << " " << hit2->GetMeas(1) << " " << h1->GetTrackID() << " " << h2->GetTrackID() << " " << zvert << " " /*<< hit1->GetUsage() << " " << hit2->GetUsage() << " "*/ << pt << endl;
- if (TMath::Abs(pt) < 0.02) continue; // skip low-Pt tracks
- if (TMath::Abs(pt) > 10.) continue; // skip high-Pt tracks
- TVector3 vertex(0,0,0);
- if (fModular) {
- track = new ((*fTrackCand)[fNTracks++]) MpdTpcKalmanTrack(hit1, hit2, vertex,
- pos1, pos2, pt, params);
- TString s1 = geo->Path (hit1->GetDetectorID());
- track->SetNodeNew(s1);
- //cout << hit1->GetDetectorID() << " " << track->GetNodeNew() << endl;
- } else {
- track = new ((*fTrackCand)[fNTracks++]) MpdTpcKalmanTrack(hit1, hit2, vertex,
- pt, params);
- }
- track->SetDirection(MpdKalmanTrack::kOutward);
- //*
- if (fModular) {
- /*TString s1, s2;
- s1 += hit1->GetDetectorID() % 1000000; // sector No.
- s2 += hit2->GetDetectorID() % 1000000; // sector No.
- if (iPass > 0) { track->SetNode(s1); track->SetNodeNew(s2); }
- else { track->SetNode(s2); track->SetNodeNew(s1); }
- */
- /*
- TString s1 = geo->Path(hit1->GetDetectorID());
- TString s2 = geo->Path(hit2->GetDetectorID());
- if (iPass > 0) track->SetNodeNew(s2); // trick
- else track->SetNodeNew(s1);
- */
- }
- //*/
- if (h1) track->SetTrackID(h1->GetTrackID());
- }
- }
- }
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
- }
- //__________________________________________________________________________
- void MpdTpcKalmanFilter::Cluster2KalmanHits()
- {
- /// Create Kalman hits from Rec. Points found from clusters
- Int_t lay, layMax = 0, nClus = fHits->GetEntriesFast();
- // !!!!!!!!! Geometry dependent values
- //Double_t rMin = 20.5, rMax = 109.5, dR = (rMax-rMin)/50; // dR = 1.78; // old version
- //Double_t rMin = 29.195, rMax = 99.81, dR = (rMax-rMin)/50; // 1.4123; // new version (with dead material)
- Double_t rMin = 34.195, rMax = 99.81, dR = (rMax-rMin)/50; // 1.1962; // new version (with dead material)
- //Double_t rMin = 35.00, rMax = 100.00, dR = (rMax-rMin)/50; // 1.1962; // new version (with dead material)
- // !!!!!!!!!
- Double_t xyz[3], rPhiErr = 0.05, zErr = 0.1;
- //TClonesArray hits("TpcLheHit");
- for (int j = 0; j < nClus; ++j ) {
- TpcCluster* clus = (TpcCluster*) fHits->UncheckedAt(j);
- xyz[0] = clus->ret_x();
- xyz[1] = clus->ret_y();
- xyz[2] = clus->ret_z();
- //TpcLheHit *hit = new (hits[j]) TpcLheHit(xyz, 0);
- //hit->SetR (TMath::Sqrt (xyz[0]*xyz[0]+xyz[1]*xyz[1]));
- Double_t rad = TMath::Sqrt (xyz[0] * xyz[0] + xyz[1] * xyz[1]);
- Double_t rPhi = TMath::ATan2(xyz[1],xyz[0]) * rad;
- lay = (Int_t) ((rad - rMin) / dR);
- lay = TMath::Max (lay, 0);
- layMax = TMath::Max (lay, layMax);
- //hit->SetLayer(lay);
- fhLays->Fill(lay+0.1);
- // Create Kalman hits
- //MpdKalmanHit *hitK = new ((*fKHits)[j]) MpdKalmanHit(rad, rPhi, xyz[2],
- // rPhiErr, zErr, lay, -1);
- //(Int_t detID, Int_t nDim, HitType hitType, Double_t *meas, Double_t *err, Double_t *cosSin, Double_t signal, Double_t dist, Int_t index)
- Double_t meas[2] = {rPhi, xyz[2]};
- Double_t err[2] = {rPhiErr, zErr};
- Double_t cossin[2] = {1., 0.};
- MpdKalmanHit *hitK = new ((*fKHits)[j]) MpdKalmanHit(0, 2, MpdKalmanHit::kFixedR, meas, err, cossin, 0., rad, -1);
- }
- cout << " Max layer = " << layMax << " " << fKHits->GetEntriesFast() << endl;
- fKHits->Sort();
- fLayPointers = new Int_t [layMax+1];
- Int_t ipos = 0;
- for (Int_t i = layMax; i >= 0; --i) {
- //cout << i << " " << fhLays->GetBinContent(i+1,0) << endl;
- //if (ipos) cout << ((TpcLheHit*)fHits->UncheckedAt(ipos))->GetLayer() << " "
- // << ((TpcLheHit*)fHits->UncheckedAt(ipos-1))->GetLayer() << endl;
- fLayPointers[i] = ipos;
- ipos += (Int_t) fhLays->GetBinContent(i+1,0);
- }
- }
- //__________________________________________________________________________
- void MpdTpcKalmanFilter::MakeKalmanHits()
- {
- /// Create Kalman hits.
- Int_t nHits = fHits->GetEntriesFast();
- cout << " Number of merged hits: " << nHits << endl;
- Int_t nKh = 0, lay = 0, layMax = 0;
- for (int j = 0; j < nHits; ++j ) {
- MpdTpcHit* hit = (MpdTpcHit*) fHits->UncheckedAt(j);
- //if (hit->GetTrackID() != 1783) continue; // keep only one track
- lay = hit->GetLayer();
- layMax = TMath::Max (lay, layMax);
- fhLays->Fill(lay+0.1);
- // Add errors
- Double_t dRPhi = 0, dZ = 0;
- gRandom->Rannor(dRPhi,dZ);
- Double_t z = hit->GetZ() + dZ * hit->GetDz(); //add error
- Double_t rPhi = hit->GetRphi() + dRPhi * hit->GetDx(); //add error
- //Double_t phi = rPhi / hit->GetR();
- //hit->SetX(hit->GetR()*TMath::Cos(phi));
- //hit->SetY(hit->GetR()*TMath::Sin(phi));
- // Create Kalman hits
- //(Int_t detID, Int_t nDim, HitType hitType, Double_t *meas, Double_t *err, Double_t *cosSin, Double_t signal, Double_t dist, Int_t index)
- Double_t meas[2] = {rPhi, z};
- Double_t err[2] = {hit->GetDx(), hit->GetDz()};
- Double_t cossin[2] = {1., 0.};
- MpdKalmanHit *hitK = new ((*fKHits)[nKh++])
- MpdKalmanHit(lay*1000000+nKh-1, 2, MpdKalmanHit::kFixedR, meas, err, cossin, hit->GetEnergyLoss()/hit->GetStep(), hit->GetR(), j);
- hitK->SetLength(hit->GetLength());
- //hitK->SetDedx (hit->GetEdep()/hit->GetStep());
- //MpdKalmanFilter::Instance()->GetGeo()->SetGlobalPos(hitK,TVector3(hit->GetX(),hit->GetY(),hit->GetZ()));
- }
- fLayPointers = new Int_t [layMax+1];
- Int_t ipos = 0;
- for (Int_t i = layMax; i >= 0; --i) {
- //cout << i << " " << fhLays->GetBinContent(i+1,0) << endl;
- //if (ipos) cout << ((TpcLheHit*)fHits->UncheckedAt(ipos))->GetLayer() << " "
- // << ((TpcLheHit*)fHits->UncheckedAt(ipos-1))->GetLayer() << endl;
- fLayPointers[i] = ipos;
- fLaySecBegPointers[i][0] = ipos;
- ipos += (Int_t) fhLays->GetBinContent(i+1,0);
- fLaySecEndPointers[i][0] = ipos;
- }
- }
- //__________________________________________________________________________
- void MpdTpcKalmanFilter::MakeKalmanHitsModul()
- {
- /// Create Kalman hits for the modular geometry of readout chambers.
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
- Int_t nHits = fHits->GetEntriesFast();
- cout << " Number of merged hits: " << nHits << endl;
- Int_t nKh = 0, lay = 0, layMax = 0;
- for (int j = 0; j < nHits; ++j ) {
- MpdTpcHit* hit = (MpdTpcHit*) fHits->UncheckedAt(j);
- lay = hit->GetLayer();
- layMax = TMath::Max (lay, layMax);
- fhLays->Fill(lay+0.1);
- // Add errors
- Double_t dX = 0, dZ = 0;
- gRandom->Rannor(dX, dZ);
- // Create Kalman hits
- //(Int_t detID, Int_t nDim, HitType hitType, Double_t *meas, Double_t *err, Double_t *cosSin, Double_t signal, Double_t dist, Int_t index)
- Double_t err[2] = {hit->GetDx(), hit->GetDz()};
- //Double_t meas[2] = {-hit->GetLocalX()-dX*err[0], hit->GetLocalZ()+dZ*err[1]};
- Double_t meas[2] = {hit->GetLocalX(), hit->GetLocalZ()};
- if (fUseMCHit) { meas[0] += dX * err[0]; meas[1] += dZ * err[1]; }
- Double_t cossin[2] = {1., 0.};
- Int_t padID = hit->GetDetectorID();
- //TpcPadID padIDobj = TpcPadID::numberToPadID(padID);
- //padID = padIDobj.sector();
- padID = fSecGeo->Sector(padID);
- Double_t xsec = (hit->GetLocalY() + fSecGeo->GetMinY()) * TMath::Tan(fSecGeo->Dphi()/2);
- Double_t xloc = hit->GetLocalX(); Double_t edge = 0.0; // distance to sector boundary
- if (xloc < 0) edge = xloc + xsec;
- else edge = xloc - xsec;
- MpdKalmanHit *hitK = new ((*fKHits)[nKh++])
- MpdKalmanHit(lay*1000000+padID, 2, MpdKalmanHit::kFixedP, meas, err, cossin, hit->GetEnergyLoss()/hit->GetStep(), hit->GetLocalY(), j, edge);
- hitK->SetLength(hit->GetLength());
- hitK->SetFlag (hitK->GetFlag() | hit->GetFlag());
- //hitK->SetDedx (hit->GetEdep()/hit->GetStep());
- //MpdKalmanFilter::Instance()->GetGeo()->SetGlobalPos(hitK,TVector3(hit->GetX(),hit->GetY(),hit->GetZ()));
- }
- fKHits->Sort();
- fLayPointers = new Int_t [layMax+1];
- Int_t ipos = 0;
- for (Int_t i = layMax; i >= 0; --i) {
- //cout << i << " " << fhLays->GetBinContent(i+1,0) << endl;
- //if (ipos) cout << ((TpcLheHit*)fHits->UncheckedAt(ipos))->GetLayer() << " "
- // << ((TpcLheHit*)fHits->UncheckedAt(ipos-1))->GetLayer() << endl;
- fLayPointers[i] = ipos;
- ipos += (Int_t) fhLays->GetBinContent(i+1,0);
- }
- // Fill indices for different sectors
- /*
- Int_t iseco = -1, isec = 0;
- for (Int_t i = layMax; i >= 0; --i) {
- Int_t nLay = GetNofHitsInLayer(lay);
- Int_t indx0 = GetHitsInLayer(lay);
- ipos = indx0;
- for (Int_t indx = 0; indx != nLay; ++indx) {
- MpdKalmanHit *hit = (MpdKalmanHit*) fKHits->UncheckedAt(indx0+indx);
- isec = hit->GetDetectorID() % 1000000;
- if (isec != iseco) {
- fLaySecBegPointers[i].insert(pair<Int_t,Int_t>(isec,ipos+indx));
- if (iseco > 0) fLaySecEndPointers[i].insert(pair<Int_t,Int_t>(iseco,ipos+indx));
- iseco = isec;
- }
- }
- }
- fLaySecEndPointers[0].insert(pair<Int_t,Int_t>(isec,nHits));
- */
- nHits = fKHits->GetEntriesFast();
- for (Int_t j = 0; j < nHits; ++j ) {
- MpdKalmanHit *hit = (MpdKalmanHit*) fKHits->UncheckedAt(j);
- lay = hit->GetLayer();
- Int_t isec = hit->GetDetectorID() % 1000000;
- if (fLaySecBegPointers[lay][isec] < 0) {
- fLaySecBegPointers[lay][isec] = j;
- fLaySecEndPointers[lay][isec] = j + 1;
- } else {
- if (fLaySecBegPointers[lay][isec] > j) fLaySecBegPointers[lay][isec] = j;
- if (fLaySecEndPointers[lay][isec] < j+1) fLaySecEndPointers[lay][isec] = j + 1;
- }
- }
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
- }
- //__________________________________________________________________________
- Double_t MpdTpcKalmanFilter::EvalPt(const MpdKalmanHit *hit1, const MpdKalmanHit *hit2, Double_t *params)
- {
- /// Evaluate signed track Pt (curvature) assuming the track coming from the
- /// primary vertex
- TVector3 pos, posLoc;
- Double_t rad1, phi1, rad2, phi2;
- if (fModular) {
- posLoc.SetXYZ(-hit1->GetMeas(0), hit1->GetDist(), 0.);
- Int_t isec = hit1->GetDetectorID() % 1000000;
- //pos = fPadPlane->fromSectorReferenceFrame(posLoc, isec);
- fSecGeo->Local2Global(isec, posLoc, pos);
- rad1 = pos.Pt();
- phi1 = pos.Phi();
- posLoc.SetXYZ(-hit2->GetMeas(0), hit2->GetDist(), 0.);
- isec = hit2->GetDetectorID() % 1000000;
- //pos = fPadPlane->fromSectorReferenceFrame(posLoc, isec);
- fSecGeo->Local2Global(isec, posLoc, pos);
- rad2 = pos.Pt();
- phi2 = pos.Phi();
- } else {
- rad1 = hit1->GetPos();
- rad2 = hit2->GetPos();
- phi1 = hit1->GetMeas(0) / rad1;
- phi2 = hit2->GetMeas(0) / rad2;
- }
- TVector2 vec1(rad1*TMath::Cos(phi1)-0.,rad1*TMath::Sin(phi1)-0.);
- TVector2 vec2(rad2*TMath::Cos(phi2)-0.,rad2*TMath::Sin(phi2)-0.);
- TVector2 vec21 = vec1 - vec2;
- Double_t dist2 = vec2.Mod();
- Double_t dist12 = vec21.Mod();
- Double_t cosAlpha = vec2 * vec21 / dist2 / dist12;
- Double_t rad = vec1.Mod() / 2. / TMath::Sin(TMath::ACos(cosAlpha));
- params[0] = TMath::ASin (dist12 / 2 / rad) * rad * TMath::Sign(2.0,rad1-rad2); // arc12 length
- params[1] = TMath::ASin (dist2 / 2 / rad) * 2 * rad; // arc02 length
- params[2] = params[0] / rad; // arc12 angle
- Double_t bz = FairRunAna::Instance()->GetField()->GetBz(0.,0.,0.);
- Double_t factor = 0.0003 * bz; // 0.3 * 0.01 * 5kG / 10
- Double_t charge = phi1 - MpdKalmanFilter::Instance()->Proxim(phi1,phi2);
- if (hit1->GetLayer() < hit2->GetLayer()) charge = -charge;
- return factor * TMath::Abs(rad) * TMath::Sign(1., -charge);
- }
- //__________________________________________________________________________
- void MpdTpcKalmanFilter::DoTracking(Int_t iPass)
- {
- /// Run Kalman tracking
- //struct timeval tvStart, tvEnd;
- //struct timezone tz;
- //gettimeofday(&tvStart, &tz);
- //fTracks->Sort(); // in descending order in Pt
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
- Int_t nTracks = fTrackCand->GetEntriesFast();
- //TH1F *hChi2 = (TH1F*) gROOT->FindObjectAny("hChi2");
- //TH1F *hNhits = (TH1F*) gROOT->FindObjectAny("hNhits");
- //#pragma omp parallel for
- for (Int_t i = 0; i < nTracks; ++i) {
- Int_t iok = 0;
- MpdTpcKalmanTrack *track = (MpdTpcKalmanTrack*) fTrackCand->UncheckedAt(i);
- MpdKalmanTrack::TrackDir dir0 = MpdKalmanTrack::kInward;
- if (((MpdKalmanHit*)track->GetHits()->First())->GetLayer() < 10) dir0 = MpdKalmanTrack::kOutward;
- //cout << " Track seed No. " << i << " " << track->GetTrackID() << " " << ((MpdKalmanHit*)track->GetHits()->First())->GetLayer() << endl;
- // << track->GetPosNew() << endl;
- //track->GetParamNew()->Print();
- // Check if the seed hit was used already
- MpdKalmanHit *hit = (MpdKalmanHit*) track->GetHits()->First();
- if (hit->GetFlag() & MpdKalmanHit::kUsed) iok = -1;
- else {
- if (track->GetDirection() == MpdKalmanTrack::kOutward) {
- //cout << " !!! " << track->GetNode() << " !!! " << track->GetNodeNew() <<endl;
- GoOutward(track);
- /*cout << i << " " << track->GetChi2() << endl;*/
- if (track->GetChi2()>1000) {
- //cout << " !!! " << " " << track->GetNofHits() << " " << track->GetNode() << " " << track->GetNodeNew() << endl;
- }
- }
- iok = RunKalmanFilter(track);
- //if (fNofEvents > 54)
- //cout << iok << " " << track->GetChi2() << " " << track->GetNofHits() << endl;
- }
- //if (iOK != 0) track->SetStatus(-1);
- if (iok != 0) {
- //fTracks->RemoveAt(i);
- //cout << track->GetNofHits() << endl;
- //#pragma omp critical
- {
- fTrackCand->RemoveAt(i);
- }
- //delete track; // !!! Object has been deleted already
- } else {
- track->SetDirection(dir0); // restore initial direction flag
- /*
- Int_t nHits = track->GetNofHits();
- //#pragma omp critical
- {
- hChi2->Fill(track->GetChi2());
- }
- //#pragma omp critical
- {
- hNhits->Fill(nHits+0.1);
- }
- */
- //if (track->GetChi2() < 1.) cout << " Strange: " << track->GetChi2() << " " << nHits << " " << track->GetDirection() << endl;
- // Mark hits as being used
- /*
- TObjArray *hits = track->GetHits();
- for (Int_t j = 0; j < nHits; ++j) {
- hit = (MpdKalmanHit*) hits->UncheckedAt(j);
- hit->SetFlag(-1);
- }
- */
- }
- if (iPass == 0 && track->GetNofHits() < 10) fTrackCand->RemoveAt(i); //AZ - 190820
- }
- //fTracks->Compress();
- fTrackCand->Compress();
- /*
- nTracks = fTrackCand->GetEntriesFast();
- for (Int_t i = 0; i < nTracks; ++i) {
- TpcLheKalmanTrack *track = (TpcLheKalmanTrack*) fTrackCand->UncheckedAt(i);
- Int_t nHits = track->GetNofHits(), iCharge = 0;
- TObjArray *hits = track->GetHits();
- TVector3 pmom, pmom1;
- for (Int_t j = 0; j < nHits; ++j) {
- TpcLheHit *hit = (TpcLheHit*) hits->UncheckedAt(j);
- if (j == nHits-1) {
- // Compare with MC track
- Int_t nMC = fLHEtracks->GetEntriesFast();
- for (Int_t itr = 0; itr < nMC; ++itr) {
- TpcLheTrack *trMC = (TpcLheTrack*) fLHEtracks->UncheckedAt(itr);
- if (trMC->GetTrackNumber() != hit->GetTrackID()) continue;
- pmom = trMC->GetMomentum();
- break;
- }
- TpcPoint *point = (TpcPoint*) fTpcPoints->UncheckedAt(hit->GetUniqueID());
- point->Momentum(pmom1);
- iCharge = -hit->GetRefIndex();
- }
- printf(" %4d",hit->GetTrackID());
- }
- printf(" hits\n");
- cout << " Momentum: " << i << " " << 1./track->GetParamNew(4) << " " << pmom1.Pt()*iCharge << " " << pmom.Pt() << " " << track->GetChi2() << endl;
- } // for (Int_t i = 0; i < nTracks;
- */
- // Remove doubles
- RemoveDoubles(fTrackCand);
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
- //gettimeofday(&tvEnd, &tz);
- //int dif = 1000000*((int)tvEnd.tv_sec - (int)tvStart.tv_sec) + (tvEnd.tv_usec - tvStart.tv_usec);
- //cout<<"TPC::DoTracking "<<dif<<" ns\n";
- }
- //__________________________________________________________________________
- Int_t MpdTpcKalmanFilter::RunKalmanFilter(MpdTpcKalmanTrack *track)
- {
- /// Run Kalman filter
- //cout << fHits->GetEntriesFast() << endl;
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
- const Double_t dStepMax = 9; // max track length distance between consecutive hits
- MpdKalmanGeoScheme *geo = MpdKalmanFilter::Instance()->GetGeo();
- Int_t layMax = ((MpdKalmanHit*)fKHits->First())->GetLayer();
- MpdKalmanHit *hitOK = (MpdKalmanHit*) track->GetHits()->Last();
- Int_t layOK = hitOK->GetLayer();
- MpdKalmanTrack::TrackDir trackDir = track->GetDirection();
- //Int_t layBeg = layOK - 1, layEnd = -1, dLay = -1;
- Int_t layBeg = layOK - 1, layEnd = ((MpdKalmanHit*)fKHits->Last())->GetLayer() - 1, dLay = -1;
- if (trackDir == MpdKalmanTrack::kOutward) {
- layBeg = layOK + 1;
- layEnd = layMax + 1;
- dLay = 1;
- }
-
- //Int_t indxOK = hits->IndexOf(hitOK);
- //Int_t nHits = hits->GetEntriesFast();
- Int_t miss = 0, iOK = 1;
- if (fModular) iOK = 0;
- TMatrixDSym pointWeight(5), pointWeightTmp(5), saveWeight(5), weightOK = *track->GetWeight();
- TMatrixD param(5,1), paramTmp(5,1), paramOK = *track->GetParamNew();
- Double_t saveR = 0., rOK = track->GetPosNew(), saveLeng = 0., lengOK = 0.;
- TString saveNode, nodeOK = track->GetNodeNew();
- //cout << " Starting hit: " << hitOK->GetLayer() << " " << hitOK->GetTrackID() << " " << hitOK->GetUsage() << endl;
-
- // Loop over layers
- Int_t nLay;
- for (Int_t lay = layBeg; lay != layEnd; lay+=dLay) {
- //#pragma omp critical
- {
- //nLay = GetNofHitsInLayer(lay);
- }
- //Int_t indx0 = GetHitsInLayer(lay);
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),"Time1");
- Double_t dChi2Min = 1.e+6, padH = lay < MpdTpcSectorGeo::Instance()->NofRowsReg(0) ?
- MpdTpcSectorGeo::Instance()->PadHeight() : MpdTpcSectorGeo::Instance()->PadHeight(1);
- MpdKalmanHit *hitMin = 0x0;
- //cout << " lay, nLay: " << lay << " " << nLay << " " << indx0 << endl;
- Int_t indxBeg = 0, indxEnd = nLay, dIndx = 1, secFirst = -1, isecHit = -1, isecDif = 0;
- if (trackDir == MpdKalmanTrack::kOutward) {
- indxBeg = nLay - 1;
- indxEnd = -1;
- dIndx = -1;
- }
- Int_t isec = 0;
- if (fModular) isec = SectorNo(track->GetNodeNew()); // sector number
- Int_t sector3[3] = {-1, -1, -1};
- MpdTpcKalmanTrack* trackBr[3] = {NULL, NULL, NULL};
- // Loop over hits of 3 sectors (central (where track is located) and +- 1)
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),"Hits");
- Int_t i3sec[3] = { 0, -1, 1};
- for (Int_t iadj = 0; iadj < 3; ++iadj) {
- Int_t isecLoop = isec + i3sec[iadj];
- //if (isecLoop < 0) isecLoop += fPadPlane->nSectors();
- //else if (isecLoop > fPadPlane->nSectors() - 1) isecLoop -= fPadPlane->nSectors();
- if (isecLoop < 0) isecLoop += fSecGeo->NofSectors();
- else if (isecLoop > fSecGeo->NofSectors() - 1) isecLoop -= fSecGeo->NofSectors();
- if (fLaySecBegPointers[lay][isecLoop] < 0) continue;
- indxBeg = fLaySecBegPointers[lay][isecLoop];
- indxEnd = fLaySecEndPointers[lay][isecLoop];
- if (trackDir == MpdKalmanTrack::kOutward) {
- Int_t tmp = indxBeg;
- indxBeg = indxEnd - 1;
- indxEnd = tmp - 1;
- }
- MpdKalmanHit *hit = NULL;
- for (Int_t indx = indxBeg; indx != indxEnd; indx+=dIndx) {
- //MpdKalmanHit *hit = (MpdKalmanHit*) fKHits->UncheckedAt(indx0+indx);
- hit = (MpdKalmanHit*) fKHits->UncheckedAt(indx);
- // Exclude used hits
- if (hit->GetFlag() & MpdKalmanHit::kUsed) continue;
- // !!! Exact ID match
- //if (GetPoint(hit)->GetTrackID() != track->GetTrackID()) continue;
- //if (GetTpcHit(hit)->GetTrackID() != track->GetTrackID()) continue;
- // Check if the hit within some window (15x9cm for the moment - check!!!)
- if (!fModular) {
- if (trackBr[0] == NULL) trackBr[0] = new MpdTpcKalmanTrack(*track);
- if (TMath::Abs(hit->GetMeas(1)-trackBr[0]->GetParamNew(1)) > 9) continue;
- if (TMath::Abs(MpdKalmanFilter::Instance()->Proxim(trackBr[0],hit)
- -trackBr[0]->GetParamNew(0)) > 15) continue;
- } else {
- isecHit = hit->GetDetectorID() % 1000000;
- isecDif = isecHit - isec;
- // Track and hit not in adjacent sectors
- //if (TMath::Abs(isecDif) > 1 && TMath::Abs(isecDif) < fPadPlane->nSectors()-1) break;
- //if (TMath::Abs(isecDif) > 1) isecDif -= TMath::Sign(fPadPlane->nSectors(),isecDif);
- if (TMath::Abs(isecDif) > 1 && TMath::Abs(isecDif) < fSecGeo->NofSectors()-1) break;
- if (TMath::Abs(isecDif) > 1) isecDif -= TMath::Sign(fSecGeo->NofSectors(),isecDif);
- ++isecDif;
- if (isecHit == secFirst && sector3[isecDif] < 0) break; // excluded sector
- if (isecHit != secFirst && trackBr[isecDif] == NULL) trackBr[isecDif] = new MpdTpcKalmanTrack(*track);
- if (TMath::Abs(hit->GetMeas(1)-trackBr[isecDif]->GetParamNew(1)) > 9) continue;
- }
- Double_t leng = trackBr[isecDif]->GetLength();
- if (fModular && isecHit != secFirst) {
- // Modular geometry
- secFirst = isecHit;
- sector3[isecDif] = 1;
- if (trackBr[isecDif]->GetNode() == "") trackBr[isecDif]->SetNodeNew("");
- //if (!MpdKalmanFilter::Instance()->PropagateToHit(&tr1,hit,kTRUE,kTRUE) ||
- // TMath::Abs(leng - tr1.GetLength()) < MpdKalmanFilter::Instance()->fgkEpsilon) {
- if (!MpdKalmanFilter::Instance()->PropagateToHit(trackBr[isecDif],hit,kTRUE,kTRUE)) {
- // Not reaching padrow
- sector3[isecDif] = -1;
- //iOK = 0;
- break; //continue;
- }
- //tr1.GetParamNew()->Print();
- //cout << " New node: " << tr1.GetNodeNew() << endl;
- iOK = 1;
- // Check for edge effect
- if (TMath::Abs(trackBr[isecDif]->GetParamNew(0)) > geo->Size(hit).X() + 0.5) {
- // Outside target sector
- sector3[isecDif] = -1;
- break; //continue;
- }
- trackBr[isecDif]->SetNode(trackBr[isecDif]->GetNodeNew());
- } else if (fModular == 0) {
- if (!MpdKalmanFilter::Instance()->PropagateToHit(trackBr[0],hit)) { iOK = 0; break; }
- sector3[0] = 1;
- } // if (fModular && isecHit != secFirst)
- // Add multiple scattering in TPC gas (TPCmixture)
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),"Scat");
- Double_t step = trackBr[isecDif]->GetLength() - leng;
- Double_t th = trackBr[isecDif]->GetParamNew(3);
- Double_t cosTh = TMath::Cos(th);
- if (1 && step > 1.e-4) {
- Double_t x0 = 13363.6; // rad. length - TPCMixture
- TMatrixDSym *cov = trackBr[isecDif]->Weight2Cov();
- Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(trackBr[isecDif], x0, step);
- //cout << " Scat: " << hit->GetLayer() << " " << step << " " << TMath::Sqrt(angle2) << endl;
- (*cov)(2,2) += (angle2 / cosTh / cosTh);
- (*cov)(3,3) += angle2;
- Int_t iok = 0;
- MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
- trackBr[isecDif]->SetWeight(*cov);
- }
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),"Scat");
- Double_t dChi2 = 999999.0;
- if (trackBr[isecDif]->GetLength() - track->GetLength() < padH / TMath::Abs(cosTh) * 1.8 * (miss + 1))
- dChi2 = MpdKalmanFilter::Instance()->FilterHit(trackBr[isecDif],hit,pointWeight,param);
- //cout << lay << " " << isecHit << " " << isec << " " << dChi2 << " " << track->GetTrackID()
- // << " " << GetTpcHit(hit)->GetTrackID() << endl;
- if (dChi2 < dChi2Min) {
- dChi2Min = dChi2;
- hitMin = hit;
- saveWeight = *trackBr[isecDif]->GetWeight();
- saveR = trackBr[isecDif]->GetPosNew();
- saveLeng = trackBr[isecDif]->GetLength();
- saveNode = trackBr[isecDif]->GetNodeNew();
- // temporary storage for the current track
- paramTmp = param;
- pointWeightTmp = pointWeight;
- /*
- cout << " New min dChi2 = " << dChi2 << " " << hitMin->GetRphi() << " " << hitMin->GetZ() << endl;
- cout << track->GetParamNew(0) << " " << track->GetParamNew(1) << endl;
- cout << hit->GetRphi() << " " << hit->GetZ() << endl;
- cout << param(0,0) << " " << param(1,0) << endl;
- */
- //paramTmp.Print();
- }
- } // for (Int_t indx = indxBeg; indx != indxEnd;
- // No hits in the sector where the track is located (or its padrow can't be reached)- try both neighbours
- //26-dec-2012 if (fModular && (trackBr[1] == NULL || sector3[1] < 0)) continue;
- if (fModular && trackBr[1] == NULL) continue;
- if (fModular && sector3[1] < 0) {
- if (trackBr[isecDif]->GetParam(0) > 0) ++iadj; // skip one step in the loop
- continue;
- }
- if (!fModular || isecHit < 0 || TMath::Abs(trackBr[isecDif]->GetParamNew(0)) < geo->Size(hit).X() - 1.0) break;
- // Close to sector edge - try adjacent sector
- //i3sec = (trackBr[isecDif]->GetParamNew(0) < 0) ? -1 : 1;
- //if (-trackBr[isecDif]->GetParamNew(0) > 0) ++iadj; // skip one step in the loop
- if (trackBr[isecDif]->GetParamNew(0) > 0) ++iadj; // skip one step in the loop
- } //for (Int_t iadj = 0; iadj < 3; ++iadj)
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),"Time1");
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),"Hits");
- if (!iOK) {
- // Curling track
- //*
- track->SetPos(rOK); //
- track->SetPosNew(rOK); //
- track->SetParam(paramOK); //
- track->SetParamNew(paramOK); //
- track->SetWeight(weightOK); //
- track->SetLength(lengOK); //
- track->SetNode(nodeOK);
- track->SetNodeNew(nodeOK);
- //*/
- for (Int_t jj = 0; jj < 3; ++jj) if (trackBr[jj]) delete trackBr[jj];
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
- return -1;
- } else if (dChi2Min < fgkChi2Cut) {
- //if (fNofEvents > 54)
- //cout << lay << " " << hitMin->GetDetectorID() % 1000000 << " " << isec << " " << dChi2Min << " "
- // << track->GetTrackID() << " " << GetTpcHit(hitMin)->GetTrackID() << " " << saveLeng-track->GetLength() << endl;
- //layOK = lay;
- hitOK = hitMin;
- track->GetHits()->Add(hitOK);
- miss = 0;
- // Restore track params at the best hit
- track->SetChi2(track->GetChi2()+dChi2Min);
- //if (track->GetTrackID() >= 0) cout << " *** Adding hit: " << track->GetTrackID() << " " << GetHitID(hitOK) << " " << hitOK->GetLayer() << " " << hitOK->GetDetectorID() << " "/*<< hitOK->GetTrackID() << " " << hitOK->GetUsage() << " "*/ << dChi2Min << " " << track->GetChi2() << " " << track->GetLength() << " " << saveLeng-track->GetLength() << endl;
- saveWeight += pointWeightTmp;
- track->SetWeight(saveWeight);
- track->SetPos(saveR);
- track->SetPosNew(saveR);
- track->SetLength(saveLeng);
- track->SetParam(paramTmp);
- track->SetParamNew(paramTmp);
- track->SetNode(saveNode);
- track->SetNodeNew(saveNode);
- if (track->GetDirection() == MpdKalmanTrack::kInward) {
- // Save track params at last hit
- track->SetLengAtHit(saveLeng);
- track->SetParamAtHit(paramTmp);
- track->SetWeightAtHit(saveWeight);
- } else {
- // 14-may-12
- track->GetSteps().insert(pair<Int_t,Double_t>(lay,saveLeng));
- }
- // Check if the accepted hit is the same as the seeded hit
- //if (hitOK->GetLayer() == f2ndHit->GetLayer() && hitOK != f2ndHit) return -1; // abandon track
- paramOK = paramTmp;
- rOK = saveR;
- lengOK = saveLeng;
- weightOK = saveWeight;
- nodeOK = saveNode;
- } else {
- ++miss;
- Int_t ibr = -1;
- for (Int_t jj = 0; jj < 3; ++jj) {
- if (sector3[jj] > 0) { ibr = jj; break; }
- }
- if (ibr < 0 || TMath::Abs(trackBr[ibr]->GetLength() - track->GetLength()) > dStepMax) {
- //cout << " !!! Too few branches " << sector3.size() << endl;
- //exit(0);
- track->SetPos(rOK); //
- track->SetPosNew(rOK); //
- track->SetParam(paramOK); //
- track->SetParamNew(paramOK); //
- track->SetWeight(weightOK); //
- track->SetLength(lengOK); //
- //27-04-2012 track->SetNode(nodeOK);
- track->SetNodeNew(nodeOK);
- } else {
- Int_t ip = ibr;
- Double_t leng = 999999.;
- for (Int_t jj = ip ; jj < 3; ++jj) {
- if (sector3[jj] < 0) continue;
- if (trackBr[jj]->GetLength() < leng) { ip = jj; leng = trackBr[ip]->GetLength(); }
- }
- // 12-may-12
- Int_t secPos1 = trackBr[ip]->GetNodeNew().Index("tpcsec");
- Int_t pos1 = trackBr[ip]->GetNodeNew().Index("/",secPos1);
- Int_t secPos2 = track->GetNodeNew().Index("tpcsec");
- Int_t pos2 = track->GetNodeNew().Index("/",secPos2);
- if (trackBr[ip]->GetNodeNew()(secPos1,pos1-secPos1) != track->GetNodeNew()(secPos2,pos2-secPos2)) {
- track->SetPos(rOK); //
- track->SetPosNew(rOK); //
- track->SetParam(paramOK); //
- track->SetParamNew(paramOK); //
- track->SetWeight(weightOK); //
- track->SetLength(lengOK); //
- track->SetNodeNew(nodeOK);
- } else {
- track->SetPos(trackBr[ip]->GetPosNew()); //
- track->SetPosNew(trackBr[ip]->GetPosNew()); //
- track->SetParam(*trackBr[ip]->GetParamNew()); //
- track->SetParamNew(*trackBr[ip]->GetParamNew()); //
- track->SetWeight(*trackBr[ip]->GetWeight()); //
- track->SetLength(trackBr[ip]->GetLength()); //
- track->SetNode(trackBr[ip]->GetNodeNew());
- track->SetNodeNew(trackBr[ip]->GetNodeNew());
- }
- }
- //if (miss > 1) {
- if (miss > 5) {
- /*
- track->SetPos(rOK); //
- track->SetPosNew(rOK); //
- track->SetParam(paramOK); //
- track->SetParamNew(paramOK); //
- */
- //return -1;
- track->SetPos(rOK); //
- track->SetPosNew(rOK); //
- track->SetParam(paramOK); //
- track->SetParamNew(paramOK); //
- track->SetWeight(weightOK); //
- track->SetLength(lengOK); //
- track->SetNodeNew(nodeOK);
- for (Int_t jj = 0; jj < 3; ++jj) if (trackBr[jj]) delete trackBr[jj];
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
- return 0; // keep short tracks
- }
- }
- //cout << " lay, miss: " << lay << " " << miss << " " << dChi2Min << " " << fChi2 << endl;
- for (Int_t jj = 0; jj < 3; ++jj) if (trackBr[jj]) delete trackBr[jj];
- } // for (Int_t lay = layBeg; lay != layEnd; lay+=dLay) {
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
- return 0;
- }
- //__________________________________________________________________________
- Int_t MpdTpcKalmanFilter::SectorNo(const char* cpath)
- {
- /// Extract sector number from TGeo path
- TString path(cpath);
- Int_t ipos = path.Index("tpcrow_");
- if (ipos < 0) return -1;
- TString sub = TString(path(ipos-3, 2)); // check for 2 digits
- if (!sub.IsDigit()) sub = TString(path(ipos-2, 1)); // 1 digit
- return sub.Atoi() - 1;
- }
- //__________________________________________________________________________
- MpdTpcHit* MpdTpcKalmanFilter::GetTpcHit(MpdKalmanHit *hit)
- {
- /// Get MpdTpcHit pointer for the Kalman hit
- return (MpdTpcHit*) fHits->UncheckedAt(hit->GetIndex());
- }
- //__________________________________________________________________________
- Int_t MpdTpcKalmanFilter::GetHitID(MpdKalmanHit *hit)
- {
- /// Get hit ID from MpdTpcHit ID
- MpdTpcHit *h = GetTpcHit(hit);
- //AZ 19-07-15 return (fUseMCHit) ? h->GetTrackID() : h->GetRefIndex();
- return h->GetTrackID();
- }
- //__________________________________________________________________________
- void MpdTpcKalmanFilter::GoOutward(MpdTpcKalmanTrack *track)
- {
- /// Propagate track in the outward direction
- TMatrixDSym weight = *track->GetWeight(); // save starting weight matrix
- Int_t ok = RunKalmanFilter(track);
- Int_t lastLay = ((MpdKalmanHit*)track->GetHits()->Last())->GetLayer();
- //cout << ok << " " << track->GetChi2() << " " << track->GetNofHits() << " " << lastLay << " " << ((MpdKalmanHit*)track->GetHits()->Last())->GetMeas(1) << endl;
- track->SetLastLay(lastLay);
- track->SetLength(0.);
- BackTrace(track, weight);
- //cout << track->GetChi2() << " " << track->GetTrackID() << " " << track->GetNofHits() << endl;
- }
- //__________________________________________________________________________
- Bool_t MpdTpcKalmanFilter::BackTrace(MpdTpcKalmanTrack *track, TMatrixDSym &weight0, Int_t iDir, Bool_t corDedx)
- {
- /// Propagate track using already found hits
- track->StartBack();
- if (iDir > 0) track->SetDirection(MpdKalmanTrack::kInward);
- else track->SetDirection(MpdKalmanTrack::kOutward);
- Int_t nHits = track->GetNofHits();
- //28-apr-2012 if (nHits == 1) track->SetWeight(weight0); // restore original weight
- //if (fModular && track->GetNode() == "") track->SetNodeNew(""); // 27-04-2012
- TMatrixD param(5,1);
- TMatrixDSym weight(5), pointWeight(5);
- Int_t ibeg = iDir > 0 ? 0 : nHits-1;
- Int_t iend = iDir > 0 ? nHits-1 : 0;
- Int_t lay0 = -1, lay = -1, ifirst = 1;
- iend += iDir;
- Bool_t ok = kTRUE;
- MpdKalmanHit *hit = 0x0;
- map<Int_t,Double_t>& stepMap = track->GetSteps();
- for (Int_t i = ibeg; i != iend; i+=iDir) {
- //for (Int_t i = 1; i < nHits; ++i) {
- if (i == ibeg && iDir > 0) track->SetLength(0.); // for correct track length
- hit = (MpdKalmanHit*) track->GetHits()->UncheckedAt(i);
- lay0 = lay;
- lay = hit->GetLayer();
- /*AZ-100620 - W.t.f. is this???
- if (!fUseMCHit && corDedx && i == ibeg) {
- Double_t cosTh = TMath::Cos(track->GetParam(3));
- hit->SetSignal(hit->GetSignal() * cosTh);
- }
- */
- if (i == ibeg) continue; // skip first hit
- Double_t leng = track->GetLength(), stepBack = -1;
- if (stepMap.find(lay) != stepMap.end() && stepMap.find(lay0) != stepMap.end())
- stepBack = TMath::Abs (stepMap[lay] - stepMap[lay0]);
- ok = MpdKalmanFilter::Instance()->PropagateToHit(track, hit, kTRUE, kTRUE, stepBack);
- if (!ok) continue;
- //cout << track->GetLength() << endl;
- // Add multiple scattering in TPC gas (TPCmixture)
- Double_t step = track->GetLength() - leng;
- if (1 && step > 1.e-4) {
- Double_t x0 = 13363.6; // rad. length
- TMatrixDSym *cov = track->Weight2Cov();
- Double_t th = track->GetParamNew(3);
- Double_t cosTh = TMath::Cos(th);
- Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, x0, step);
- //cout << " Scat: " << hit->GetLayer() << " " << step << " " << TMath::Sqrt(angle2) << endl;
- (*cov)(2,2) += (angle2 / cosTh / cosTh);
- (*cov)(3,3) += angle2;
- Int_t iok = 0;
- MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
- track->SetWeight(*cov);
- }
- Double_t dChi2 = MpdKalmanFilter::Instance()->FilterHit(track,hit,pointWeight,param);
- if (ifirst && dChi2 > 1000) {
- ifirst = 0;
- //cout << " !!! Back: " << dChi2 << " " << nHits << " " << hit->GetLayer() << " " << ok << " " << step << " " << track->GetNode() << " " << track->GetNodeNew() << endl;
- }
- track->SetChi2(track->GetChi2()+dChi2);
- weight = *track->GetWeight();
- weight += pointWeight;
- track->SetWeight(weight);
- track->SetParamNew(param);
- //if (i == 0 && iDir > 0) track->SetLength(0.); // for correct track length
- // Correct dE/dx for step length (for clusters) - due to track angles
- /* Angle corrections moved to Dedx task
- if (!fUseMCHit && corDedx) {
- Double_t cosTh = TMath::Cos(track->GetParamNew(3));
- Double_t phi = track->GetParamNew(2);
- //Double_t secPhi = fSecGeo->SectorAngle(fSecGeo->Sector(hit->GetDetectorID()));
- Double_t secPhi = phi; // !!! it seems to work better !!!
- Double_t cosPhi = TMath::Cos(phi-secPhi);
- if (cosPhi < 0.1) { cout << " oops " << secPhi << " " << phi << " " << track->Pt() << endl; cosPhi = 1; }
- hit->SetSignal(hit->GetSignal() * cosTh * cosPhi);
- }
- */
- }
- // Save track params at last hit
- track->SetLengAtHit(track->GetLength());
- track->SetParamAtHit(*track->GetParamNew());
- track->SetWeightAtHit(*track->GetWeight());
- track->GetHits()->UnSort(); // 20-02-2012
- track->GetHits()->Sort(); // 20-02-2012
- return ok;
- }
- //__________________________________________________________________________
- void MpdTpcKalmanFilter::GoOut()
- {
- /// Backpropagate tracks in the outward direction for matching with outer detectors
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Start(Class()->GetName(),__FUNCTION__);
- Int_t nReco = fTracks->GetEntriesFast(), isec = 0;
- TMatrixDSym tmp;
- TVector3 posLoc, pos;
- for (Int_t i = 0; i < nReco; ++i) {
- MpdTpcKalmanTrack *track = (MpdTpcKalmanTrack*) fTracks->UncheckedAt(i);
- //track->GetSteps().clear();
- MpdTpcKalmanTrack tr = *track;
- //Int_t id = tr.GetTrackID();
- tr.SetParam(*tr.GetParamAtHit());
- tr.SetParamNew(*tr.GetParamAtHit());
- tr.SetWeight(*tr.GetWeightAtHit());
- tr.SetLength(tr.GetLengAtHit());
- MpdKalmanHit *hit = (MpdKalmanHit*) tr.GetTrHits()->Last();
- tr.SetPos(hit->GetPos());
- if (fModular) {
- tr.SetNode(MpdKalmanFilter::Instance()->GetGeo()->Path(hit->GetDetectorID()));
- tr.SetNodeNew(tr.GetNode());
- //posLoc.SetXYZ(-hit->GetMeas(0), hit->GetDist(), hit->GetMeas(1));
- posLoc.SetXYZ(hit->GetMeas(0), hit->GetDist(), hit->GetMeas(1));
- isec = hit->GetDetectorID() % 1000000;
- //pos = fPadPlane->fromSectorReferenceFrame(posLoc, isec);
- fSecGeo->Local2Global(isec, posLoc, pos);
- gGeoManager->cd(tr.GetNodeNew());
- gGeoManager->CdUp();
- Double_t v7[3], v77[3];
- pos.GetXYZ(v7);
- gGeoManager->MasterToLocal(v7,v77);
- tr.SetPos(v77[2]);
- }
- tr.SetPosNew(tr.GetPos());
- //cout << tr.GetTrackID() << " " << tr.GetLengAtHit() << " " << tr.GetPos() << endl;
- Bool_t ok = BackTrace(&tr,tmp,-1,kTRUE);
- //cout << " aaaaa " << tr.Momentum3().Eta() << " " << tr.GetLength() << endl;
- if (ok) {
- MpdKalmanHit hitTmp;
- hitTmp.SetType(MpdKalmanHit::kFixedR);
- if (fModular) {
- // Transform to global frame
- /*
- hit = (MpdKalmanHit*) tr.GetTrHits()->First();
- posLoc.SetXYZ(-hit->GetMeas(0), hit->GetDist(), hit->GetMeas(1));
- isec = hit->GetDetectorID() % 1000000;
- pos = fPadPlane->fromSectorReferenceFrame(posLoc, isec);
- */
- //Double_t v7[3] = {tr.GetParamNew(0), tr.GetParamNew(1), tr.GetPosNew()}, v77[3];
- Double_t v7[3] = {-tr.GetParamNew(0), tr.GetParamNew(1), tr.GetPosNew()}, v77[3];
- gGeoManager->cd(tr.GetNodeNew());
- gGeoManager->CdUp();
- gGeoManager->LocalToMaster(v7,v77);
- //MpdKalmanHit hitTmp;
- hitTmp.SetDist(TMath::Sqrt(v77[0]*v77[0]+v77[1]*v77[1]) + 0.2);
- //hitTmp.SetType(MpdKalmanHit::kFixedR);
- } else hitTmp.SetDist(tr.GetPosNew() + 0.2);
- ok = MpdKalmanFilter::Instance()->PropagateToHit(&tr,&hitTmp);
- }
- if (!ok) track->SetFlag(MpdKalmanTrack::kNoOut); //AZ
- track->SetLengAtHit(tr.GetLength());
- track->SetParamAtHit(*tr.GetParamNew());
- track->SetWeightAtHit(*tr.GetWeight());
- track->SetPosAtHit(tr.GetPosNew());
- }
- // Below is just for debugging
- if (0) {
- TClonesArray *tofPoints = (TClonesArray*) FairRootManager::Instance()->GetObject("TOFPoint");
- Int_t nTof = tofPoints->GetEntriesFast();
- //cout << nTof << endl;
- for (Int_t i = 0; i < nReco; ++i) {
- MpdTpcKalmanTrack *track = (MpdTpcKalmanTrack*) fTracks->UncheckedAt(i);
- MpdTpcKalmanTrack tr = *track;
- //cout << i << " " << tr.GetTrackID() << endl;
- tr.SetParam(*tr.GetParamAtHit());
- tr.SetParamNew(*tr.GetParamAtHit());
- tr.SetWeight(*tr.GetWeightAtHit());
- tr.SetLength(tr.GetLengAtHit());
- tr.SetPos(((MpdKalmanHit*)tr.GetTrHits()->First())->GetPos());
- tr.SetPosNew(tr.GetPos());
- for (Int_t j = 0; j < nTof; ++j) {
- FairMCPoint *p = (FairMCPoint*) tofPoints->UncheckedAt(j);
- if (p->GetTrackID() != tr.GetTrackID()) continue;
- p->Position(pos);
- MpdKalmanHit hitTmp;
- hitTmp.SetDist(pos.Pt());
- hitTmp.SetType(MpdKalmanHit::kFixedR);
- MpdKalmanFilter::Instance()->PropagateToHit(&tr,&hitTmp);
- Double_t dl = tr.GetLength() - p->GetLength();
- //cout << dl << " " << tr.GetPosNew()-pos.Pt() << " " << tr.GetParamNew(1)-pos.Z() << " " << pos.Z() << endl;
- if (lunTpc) fprintf(lunTpc, "%10.3f %10.3f %10.3f %10.3f \n",
- dl, tr.GetPosNew()-pos.Pt(), tr.GetParamNew(1)-pos.Z(), pos.Z());
- break;
- }
- }
- }
- if (MpdCodeTimer::Active()) MpdCodeTimer::Instance()->Stop(Class()->GetName(),__FUNCTION__);
- }
- //__________________________________________________________________________
- void MpdTpcKalmanFilter::RemoveDoubles(TClonesArray *tracks)
- {
- /// Remove double tracks (keep the ones with better quality)
- Int_t ntracks = tracks->GetEntriesFast();
- cout << " Total tracks: " << ntracks << endl;
- // Fill map with track quality
- multimap<Double_t,MpdTpcKalmanTrack*> qualMap;
- multimap<Double_t,MpdTpcKalmanTrack*>::iterator mit1, mit2;
- for (Int_t i = 0; i < ntracks; i++) {
- MpdTpcKalmanTrack *tr = (MpdTpcKalmanTrack*) tracks->UncheckedAt(i);
- //AZ-190820 Double_t quality = tr->GetNofHits() + TMath::Min(tr->GetChi2(),999.0) / 1000.0;
- Double_t quality = tr->GetNofHits() + (1 - TMath::Min(tr->GetChi2(),999.0) / 1000.0);
- qualMap.insert(pair<Double_t,MpdTpcKalmanTrack*>(-quality,tr));
- }
- MpdTpcKalmanTrack *tr1, *tr2;
- for (mit1 = qualMap.begin(); mit1 != qualMap.end(); ++mit1) {
- tr1 = mit1->second;
- if (tr1 == NULL) continue; // killed track
- mit2 = mit1;
- ++mit2;
- for ( ; mit2 != qualMap.end(); ++mit2) {
- tr2 = mit2->second;
- if (tr2 == NULL) continue; // killed track
- //Int_t nHitsCommon = GetNofCommonHits(tr1, tr2);
- //if ((float)nHitsCommon / TMath::Min(tr1->GetNofHits(),tr2->GetNofHits()) < 0.5) continue;
- if (!AreTracksDoubles(tr1, tr2)) continue;
- if (mit1->first < mit2->first) { tracks->Remove(tr2); mit2->second = NULL; }
- else {
- if (tr1->GetChi2() <= tr2->GetChi2()) { tracks->Remove(tr2); mit2->second = NULL; }
- else { tracks->Remove(tr1); mit1->second = NULL; break; }
- }
- }
- }
- tracks->Compress();
- fNTracks = tracks->GetEntriesFast();
- }
- //__________________________________________________________________________
- Int_t MpdTpcKalmanFilter::GetNofCommonHits(MpdKalmanTrack *tr1, MpdKalmanTrack *tr2)
- {
- /// Compute number of common hits in 2 tracks
- TObjArray *hits1 = tr1->GetHits(), *hits2 = tr2->GetHits();
- Int_t nh1 = hits1->GetEntriesFast(), nh2 = hits2->GetEntriesFast();
- Int_t nCom = 0;//, id = -1;
- for (Int_t i = 0; i < nh1; ++i) {
- MpdKalmanHit *hit1 = (MpdKalmanHit*) hits1->UncheckedAt(i);
- for (Int_t j = 0; j < nh2; ++j) {
- MpdKalmanHit *hit2 = (MpdKalmanHit*) hits2->UncheckedAt(j);
- if (hit1 == hit2) {
- ++nCom;
- //id = hit1->GetTrackID();
- break;
- }
- if (hit2->GetLayer() < hit1->GetLayer()) break; // already closer to beam
- }
- }
- /*
- if (nCom > 0) cout << " nCom: " << nCom << " " << nh1 << " " << nh2 << " " << id
- << " " << 1./tr1->GetParamNew(4) << " " << 1./tr2->GetParamNew(4)
- << " " << tr1->GetChi2() << " " << tr2->GetChi2() << endl;
- */
- return nCom;
- }
- //__________________________________________________________________________
- Bool_t MpdTpcKalmanFilter::AreTracksDoubles(MpdKalmanTrack *tr1, MpdKalmanTrack *tr2)
- {
- /// Searching common hits in 2 tracks to determine doubles
- //track1 consists of less count of hits than track2
- MpdKalmanTrack *track1 = tr1, *track2 = tr2;
- if (tr1->GetNofHits() > tr2->GetNofHits())
- track1 = tr2, track2 = tr1;
- Int_t limCommonPoint = (track1->GetNofHits()+1) / 2; // at least so many common hits should be found
- TObjArray *hits1 = track1->GetHits(), *hits2 = track2->GetHits();
- Int_t nh1 = hits1->GetEntriesFast(), nh2 = hits2->GetEntriesFast(), nHitsCommon = 0, j = 0;
- for (Int_t i = 0; i < nh1; i++){
- MpdKalmanHit *hit1 = (MpdKalmanHit*) hits1->UncheckedAt(i);
- for ( ; j < nh2; j++){
- MpdKalmanHit *hit2 = (MpdKalmanHit*) hits2->UncheckedAt(j);
- // Is the hit common for two tracks?
- if (hit1 == hit2) {
- nHitsCommon++;
- if (nHitsCommon == limCommonPoint) return kTRUE; // already enough common hits
- break;
- }
- if (hit2->GetLayer() < hit1->GetLayer()) break; // already closer to beam
- }
- if (i+limCommonPoint-nHitsCommon > nh1) return kFALSE; // there'll be not enough common hits already
- }
- if (nHitsCommon < limCommonPoint) return kFALSE; // not too many common hits
- return kTRUE;
- }
- //__________________________________________________________________________
- void MpdTpcKalmanFilter::ExcludeHits()
- {
- /// Exclude hits, already used for tracking, from consideration during the next passes
- Int_t nReco = fTracks->GetEntriesFast();
- cout << " nReco: " << nReco << endl;
- for (Int_t i = 0; i < nReco; ++i) {
- MpdKalmanTrack *track = (MpdKalmanTrack*) fTracks->UncheckedAt(i);
- Int_t nhitsKF = track->GetNofHits();
- TObjArray *hits = track->GetHits();
- for (Int_t j = 0; j < nhitsKF; ++j) {
- MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(j);
- hit->SetFlag (hit->GetFlag() | MpdKalmanHit::kUsed);
- }
- }
- }
- //__________________________________________________________________________
- Double_t MpdTpcKalmanFilter::CorrectForLoss(Double_t pt, Double_t the, Int_t id)
- {
- /// Apply momentum correction due to ionization loss in the beam pipe (1 mm of Be)
- /// and TPC inner shell (tpc_v6) - for pions and electrons
- /// Attention!!! Neither pion nor muon with p = 0.050 GeV/c survive - they decay
- /// before entering TPC
- const Double_t piMass = 0.13957, eMass = 0.00051;
- Double_t mass;
- Int_t elec = GetParticleId(id);
- if (elec) mass = eMass;
- else mass = piMass;
- return CorrectForLoss(TMath::Abs(pt), the, mass, 1); // ionization loss correction
- }
- //__________________________________________________________________________
- /*
- Double_t MpdTpcKalmanFilter::CorrectForLoss(Double_t pt, Double_t the, Int_t id)
- {
- /// Apply momentum correction due to ionization loss in the beam pipe (1 mm of Be)
- /// and TPC inner shell (tpc_v6) - for pions and electrons
- /// Attention!!! Neither pion nor muon with p = 0.050 GeV/c survive - they decay
- /// before entering TPC
- const Int_t nDim = 13;
- const Double_t mom[nDim]= {0.075, 0.100, 0.200, 0.300, 0.400, 0.500, 0.700, 0.900,
- 1.200, 1.500, 2.000, 2.500, 3.500};
- const Double_t tkin[nDim]= {0.019, 0.032, 0.104, 0.191, 0.284, 0.380, 0.574, 0.771,
- 1.069, 1.367, 1.865, 2.364, 3.363};
- const Double_t dedxp[nDim]={8.700, 4.900, 2.350, 1.950, 1.850, 1.850, 1.750, 1.750,
- 1.750, 1.750, 1.850, 1.850, 1.850}; // peak
- const Double_t dedxm[nDim]={8.773, 5.002, 2.597, 2.259, 2.175, 2.143, 2.129, 2.131,
- 2.152, 2.154, 2.182, 2.193, 2.239}; // mean of fit Gauss*Landau
- const Double_t piMass = 0.13957, eMass = 0.00051;
- Double_t dt, mass;
- Int_t elec = GetParticleId(id);
- if (elec) mass = eMass;
- else mass = piMass;
- Double_t mass2 = mass * mass;
- Double_t p = pt / TMath::Sin(the);
- Double_t t = TMath::Sqrt (p*p + mass2) - mass;
- if (elec) dt = 2.9; // 2.9 MeV loss for electrons
- else {
- if (t < tkin[0]) dt = dedxm[0];
- //if (t < tkin[0]) {
- //dt = dedxm[0] + (dedxm[1]-dedxm[0])/(tkin[1]-tkin[0]) * (t-tkin[0]);
- //dt = TMath::Min (dt,10.);
- //}
- else if (t > tkin[nDim-1]) dt = dedxm[nDim-1];
- else {
- // Binary search
- Int_t i1 = 0, i2 = nDim-1, i = i2;
- do {
- i = i1 + (i2-i1) / 2;
- if (t > tkin[i]) i1 = i;
- else i2 = i;
- //cout << i << " " << i1 << " " << i2 << " " << t << " " << tkin[i1] << " " << tkin[i2] << endl;
- } while (i2 - i1 > 1);
- // Linear interpolation
- dt = dedxm[i1] + (dedxm[i2]-dedxm[i1])/(tkin[i2]-tkin[i1]) * (t-tkin[i1]);
- }
- }
- dt /= TMath::Sin(the);
- t += dt * 1.e-3;
- Double_t e = mass + t;
- p = TMath::Sqrt (e*e - mass2);
- pt = p * TMath::Sin(the);
- return pt;
- }
- */
- /*
- //__________________________________________________________________________
- Double_t MpdTpcKalmanFilter::CorrectForLoss(Double_t pt, Double_t the, Int_t id)
- {
- /// Apply momentum correction due to ionization loss in the beam pipe (0.5 mm of Al)
- /// and TPC inner shell
- //return pt; ///
- const Int_t nDim = 13;
- const Double_t mom[nDim]={0.075,0.100,0.200,0.300,0.400,0.500,0.700,0.900,
- 1.200,1.500,2.000,2.500,3.500};
- const Double_t tkin[nDim]={0.019,0.032,0.104,0.191,0.284,0.380,0.574,0.771,
- 1.069,1.367,1.865,2.364,3.363};
- const Double_t dedxp[nDim]={3.125,1.950,0.975,0.800,0.800,0.750,0.750,0.750,
- 0.750,0.750,0.750,0.775,0.775}; // peak
- const Double_t dedxm[nDim]={3.195,2.085,1.100,0.952,0.918,0.898,0.898,0.890,
- 0.887,0.886,0.901,0.904,0.913}; // mean
- const Double_t piMass = 0.13957, eMass = 0.00051;
- Double_t dt, mass;
-
- Int_t elec = GetParticleId(id);
- if (elec) mass = eMass;
- else mass = piMass;
- Double_t mass2 = mass * mass;
- Double_t p = pt / TMath::Sin(the);
- Double_t t = TMath::Sqrt (p*p + mass2) - mass;
- if (elec) dt = 1.45; // 1.45 MeV loss for electrons
- else {
- //if (t < tkin[0]) dt = dedxm[0];
- if (t < tkin[0]) {
- dt = dedxm[0] + (dedxm[1]-dedxm[0])/(tkin[1]-tkin[0]) * (t-tkin[0]);
- dt = TMath::Min (dt,10.);
- }
- else if (t > tkin[nDim-1]) dt = dedxm[nDim-1];
- else {
- // Binary search
- Int_t i1 = 0, i2 = nDim-1, i = i2;
- do {
- i = i1 + (i2-i1) / 2;
- if (t > tkin[i]) i1 = i;
- else i2 = i;
- //cout << i << " " << i1 << " " << i2 << " " << t << " " << tkin[i1] << " " << tkin[i2] << endl;
- } while (i2 - i1 > 1);
- // Linear interpolation
- dt = dedxm[i1] + (dedxm[i2]-dedxm[i1])/(tkin[i2]-tkin[i1]) * (t-tkin[i1]);
- }
- }
- dt /= TMath::Sin(the);
- t += dt * 1.e-3;
- Double_t e = mass + t;
- p = TMath::Sqrt (e*e - mass2);
- pt = p * TMath::Sin(the);
- return pt;
- }
- */
- //__________________________________________________________________________
- Int_t MpdTpcKalmanFilter::GetParticleId(Int_t id)
- {
- /// Particle ID:
- /// !!! based on MC information at the moment !!!
- MpdMCTrack* mcTr = (MpdMCTrack*) fMCtracks->UncheckedAt(id);
- Int_t pdg = mcTr->GetPdgCode();
- if (TMath::Abs(pdg) == 11) return 1;
- return 0;
- }
- //__________________________________________________________________________
- //Bool_t MpdTpcKalmanFilter::SameOrigin(TpcLheHit *hit, Int_t idKF, Int_t *mcTracks)
- Bool_t MpdTpcKalmanFilter::SameOrigin(TpcPoint *hit, Int_t idKF, Int_t *mcTracks)
- {
- /// Check if two hits belong to tracks from the same mother
- //cout << " Wrongs: " << hit->GetTrackID() << " " << idKF << " " << mcTracks[hit->GetTrackID()] << " " << mcTracks[idKF] << endl;
- if (mcTracks[hit->GetTrackID()] <= 0 && mcTracks[idKF] <= 0) return kFALSE; // both hits are from primaries
- else if (mcTracks[hit->GetTrackID()] > 0) {
- // Check mother
- MpdMCTrack *track = (MpdMCTrack*) fMCtracks->UncheckedAt(hit->GetTrackID());
- while (track->GetMotherId() > 0) {
- //cout << " Wrongs1: " << track->GetMotherId() << endl;
- if (track->GetMotherId() == idKF) return kTRUE;
- track = (MpdMCTrack*) fMCtracks->UncheckedAt(track->GetMotherId());
- }
- } else {
- MpdMCTrack *track = (MpdMCTrack*) fMCtracks->UncheckedAt(idKF);
- while (track->GetMotherId() > 0) {
- //cout << " Wrongs2: " << track->GetMotherId() << endl;
- if (track->GetMotherId() == hit->GetTrackID()) return kTRUE;
- track = (MpdMCTrack*) fMCtracks->UncheckedAt(track->GetMotherId());
- }
- }
- return kFALSE;
- }
- //__________________________________________________________________________
- Bool_t MpdTpcKalmanFilter::Refit(MpdKalmanTrack *track, Double_t mass, Int_t charge, Bool_t skip, Int_t iDir,
- Bool_t exclude, std::map<Double_t,matrix4> *cache)
- // Bool_t exclude, std::map<Double_t,std::tuple<TMatrixD,TMatrixD,TMatrixD,TMatrixD> > *cache)
- {
- /// Refit track in TPC using track hits (toward beam line) for some
- /// particle mass and charge hypothesis
- //#pragma omp critical
- {
- track->GetHits()->UnSort();
- track->GetHits()->Sort();
- }
- track->SetChi2(0.);
- if (iDir == 1) track->SetDirection(MpdKalmanTrack::kInward);
- //if (GetNofHits() == 1) return;
- track->SetWeight(*track->GetWeightAtHit());
- TMatrixDSym *w = track->GetWeight();
- Int_t nHits = track->GetNofHits(), nHits2 = nHits * nHits;;
- for (Int_t i = 0; i < 5; ++i) {
- for (Int_t j = i; j < 5; ++j) {
- //if (j == i) (*w)(i,j) /= 1000.;
- if (j == i) (*w)(i,j) /= nHits2;
- //if (j == i) (*w)(i,j) /= 1;
- else (*w)(i,j) = (*w)(j,i) = 0.;
- }
- }
- TMatrixD param(5,1);
- TMatrixDSym weight(5), pointWeight(5);
- //Int_t iDir = 1;
- Int_t ibeg = (iDir > 0) ? 0 : nHits-1;
- Int_t iend = (iDir > 0) ? nHits-1 : 0;
- iend += iDir;
- MpdKalmanHit *hit = 0x0;
- track->SetPosNew(track->GetPosAtHit()+0.2);
- track->SetParamNew(*track->GetParamAtHit());
- track->SetNodeNew("");
- Double_t leng = 0.0;
- track->SetLength(leng);
- TString mass2 = "";
- mass2 += mass * mass;
- for (Int_t i = ibeg; i != iend; i+=iDir) {
- hit = (MpdKalmanHit*) track->GetHits()->UncheckedAt(i);
- if (hit->GetUniqueID()) continue; // 5-jan-2018 - skip ITS hits
- //if (!MpdKalmanFilter::Instance()->PropagateToHit(track, hit, kTRUE, kTRUE)) return kFALSE;
- if (!MpdKalmanFilter::Instance()->PropagateToHit(track, hit, kTRUE, kTRUE)) {
- if (fVerbose > 1) cout << " III " << i << " " << hit->GetDist() << " " << track->Pt() << " " << track->GetDirection() << endl;
- return kFALSE;
- }
- Double_t step = track->GetLength() - leng;
- Double_t th = track->GetParamNew(3);
- Double_t cosTh = TMath::Cos(th);
- if (step > 1.e-4) {
- Double_t x0 = 13363.6; // rad. length - TPCMixture
- TMatrixDSym *cov = track->Weight2Cov();
- Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, x0, step, mass2, charge);
- //cout << " Scat: " << hit->GetLayer() << " " << step << " " << TMath::Sqrt(angle2) << endl;
- (*cov)(2,2) += (angle2 / cosTh / cosTh);
- (*cov)(3,3) += angle2;
- Int_t iok = 0;
- MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
- track->SetWeight(*cov);
- }
- Double_t dChi2 = MpdKalmanFilter::Instance()->FilterHit(track, hit, pointWeight, param);
- leng = track->GetLength();
- if (exclude && dChi2 > fgkChi2Cut) {
- track->GetHits()->Remove(hit);
- continue; // too high Chi2
- }
- track->SetChi2(track->GetChi2()+dChi2);
- weight = *track->GetWeight();
- weight += pointWeight;
- if (cache) {
- // For smoother
- TMatrixD fc(MpdKalmanFilter::Instance()->GetJacob(),TMatrixD::kMult,*track->GetWeight()); // DW
- (*cache)[track->GetPosNew()] = matrix4(*track->GetParamNew(),param,fc,weight);
- }
- track->SetWeight(weight);
- track->SetParamNew(param);
- //cout << i << " " << dChi2 << " " << 1./track->GetParamNew(4) << endl;
- }
- if (exclude) track->GetHits()->Compress();
- if (iDir < 0 || skip) return kTRUE; // going outward or do not go to the beam line
- // Go to the beam line
- const Int_t nR = 7;
- const Double_t rad[nR] = {27.00, 27.01, 27.31, 27.32, 33.82, 33.83, 34.13};
- const Double_t dx[nR] = {7.8e-4, 1.194e-2, 7.8e-4, 3.3e-4, 7.8e-4, 1.194e-2, 9.6e-4};
- Bool_t ok = 0;
- MpdKalmanHit hitK;
- hitK.SetType(MpdKalmanHit::kFixedR);
- for (Int_t j = nR-1; j >= 0; --j) {
- hitK.SetPos(rad[j]);
- ok = MpdKalmanFilter::Instance()->PropagateToHit(track,&hitK,kTRUE);
- // Add multiple scattering
- TMatrixDSym *cov = track->Weight2Cov();
- Double_t th = track->GetParamNew(3);
- Double_t cosTh = TMath::Cos(th);
- Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(track, dx[j], mass2, charge);
- (*cov)(2,2) += (angle2 / cosTh / cosTh);
- (*cov)(3,3) += angle2;
- //cov->Print();
- Int_t iok = 0;
- MpdKalmanFilter::Instance()->MnvertLocal(cov->GetMatrixArray(), 5, 5, 5, iok);
- track->SetWeight(*cov);
- }
- //if (!ok) continue;
- /* TMatrixDSym *cov = */track->Weight2Cov();
- // Update track
- track->SetParam(*track->GetParamNew());
- track->SetPos(track->GetPosNew());
- // Correct for energy loss
- //*
- Double_t theta = TMath::PiOver2() - track->GetParam(3);
- Double_t pt = 1. / track->GetParam(4);
- Double_t ptCor = CorrectForLoss(TMath::Abs(pt), theta, mass, charge); // ionization loss correction
- track->SetParam(4,TMath::Sign(1./ptCor,pt));
- track->SetParamNew(*track->GetParam());
- //*/
-
- //*
- hitK.SetPos(0.);
- hitK.SetMeas(0,track->GetParam(2));
- //hit.SetRphi(track->GetParam(2)); // track Phi - check if it is necessary !!!!!!!
- Double_t pos = track->GetPosNew();
- MpdKalmanFilter::Instance()->PropagateToHit(track, &hitK, kTRUE);
- track->SetPos(pos); // save position after passing inner shell
- track->SetParam(*track->GetParamNew()); // !!! track params at PCA
- //track->SetLengAtHit(track->GetLength()-track->GetLengAtHit()); // track length from PCA to the nearest hit
- if (0) {
- TVector3 pos3 = MpdKalmanFilter::Instance()->GetGeo()->GlobalPos(hit);
- cout << " Hit: " << pos3.Pt()*pos3.Phi() << " " << pos3.Pt() << " " << pos3.Z() << endl;
- cout << " Track: " << nHits << " " << track->GetParamNew(0) << " " << track->GetPosNew() << " " << track->GetParamNew(1) << " " << track->GetChi2() << endl;
- }
- return kTRUE;
- }
- //__________________________________________________________________________
- Double_t MpdTpcKalmanFilter::CorrectForLossFluct(Double_t pt, Double_t the, Double_t mass, Int_t charge)
- {
- // Ionization loss correction fluctuations - only for low momentum pions
- // Pion
- const Int_t nP = 4, nThe = 5;
- const Double_t moms[nP] = {0.1, 0.2, 0.3, 0.4}; // p
- const Double_t thes[nThe] = {0.015, 0.315, 0.615, 0.915, 1.215}; // dip angle
- static vector<vector<Double_t> > sigmas(nP,vector<Double_t>(nThe));
- static Int_t first = 1;
- if (first) {
- first = 0;
- sigmas[0][0] = 3.2; sigmas[0][1] = 3.3; sigmas[0][2] = 3.6; sigmas[0][3] = 5.3; sigmas[0][4] = 10.7;
- sigmas[1][0] = 2.8; sigmas[1][1] = 2.4; sigmas[1][2] = 2.1; sigmas[1][3] = 2.3; sigmas[1][4] = 2.8;
- sigmas[2][0] = 1.6; sigmas[2][1] = 1.7; sigmas[2][2] = 1.7; sigmas[2][3] = 1.5; sigmas[2][4] = 2.0;
- sigmas[3][0] = 2.5; sigmas[3][1] = 2.4; sigmas[3][2] = 2.2; sigmas[3][3] = 1.8; sigmas[3][4] = 1.9;
- for (Int_t i = 0; i < nP; ++i) {
- for (Int_t j = 0; j < nThe; ++j) sigmas[i][j] *= 1.e-4;
- }
- }
- if (mass < 0.1 || mass > 0.2) return 0;
- Double_t p = pt / TMath::Sin(the) * charge;
- if (p > 0.4) return 0;
- Double_t e = TMath::Sqrt (p*p + mass*mass), sigt = 0.0;
- sigt = Interp2d(moms, thes, sigmas, p, TMath::Abs(TMath::PiOver2()-the));
- //fprintf(lunTpc,"%f %f %f\n",p,TMath::Abs(TMath::PiOver2()-the),sigt);
- sigt /= TMath::Sin(the);
- return e * sigt / p / p / p / TMath::Sin(the); // Delta (1/pt)
- }
- //__________________________________________________________________________
- Double_t MpdTpcKalmanFilter::CorrectForLoss(Double_t pt, Double_t the, Double_t mass, Int_t charge)
- {
- // Ionization loss correction
- // Pion
- const Int_t nPi = 13; // pi-
- //const Double_t momPi[nPi] = {0.075, 0.100, 0.200, 0.300, 0.400, 0.500, 0.700, 0.900,
- // 1.200, 1.500, 2.000, 2.500, 3.500};
- const Double_t momPi[nPi] = {0.055, 0.091, 0.197, 0.297, 0.398, 0.498, 0.698, 0.898,
- 1.198, 1.498, 1.998, 2.498, 3.498};
- //const Double_t dedxPi[nPi] = {8.773, 5.002, 2.597, 2.259, 2.175, 2.143, 2.129, 2.131,
- const Double_t dedxPi[nPi] = {8.476, 5.050, 2.625, 2.304, 2.175, 2.143, 2.129, 2.131,
- 2.152, 2.154, 2.182, 2.193, 2.239}; // mean of fit Gauss*Landau
- //const Double_t sigmaPi[nPi] = {4.2e-4, 3.4e-4, 2.2e-4, 2.0e-4}; // gaussian fit
- // Kaon
- const Int_t nK = 12; // K-
- //const Double_t momK[nK] = {0.200, 0.250, 0.300, 0.400, 0.500, 0.700, 0.900, 1.200,
- // 1.500, 2.000, 2.500, 3.500};
- const Double_t momK[nK] = {0.167, 0.233, 0.289, 0.393, 0.495, 0.697, 0.897, 1.198,
- 1.498, 1.998, 2.498, 3.498};
- const Double_t dedxK[nK] = {11.40, 7.584, 5.858, 4.128, 3.307, 2.635, 2.394, 2.253,
- 2.314, 2.256, 2.253, 2.249}; // mean of fit Gauss*Landau
- // Proton
- const Int_t nP = 11;
- //const Double_t momP[nP] = {0.300, 0.350, 0.400, 0.500, 0.700, 0.900, 1.200, 1.500, 2.000,
- // 2.500, 3.500};
- const Double_t momP[nP] = {0.239, 0.313, 0.375, 0.485, 0.693, 0.895, 1.196, 1.498, 1.997,
- 2.498, 3.498};
- const Double_t dedxP[nP] = {16.87, 12.31, 9.751, 6.948, 4.460, 3.483, 2.808, 2.508, 2.295,
- 2.221, 2.164}; // mean of fit Gauss*Landau
- // Deuteron
- const Int_t nD = 10;
- //const Double_t momD[nD] = {0.450, 0.500, 0.600, 0.700, 0.900, 1.200, 1.500, 2.000,
- // 2.500, 3.500};
- const Double_t momD[nD] = {0.303, 0.406, 0.547, 0.666, 0.882, 1.190, 1.493, 1.996,
- 2.497, 3.497};
- const Double_t dedxD[nD] = {28.85, 22.02, 15.40, 11.73, 7.962, 5.369, 4.159, 3.178,
- 2.705, 2.298}; // mean of fit Gauss*Landau
- // Triton
- const Int_t nT = 10;
- //const Double_t momT[nT] = {0.600, 0.650, 0.700, 0.800, 1.000, 1.200, 1.500, 2.000,
- // 2.500, 3.500};
- const Double_t momT[nT] = {0.400, 0.512, 0.594, 0.730, 0.962, 1.176, 1.485, 1.992,
- 2.494, 3.496};
- const Double_t dedxT[nT] = {35.02, 27.97, 23.83, 18.43, 12.58, 9.341, 6.835, 4.730,
- 3.662, 2.814}; // mean of fit Gauss*Landau
- // He3
- const Int_t nHe3 = 11;
- //const Double_t momHe3[nHe3] = {0.850, 0.900, 1.000, 1.100, 1.200, 1.400, 1.600, 2.000, 2.500,
- // 3.000, 3.500};
- const Double_t momHe3[nHe3] = {0.455, 0.623, 0.819, 0.967, 1.096, 1.330, 1.548, 1.967, 2.478,
- 2.983, 3.485};
- const Double_t dedxHe3[nHe3] = {89.13, 72.42, 55.85, 45.98, 39.47, 30.53, 25.15, 19.07, 14.73,
- 12.76, 11.42}; // mean of fit Gauss*Landau
- // He4
- const Int_t nHe4 = 9;
- //const Double_t momHe4[nHe4] = {1.050, 1.100, 1.200, 1.400, 1.600, 2.000, 2.500,
- // 3.000, 3.500};
- const Double_t momHe4[nHe4] = {0.606, 0.759, 0.962, 1.256, 1.500, 1.941, 2.464,
- 2.974, 3.479};
- const Double_t dedxHe4[nHe4] = {96.18, 82.70, 66.24, 48.44, 38.43, 27.48, 20.09,
- 16.48, 14.14}; // mean of fit Gauss*Landau
- charge = TMath::Abs(charge);
- Double_t p = pt / TMath::Sin(the) * charge, mass2 = mass * mass;
- Double_t t = TMath::Sqrt (p*p + mass2) - mass, dt = 0.0;
- if (mass < 0.1) dt = 2.9; // 2.9 MeV loss for electrons
- else if (mass < 0.2) {
- if (p < 0.45) {
- dt = 0.0152867 - 0.0905572*p - 1.98303*p*p + 35.3591*p*p*p - 245.426*p*p*p*p + 915.386*p*p*p*p*p
- - 1926.86*p*p*p*p*p*p + 2157.13*p*p*p*p*p*p*p - 999.496*p*p*p*p*p*p*p*p;
- dt *= 1000;
- } else dt = MpdKalmanFilter::Instance()->Interp(nPi, momPi, dedxPi, p); // Pion
- }
- else if (mass < 0.6) dt = MpdKalmanFilter::Instance()->Interp(nK, momK, dedxK, p); // Kaon
- else if (mass < 1.0) dt = MpdKalmanFilter::Instance()->Interp(nP, momP, dedxP, p); // Proton
- else if (mass < 2.0) dt = MpdKalmanFilter::Instance()->Interp(nD, momD, dedxD, p); // Deuteron
- else if (mass < 2.9 && charge == 1) dt = MpdKalmanFilter::Instance()->Interp(nT, momT, dedxT, p); // Triton
- else if (mass < 2.9) dt = MpdKalmanFilter::Instance()->Interp(nHe3, momHe3, dedxHe3, p); // He3
- else if (mass < 3.9) dt = MpdKalmanFilter::Instance()->Interp(nHe4, momHe4, dedxHe4, p); // He4
- dt /= TMath::Sin(the);
- t += dt * 1.e-3;
- Double_t e = mass + t;
- p = TMath::Sqrt (e*e - mass2);
- pt = p * TMath::Sin(the);
- if (charge) pt /= charge;
- return pt;
- }
- //__________________________________________________________________________
- void MpdTpcKalmanFilter::MergeTracks(Int_t ipass)
- {
- // Merge tracks
- multimap<Int_t,Int_t> secMaps[24]; // first hit row# - track index
- Int_t nReco = fTracks->GetEntriesFast();
- for (Int_t i = 0; i < nReco; ++i) {
- MpdTpcKalmanTrack *track = (MpdTpcKalmanTrack*) fTracks->UncheckedAt(i);
- MpdKalmanHit *hit = (MpdKalmanHit*) track->GetTrHits()->Last();
- Int_t isec = fSecGeo->Sector(hit->GetDetectorID());
- if (hit->GetMeas(1) < 0) isec += 12;
- Int_t lay = hit->GetLayer();
- secMaps[isec].insert(pair<Int_t,Int_t>(lay,i));
- }
- // Loop over tracks and try to merge them: at first within one sector,
- // then in adjacent sectors
- TVector3 mom;
- multimap<Int_t,Int_t>::iterator it, it1;
- // Loop over sectors
- for (Int_t isec = 0; isec < 24; ++isec) {
- //cout << " Sector No: " << isec << ", tracks: " << secMaps[isec].size() << endl;
- multimap<Int_t,Int_t> ids;
- Int_t isec1 = isec;
- if (ipass == 2) {
- // Merge sectors
- isec1 = isec + 1;
- if (isec1 == 12) isec1 = 0;
- else if (isec1 == 24) isec1 = 12;
- }
- // Loop over tracks in the sector
- for (it = secMaps[isec].begin(); it != secMaps[isec].end(); ++it) {
- MpdTpcKalmanTrack *tr = (MpdTpcKalmanTrack*) fTracks->UncheckedAt(it->second);
- if (tr->GetTrackID() < 0) continue; // already processed
- Int_t hitLays[70] = {0}, nh = tr->GetNofTrHits();
- if (nh >= 50) continue; // exclude long tracks
- if (ipass == 0 && nh < 10) continue; // at first exclude short tracks
- TClonesArray *hits = tr->GetTrHits();
- for (Int_t ih = 0; ih < nh; ++ih) {
- MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(ih);
- hitLays[hit->GetLayer()] = ih + 1;
- }
- // Loop over tracks in the sector
- it1 = it;
- if (isec1 == isec) ++it1;
- else it1 = secMaps[isec1].begin();
- for ( ; it1 != secMaps[isec1].end(); ++it1) {
- //if (it1 == it) continue;
- MpdTpcKalmanTrack *tr1 = (MpdTpcKalmanTrack*) fTracks->UncheckedAt(it1->second);
- if (tr1->GetTrackID() < 0) continue; // already processed
- Int_t hitLays1[70] = {0}, nh1 = tr1->GetNofTrHits(), nSameLay = 0;
- if (nh1 >= 50) continue; // exclude long tracks
- if (ipass == 0 && nh1 < 10) continue; // at first exclude short tracks
- //if (TMath::Abs(tr->Momentum3().Eta()) < 1.5) {
- if (TMath::Abs(tr->Momentum3().Eta()) < 1.3) {
- if (nh < 10 && nh1 < 10) continue; // both tracks too short
- } else if (nh < 5 && nh1 < 5) continue; // both tracks too short
- TClonesArray *hits1 = tr1->GetTrHits();
- // Do not merge tracks with too different angles or going in the same direction or having different charges
- Double_t dphi = tr->GetParam(2) - tr1->GetParam(2);
- if (TMath::Abs(dphi) > TMath::Pi()) dphi -= TMath::Sign (TMath::TwoPi(),1.*dphi);
- if (ipass == 0 && (TMath::Abs(tr->GetParam(3) - tr1->GetParam(3)) > 0.05 ||
- tr->GetDirection() + tr1->GetDirection() != 1 ||
- TMath::Abs(dphi) > 0.5 || tr1->Charge() - tr->Charge() != 0)) continue;
- else if (ipass == 1 && (TMath::Abs(tr->GetParam(3) - tr1->GetParam(3)) > 0.1 ||
- tr->GetDirection() + tr1->GetDirection() != 1 ||
- TMath::Abs(dphi) > 0.75)) continue;
- else if (ipass == 2 && (TMath::Abs(tr->GetParam(3) - tr1->GetParam(3)) > 0.1 ||
- tr->GetDirection() + tr1->GetDirection() != 1 ||
- TMath::Abs(dphi) > 0.75)) continue;
- for (Int_t ih = 0; ih < nh1; ++ih) {
- MpdKalmanHit *hit = (MpdKalmanHit*) hits1->UncheckedAt(ih);
- hitLays1[hit->GetLayer()] = ih + 1;
- if (hitLays[hit->GetLayer()]) nSameLay += 1; // number of hits on the same rows
- }
- // If there are hits on the same rows, compute average distance between them
- if (nSameLay) {
- // Compute average distance
- Double_t dist = 0;
- TVector3 posLoc, pos1, pos2;
- for (Int_t ih = 0; ih < 70; ++ih) {
- if (hitLays[ih] == 0 || hitLays1[ih] == 0) continue;
- MpdKalmanHit *hit = (MpdKalmanHit*) hits->UncheckedAt(hitLays[ih]-1);
- MpdKalmanHit *hit1 = (MpdKalmanHit*) hits1->UncheckedAt(hitLays1[ih]-1);
- if (isec1 == isec) {
- Double_t distx = hit->GetMeas(0) - hit1->GetMeas(0);
- distx *= distx;
- Double_t errx2 = 1; //hit->GetErr(0) * hit->GetErr(0) + hit1->GetErr(0) * hit1->GetErr(0);
- Double_t distz = hit->GetMeas(1) - hit1->GetMeas(1);
- distz *= distz;
- Double_t errz2 = 1; //hit->GetErr(1) * hit->GetErr(1) + hit1->GetErr(1) * hit1->GetErr(1);
- //chi2 += (distx / errx2 + distz / errz2);
- dist += TMath::Sqrt (distx / errx2 + distz / errz2);
- } else {
- posLoc.SetXYZ(hit->GetMeas(0), hit->GetDist(), hit->GetMeas(1));
- fSecGeo->Local2Global(isec, posLoc, pos1);
- posLoc.SetXYZ(hit1->GetMeas(0), hit1->GetDist(), hit1->GetMeas(1));
- fSecGeo->Local2Global(isec1, posLoc, pos2);
- pos2 -= pos1;
- dist += pos2.Mag();
- }
- }
- dist /= nSameLay;
- if (dist > 5) continue; // do not merge tracks
- }
- // Try to merge tracks
- // Take track with more hits and try to add the second one
- MpdKalmanHit *hit, *hits2[70] = {NULL};
- MpdTpcKalmanTrack *trL = tr, *trS = tr1;
- Int_t *hitLaysL = hitLays, *hitLaysS = hitLays1;
- if (nh < nh1) { trL = tr1; trS = tr; hitLaysL = hitLays1; hitLaysS = hitLays; }
- for (Int_t ih = 0; ih < 70; ++ih) {
- if (hitLaysL[ih]) hit = (MpdKalmanHit*) trL->GetTrHits()->UncheckedAt(hitLaysL[ih]-1);
- else if (hitLaysS[ih]) hit = (MpdKalmanHit*) trS->GetTrHits()->UncheckedAt(hitLaysS[ih]-1);
- else hit = NULL;
- hits2[ih] = hit;
- }
- //cout << " ---------------- " << trL->GetNofTrHits() << " " << trS->GetNofTrHits() << endl;
- MpdTpcKalmanTrack track = *trL;
- MpdKalmanHit hitTmp;
- hitTmp.SetType(MpdKalmanHit::kFixedR);
- //if (track.GetDirection() == MpdKalmanTrack::kInward) MpdKalmanFilter::Instance()->Refit(&track, 1);
- if (track.GetDirection() == MpdKalmanTrack::kInward) {
- track.SetDirection(MpdKalmanTrack::kOutward);
- MpdKalmanFilter::Instance()->Refit(&track, -1);
- }
- //else if (track.GetDirection() == MpdKalmanTrack::kOutward) MpdKalmanFilter::Instance()->Refit(&track, -1);
- //cout << " Refit chi2: " << track.GetChi2() << " " << track.GetNofHits() << " " << trL->GetDirection() << endl;
- if (track.GetChi2() > 1000) continue; // do not merge tracks
- Int_t ok = 0;
- if (trL->GetDirection() == MpdKalmanTrack::kInward) {
- //hitTmp.SetDist(fSecGeo->GetMaxY()/TMath::Cos(fSecGeo->Dphi()/2));
- Double_t radMax = TMath::Max (((MpdKalmanHit*)trL->GetHits()->First())->GetDist(),
- ((MpdKalmanHit*)trS->GetHits()->First())->GetDist());
- radMax += fSecGeo->GetMinY() + 5.0;
- hitTmp.SetDist(radMax);
- ok = MpdKalmanFilter::Instance()->PropagateToHit(&track,&hitTmp,kFALSE);
- } else {
- hitTmp.SetDist(fSecGeo->GetMinY());
- ok = MpdKalmanFilter::Instance()->PropagateToHit(&track,&hitTmp,kFALSE);
- }
- if (!ok) continue; // do not merge tracks
- // Save info for further use
- track.SetWeightAtHit(*track.GetWeight());
- track.SetParamAtHit(*track.GetParamNew());
- track.SetPosAtHit(track.GetPosNew());
- TObjArray *hitss = track.GetHits();
- hitss->Clear();
- for (Int_t ih = 0; ih < 70; ++ih) {
- if (hits2[ih]) hitss->Add(hits2[ih]);
- //if (hits2[ih]) { hitss->Add(hits2[ih]); cout << " dist: " << hits2[ih]->GetLayer() << " " << hits2[ih]->GetDist() << " "
- // << hits2[ih]->GetMeas(0) << " " << hits2[ih]->GetMeas(1) << endl; }
- }
- // Refit with multiple scattering
- if (trL->GetDirection() == MpdKalmanTrack::kInward) {
- if (!Refit(&track, 0.13957, 1, kTRUE, 1, kTRUE)) continue; // do not merge tracks
- //cout << " *** RefitIn ***: " << trL->GetNofTrHits() << " " << track.GetNofHits() << " " << trL->GetChi2()
- // << " " << track.GetChi2() << " " << trL->GetTrackID() << " " << trS->GetTrackID() << endl;
- if (1.0*(track.GetNofHits() - trL->GetNofTrHits()) / trS->GetNofTrHits() < 0.3) continue; // do not merge tracks
- } else if (trL->GetDirection() == MpdKalmanTrack::kOutward) {
- // Refit outward
- if (!Refit(&track, 0.13957, 1, kTRUE, -1, kTRUE) || track.GetNofHits() == 0) continue; // do not merge tracks
- //cout << " *** RefitOut ***: " << trL->GetNofTrHits() << " " << track.GetNofHits() << " " << trL->GetChi2()
- //<< " " << track.GetChi2() << " " << trL->GetTrackID() << " " << trS->GetTrackID() << endl;
-
- //hitTmp.SetDist(fSecGeo->GetMaxY()/TMath::Cos(fSecGeo->Dphi()/2));
- hitTmp.SetDist(((MpdKalmanHit*)track.GetHits()->First())->GetDist() + fSecGeo->GetMinY() + 5.0);
- MpdKalmanFilter::Instance()->PropagateToHit(&track,&hitTmp,kFALSE);
- track.SetPosAtHit(track.GetPosNew());
- track.SetWeightAtHit(*track.GetWeight());
- track.SetParamAtHit(*track.GetParamNew());
- track.SetDirection(MpdKalmanTrack::kInward);
- Refit(&track, 0.13957, 1, kTRUE);
- //cout << " *** RefitOut1 ***: " << trL->GetNofTrHits() << " " << track.GetNofHits() << " " << trL->GetChi2()
- // << " " << track.GetChi2() << " " << trL->GetTrackID() << " " << trS->GetTrackID() << endl;
- if (1.0*(track.GetNofHits() - trL->GetNofTrHits()) / trS->GetNofTrHits() < 0.3) continue; // do not merge tracks
- }
- // Merge tracks
- //cout << track.GetNofTrHits() << " " << track.GetNofHits() << " " << trL->GetNofWrong() << " " << endl;
- MpdTpcKalmanTrack *newTrack = new ((*fTracks)[fTracks->GetEntriesFast()]) MpdTpcKalmanTrack(track);
- //newTrack->SetLengAtHit(saveLeng);
- newTrack->SetParamAtHit(*track.GetParamNew());
- newTrack->SetWeightAtHit(*track.GetWeight());
- newTrack->SetLengAtHit(track.GetLength());
- AddHits(fTracks->GetEntriesFast()-1);
- //cout << newTrack->GetNofTrHits() << " " << newTrack->GetNofWrong() << " " << newTrack->GetNofHits() << endl;
- trS->SetTrackID(-trS->GetTrackID()-1);
- trL->SetTrackID(-trL->GetTrackID()-1);
- break;
- } // for (it1 = it; it1 != secMaps[isec1].end();
- } // for (it = secMaps[isec].begin(); it != secMaps[isec].end();
- } // for (Int_t isec = 0; isec < 24;
- for (Int_t i = 0; i < nReco; ++i) {
- MpdTpcKalmanTrack *track = (MpdTpcKalmanTrack*) fTracks->UncheckedAt(i);
- if (track->GetTrackID() < 0) fTracks->Remove(track);
- }
- fTracks->Compress();
- }
- //__________________________________________________________________________
- Double_t MpdTpcKalmanFilter::Interp2d(const Double_t *moms, const Double_t *thes, vector<vector<Double_t> > &sigmas,
- Double_t p, Double_t dip)
- {
- // 2-d linear interpolation
- Int_t nx = sigmas.size(), ny = sigmas[0].size(), i0 = 0, j0 = 0, i1 = 0, j1 = 0;
- for ( ; i1 < nx; ++i1) if (moms[i1] >= p) break;
- if (i1 == nx) --i1;
- if (i1 == 0) i0 = 1;
- else i0 = i1 - 1;
- Double_t u = (p - moms[i0]) / (moms[i1] - moms[i0]);
- for ( ; j1 < ny; ++j1) if (thes[j1] >= dip) break;
- if (j1 == ny) --j1;
- if (j1 == 0) j0 = 1;
- else j0 = j1 - 1;
- Double_t v = (dip - thes[j0]) / (thes[j1] - thes[j0]);
- Double_t uv = u * v;
- Double_t res = sigmas[i0][j0] * (1 - u - v + uv) + sigmas[i1][j0] * (u - uv) +
- sigmas[i0][j1] * (v - uv) + sigmas[i1][j1] * uv;
- return res;
- }
- //__________________________________________________________________________
- void MpdTpcKalmanFilter::Smooth(MpdTpcKalmanTrack *track, std::vector<std::pair<Double_t,TMatrixD> >& vecSmooth)
- {
- // Smooth track
- fCache->clear();
- Refit(track, 0.13957, 1, kTRUE, 1, kFALSE, fCache); // refit toward beam line
- map<Double_t,matrix4>::iterator mit = fCache->begin(), mit1 = mit;
- for ( ; mit != fCache->end(); ++mit) {
- if (mit == fCache->begin()) {
- vecSmooth.push_back(pair<Double_t,TMatrixD>(mit->first,mit->second.xfilt));
- mit1 = mit;
- continue;
- }
- //*
- TMatrixD wfilt(mit->second.wfilt); // filtered weight
- wfilt.Invert(); // covar
- TMatrixD ak(wfilt,TMatrixD::kMult,mit1->second.jacob); // smoother gain matrix
- TMatrixD dx(vecSmooth.back().second); // smoothed previous
- dx -= (mit1->second.xextr); // -extrap. previous
- TMatrixD xsmooth(ak,TMatrixD::kMult,dx);
- xsmooth += (mit->second.xfilt);
- vecSmooth.push_back(pair<Double_t,TMatrixD>(mit->first,xsmooth));
- //*/
- mit1 = mit;
- }
-
- }
- //__________________________________________________________________________
- ClassImp(MpdTpcKalmanFilter)
|