|
//____________________________________________________________________ // // This is the base class for doing all tracking in the Tpcs // The methods are used by the 2 track classes // BrTpcTrackFollowModule // BrTpcTrackStringModule // //____________________________________________________________________ // // $Id: BrTpcTrackModule.cxx,v 1.10 2002/06/13 18:49:10 videbaek Exp $ // $Author: videbaek $ // $Date: 2002/06/13 18:49:10 $ // $Copyright: (C) 2001 BRAHMS Collaboration <brahmlib@rhic.bnl.gov> // #ifndef BRAT_BrTpcTrackModule #include "BrTpcTrackModule.h" #endif #ifndef ROOT_TDirectory #include "TDirectory.h" #endif #ifndef BRAT_BrTpcHit #include "BrTpcHit.h" #endif #ifndef BRAT_BrParameterDbManager #include "BrParameterDbManager.h" #endif #ifndef BRAT_BrTableNames #include "BrTableNames.h" #endif #ifndef ROOT_TF1 #include "TF1.h" #endif #ifndef ROOT_TH2 #include "TH2.h" #endif #ifndef BRAT_BrTrackResidual #include "BrTrackResidual.h" #endif #ifndef WIN32 #include <iostream> #include <iomanip> #else #include <iostream.h> #include <iomanip.h> #endif //____________________________________________________________________ ClassImp(BrTpcTrackModule); //____________________________________________________________________ BrTpcTrackModule::BrTpcTrackModule() : BrModule() { // Default constructor. DO NOT USE SetState(kSetup); fTrackGroups = 0; fNTrackGroups = 0; fHitTable = 0; fTrackTable = 0; fParams = 0; fHitsInRow = 0; fStatus = 0; // set default parameters SetGhostCutOff(); SetMaxTrackCandidates(); SetUseCL(); } //____________________________________________________________________ BrTpcTrackModule::BrTpcTrackModule(const Char_t* name, const Char_t* title) : BrModule(name, title) { // Named Constructor SetState(kSetup); fTrackGroups = 0; fNTrackGroups = 0; fHitTable = 0; fTrackTable = 0; fParams = 0; fHitsInRow = 0; fStatus = 0; // set default parameters SetGhostCutOff(); SetMaxTrackCandidates(); SetUseCL(); } //____________________________________________________________________ BrTpcTrackModule::~BrTpcTrackModule() { delete fTrackGroups; delete fTrackCandidates; } //____________________________________________________________________ void BrTpcTrackModule::DefineHistograms() { // // The dir should be choosen elsewhere fx BrTpcTrackFollow/String. // Char_t HistName[32]; Char_t HistTitle[32]; sprintf(HistName,"hNtracks_%s",GetName()); sprintf(HistTitle,"%s: # local tracks per event",GetName()); hNtracks = new TH1F(HistName, HistTitle, 60, -0.5, 59.5); hNtracks->SetXTitle( "N_{tracks}" ); hNtracks->SetYTitle( "dN/dN_{tracks}" ); sprintf(HistName,"hNhits_%s",GetName()); sprintf(HistTitle,"%s: # hits per track",GetName()); hNhits = new TH1F(HistName, HistTitle, 13, -0.5, 12.5); hNhits->SetXTitle( "N_{hits}" ); hNhits->SetYTitle( "dN/dN_{hits}" ); sprintf(HistName,"hChi2_%s",GetName()); sprintf(HistTitle,"%s: #chi^{2}/DOF per track",GetName()); hChi2 = new TH1F(HistName, HistTitle, 160, 0.0, 40.0); hChi2->SetXTitle( "#chi^{2}/DOF" ); hChi2->SetYTitle( "dN/d(#chi^{2}/DOF)" ); sprintf(HistName,"hCL_%s",GetName()); sprintf(HistTitle,"%s: CL per track",GetName()); hCL = new TH1F(HistName, HistTitle, 100, 0.0, 1.0); hCL->SetXTitle( "Confidence Level" ); hCL->SetYTitle( "dN/d(CL)" ); sprintf(HistName,"hStatus_%s",GetName()); sprintf(HistTitle,"%s: status (0=Ok, 1=Fail)",GetName()); hStatus = new TH1F(HistName, HistTitle, 2, 0.0, 2.0); hStatus->SetXTitle( "status" ); hStatus->SetYTitle( "dN/d(status)" ); sprintf(HistName,"hVecErrX_%s",GetName()); sprintf(HistTitle,"%s: error on track x slope",GetName()); hVecErrX = new TH1F(HistName, HistTitle, 100, 0, 0.1); hVecErrX->SetXTitle( "#sigma_{dx/dz}" ); hVecErrX->SetYTitle( "dN/d#sigma_{dx/dz}" ); sprintf(HistName,"hVecErrY_%s",GetName()); sprintf(HistTitle,"%s: error on track y slope",GetName()); hVecErrY = new TH1F(HistName, HistTitle, 100, 0, 0.1); hVecErrY->SetXTitle( "#sigma_{dy/dz}" ); hVecErrY->SetYTitle( "dN/d#sigma_{dy/dz}" ); Float_t xmin = fParams->PadToIntrinsic( fParams->GetPadsprow() - 1 ); Float_t xmax = fParams->PadToIntrinsic( 0 ); Float_t ymin = fParams->TimeToIntrinsic( fParams->GetNoBuckets() - 1 ); Float_t ymax = fParams->TimeToIntrinsic( 0 ); sprintf(HistName, "hX_%s", GetName()); sprintf(HistTitle, "%s: track x distribution", GetName()); hX = new TH1F(HistName, HistTitle, 100, xmin, xmax ); hX->SetXTitle( "x [cm]" ); hX->SetYTitle( "dN/dx" ); sprintf(HistName, "hY_%s", GetName()); sprintf(HistTitle, "%s: track y distribution", GetName()); hY = new TH1F(HistName, HistTitle, 100, ymin, ymax ); hY->SetXTitle( "y [cm]" ); hY->SetYTitle( "dN/dx" ); const Float_t low = -1; const Float_t high = 1; sprintf(HistName,"hdevx_%s",GetName()); sprintf(HistTitle,"%s: x residuals all rows",GetName()); hdevx = new TH1F(HistName, HistTitle, 200,low,high); hdevx->SetXTitle( "#delta_{x} [cm]" ); hdevx->SetYTitle( "dN/d#delta_{x}" ); const Int_t nRows = fParams->GetNumberOfRows(); sprintf(HistName,"hdevx_vsrow%s",GetName()); sprintf(HistTitle,"%s: x residuals vs rows",GetName()); hdevxvsrow = new TH2F(HistName, HistTitle, nRows, 0.5, nRows+1, 50,low,high); hdevxvsrow->SetXTitle( "row Number" ); sprintf(HistName,"hdevy_vsrow%s",GetName()); sprintf(HistTitle,"%s: y residuals vs rows",GetName()); hdevyvsrow = new TH2F(HistName, HistTitle, nRows, 0.5, nRows+1, 50,low,high); hdevyvsrow->SetXTitle( "row Number" ); sprintf(HistName,"hdevy_%s",GetName()); sprintf(HistTitle,"%s: y residuals all rows",GetName()); hdevy = new TH1F(HistName, HistTitle, 200,low,high); hdevy->SetXTitle( "#delta_{y} [cm]" ); hdevy->SetYTitle( "dN/d#delta_{y}" ); for(Int_t ir = 0; ir < nRows; ir++){ sprintf(HistName,"hdevRowx%d_%s", ir+1, GetName()); sprintf(HistTitle,"%s: x residuals, row %d", GetName(), ir+1); hdevRowx[ir] = new TH1F(HistName, HistTitle, 100, low, high); hdevRowx[ir]->SetXTitle( "#delta_{x} [cm]" ); hdevRowx[ir]->SetYTitle( "dN/d#delta_{x}" ); sprintf(HistName,"hdevRowy%d_%s", ir+1, GetName()); sprintf(HistTitle,"%s: y residuals, row %d", GetName(), ir+1); hdevRowy[ir] = new TH1F(HistName, HistTitle, 100, low, high); hdevRowy[ir]->SetXTitle( "#delta_{y} [cm]" ); hdevRowy[ir]->SetYTitle( "dN/d#delta_{y}" ); } } //____________________________________________________________________ void BrTpcTrackModule::Init() { // Job-level initialisation SetState(kInit); BrParameterDbManager* paramDb = BrParameterDbManager::Instance(); fParams = (BrDetectorParamsTPC*) paramDb->GetDetectorParameters("BrDetectorParamsTPC", GetName()); fTrackCandidates = new TObjArray(fMaxTrackCandidates); fTrackGroups = new TObjArray(500); const Int_t nRows = fParams->GetNumberOfRows(); fHitsInRow = new TObjArray [nRows]; } //____________________________________________________________________ BrTpcTrackCandidate* BrTpcTrackModule::AddNewTrackCandidate() { // Add new track to candidate list if max is not exceeded if(fNTrackCandidates > fMaxTrackCandidates) { Warning("AddNewTrack","Max track candidates exceeded"); fStatus = 1; return 0; } BrTpcTrackCandidate *track = new BrTpcTrackCandidate(); fTrackCandidates->Add(track); fNTrackCandidates++; Assert(fTrackCandidates->GetEntriesFast() == fNTrackCandidates); return track; } //____________________________________________________________________ BrTpcTrackCandidate* BrTpcTrackModule::AddNewTrackCandidate(BrTpcTrackCandidate *original) { // Add new copy of track to candidate list if max is not exceeded if(fNTrackCandidates > fMaxTrackCandidates) { Warning("AddNewTrack","Max track candidates exceeded"); fStatus = 1; return 0; } BrTpcTrackCandidate *track = new BrTpcTrackCandidate(*original); fTrackCandidates->Add(track); fNTrackCandidates++; Assert(fTrackCandidates->GetEntriesFast() == fNTrackCandidates); return track; } //_____________________________________________________________________________ void BrTpcTrackModule::DeleteBadTracks() { // // Find tracks marked in track table as bad, delete them, // and compress the track table const Int_t nTracks = fNTrackCandidates; for(Int_t i = 0;i < nTracks; i++) { BrTpcTrackCandidate *track = (BrTpcTrackCandidate*)fTrackCandidates->At(i); if( track->IsBad()) { fTrackCandidates->RemoveAt(i); fNTrackCandidates--; delete track; } } fTrackCandidates->Compress(); Assert(fTrackCandidates->GetEntries() == fNTrackCandidates); } //_________________________________________________________________________ void BrTpcTrackModule::Clear(){ // // Delete entries of BrClonesArray(s) owned by this module. // fTrackGroups->Delete(); fNTrackGroups = 0; fTrackCandidates->Delete(); fNTrackCandidates = 0; const Int_t nRows = fParams->GetNumberOfRows(); for(Int_t i = 0; i < nRows; i++) fHitsInRow[i].Clear(); } //_____________________________________________________________________________ void BrTpcTrackModule::FillTrackGroups() { // Create the different groups of tracks group. i.e. those localtracks // that has overlapping hits. Tracks in a track group are considered // identical (only 1 is real, the rest being ghosts). // The ghost busting is done at a later stage when the DetectorTracks // are filled. // For this method to work the ghost cut off should be related // to the number of planes etc. // // The algorithm could be refined since it is only the first track // in a track group that is matched to the rest. // // Loop over all tracks for(Int_t i = 0; i < fNTrackCandidates; i++) { BrTpcTrackCandidate *track1 = (BrTpcTrackCandidate*)fTrackCandidates->At(i); if(track1->GetTrackGroup() != 0) continue; // Is already in a track group // Create a new track group BrTrackGroup *trackGroup = new BrTrackGroup(); fTrackGroups->Add(trackGroup); fNTrackGroups++; // Add track group pointer to track trackGroup->AddTrack(track1); trackGroup->SetID((Int_t)trackGroup); track1->SetTrackGroup(trackGroup); // Now we loop over the remaining tracks to see if there are other // tracks belonging in the track group const Int_t nHitsTrack1 = track1->GetNhit(); // loop over remaining tracks for(Int_t j = i+1; j < fNTrackCandidates; j++) { BrTpcTrackCandidate *track2 = (BrTpcTrackCandidate*)fTrackCandidates->At(j); if(track2->GetTrackGroup() != 0 ) continue; const Int_t nHitsTrack2 = track2->GetNhit(); // Compare hits for ghost Int_t nSame = 0; for(Int_t n1 = 0; n1 < nHitsTrack1; n1++) { BrTpcHit *hit1 = (BrTpcHit*) track1->GetTpcHitAt(n1); for(Int_t n2 = 0; n2 < nHitsTrack2; n2++) { BrTpcHit *hit2 = (BrTpcHit*) track2->GetTpcHitAt(n2); // If the pointer addresses are identical the hit are identical if(hit1 == hit2) nSame++; } } // Decide on sameness of tracks if(nSame >= fGhostCutOff) { track2->SetTrackGroup(trackGroup); trackGroup->AddTrack(track2); } if(nSame == nHitsTrack1 && nSame == nHitsTrack2) { Error("FillTrackGroups", "Problem : two tracks are identical"); } } } if(DebugLevel() > 2) cout << "Number of track groups : " << fNTrackGroups << endl; Assert(fNTrackGroups == fTrackGroups->GetEntries()); } //_____________________________________________________________________________ void BrTpcTrackModule::SelectBestTracks(BrEventNode *outNode) { // // Method to add the best tracks to the output node // This method assumes that the tracks have been grouped into track // groups of similar tracks (tracks sharing many (>= fGhostCutOff) // hits with FillTrackGroups // Among these the track with the lowest chi**2 is selected and // added to the out table // Create the output table and add it to the outNode fTrackTable = new BrDataTable(Form("%s %s", BRTABLENAMES kTpcTrack, GetName())); outNode->AddDataTable(fTrackTable); if(DebugLevel()) cout << "SelectBestTracks - No Track Groups = " << fNTrackGroups << endl; for(Int_t i = 0; i < fNTrackGroups; i++) { BrTrackGroup *trackGroup = (BrTrackGroup*) fTrackGroups->At(i); // Look in the list of loctra connected to the vtrack and pick // the "best" Int_t bestTrackNo = -1; Float_t chiMin = 100000000.; const Int_t nTracks = trackGroup->GetEntries(); for(Int_t j = 0; j < nTracks; j++) { BrTpcTrackCandidate *track = (BrTpcTrackCandidate*)trackGroup->GetTrackAt(j); if(track->GetQuality() < chiMin) { bestTrackNo = j; chiMin = track->GetQuality(); } } if(bestTrackNo == -1) { Error("SelectBestTracks", "No best track"); return; } BrTpcTrackCandidate *bestTrack = (BrTpcTrackCandidate*)trackGroup->GetTrackAt(bestTrackNo); bestTrack->FillResiduals(); fTrackTable->Add(bestTrack); // Remove it from the clones array container fTrackCandidates->Remove(bestTrack); if(HistOn()) { hNhits->Fill(bestTrack->GetNhit()); hChi2->Fill(bestTrack->GetQuality()); TF1 func; const Int_t NDF = (bestTrack->GetNhit()-2)*2; func.SetNumberFitPoints(NDF); func.SetChisquare(NDF*bestTrack->GetQuality()); hCL->Fill(func.GetProb()); hVecErrX->Fill(bestTrack->GetSlopeError(0)); hVecErrY->Fill(bestTrack->GetSlopeError(1)); const BrVector3D pos = bestTrack->GetPos(); hX->Fill( pos[0] ); hY->Fill( pos[1] ); } if(fUseCL) { TF1 func; const Int_t NDF = (bestTrack->GetNhit()-2)*2; func.SetNumberFitPoints(NDF); func.SetChisquare(NDF*bestTrack->GetQuality()); bestTrack->SetQuality(func.GetProb()); } } if(DebugLevel() > 0) { cout<<"Number of Tracks for " << GetName() << " is : " << fNTrackGroups << endl; if(DebugLevel() > 2) { for(Int_t i = 0; i < fNTrackGroups; i++) { BrTpcTrackCandidate *track = (BrTpcTrackCandidate*)fTrackTable->At(i); track->Print(); } } } // !! The rest of the method only fills histograms !! if(HistOn()) { hNtracks->Fill(fNTrackGroups); for(Int_t i = 0; i < fNTrackGroups; i++) { BrTpcTrackCandidate *track = (BrTpcTrackCandidate*)fTrackTable->At(i); Int_t nHit = track->GetNhit(); for(Int_t j = 0; j < nHit; j++) { BrTrackResidual *res = track->GetResidualAt(j); hdevx->Fill(res->GetDx()); hdevy->Fill(res->GetDy()); hdevxvsrow->Fill(res->GetRow(),res->GetDx()); hdevyvsrow->Fill(res->GetRow(),res->GetDy()); hdevRowx[res->GetRow()-1]->Fill(res->GetDx()); hdevRowy[res->GetRow()-1]->Fill(res->GetDy()); } } } } //_____________________________________________________________________________ void BrTpcTrackModule::SortHitsIntoRowTable() { // Organize the list of sorted hits into the TobjArray of // hits per row. // // Algorithm: // Loop over all Hits. Lookup active row number from actual number. // Add Hits to the array of fRowHits // Initial implementation uses rownumbers not active row numbers. This // implies a few more (empty) TObArray that are created but bot used. // const Int_t nHits = fHitTable->GetEntries(); for(Int_t i = 0; i < nHits; i++) { BrTpcHit *hit = (BrTpcHit*)fHitTable->At(i); Int_t rowIndex = hit->GetRow()-1; fHitsInRow[rowIndex].Add(hit); } } //____________________________________________________________________ void BrTpcTrackModule::Print(Option_t* option) const { // Print module information // See BrModule::Print for options. // In addition this module defines the Option: // <fill in here> cout << " Status = " << setw(8) << fStatus << endl; cout << " GhostCutOff = " << setw(8) << fGhostCutOff << endl; cout << " MaxTrackCandidates = " << setw(8) << fMaxTrackCandidates << endl; } //____________________________________________________________________ // // $Log: BrTpcTrackModule.cxx,v $ // Revision 1.10 2002/06/13 18:49:10 videbaek // Add 2d summery devx/y histogram // // Revision 1.9 2002/01/03 19:53:46 cholm // Prepared to use BrTableNames class (or perhaps BrDetectorList) for table names // // Revision 1.8 2001/11/13 10:22:34 pchristi // Changed definition of deviation histograms so it fills for all rows // instead of just instrumented. // // Revision 1.7 2001/11/02 13:54:05 pchristi // Added changes to BrTpcTrackFollowModule that reflects the use of the // new class BrTrackResidual. // // Revision 1.6 2001/10/12 10:57:34 pchristi // Added ned histogram with confidence levels of best track. I also // added option to get CL as quality in best track instead of chisquared. // To use you have to call the method SetUseCL(kTRUE). // // Revision 1.5 2001/08/23 22:50:21 pchristi // Made a BrVector3D const to fix Stevs compile problems. // // Revision 1.4 2001/08/17 16:33:12 jens // Added histograms for track x and y distributions // // Revision 1.3 2001/08/16 15:45:17 jens // // BrTpcClusterModule // BrTpcDeconvoluteClusterModule // BrTpcTrackModule // Changed some variable names for histograms and the title scheme for histogram. // Now histograms have title like e.g. "<TPC name>: max ADC in clusters". // // BrTpcClusterModule // Added the possibility to reset sequences in fSeqTable. BrTpcSequence::fClustNum // is altered in BrTpcClusterModule::FindClusters, and thus made it impossible to // reanalyze an event. I.e. previously a new instance of BrEvent had to be // allocated and filled from input. // // Revision 1.2 2001/07/06 10:24:50 pchristi // Changed assert to the ROOT Assert in TError.h // // Revision 1.1.1.1 2001/06/21 14:55:13 hagel // Initial revision of brat2 // // Revision 1.1 2001/06/17 17:41:51 pchristi // The new tracking classes. // // |
||||||
|
This page automatically generated by script docBrat by Christian Holm |
Copyright ; 2002 BRAHMS Collaboration
<brahmlib@rcf.rhic.bnl.gov>
|