|
#include "BrGeantToTpcTrackCandidate.h" ClassImp(BrGeantToTpcTrackCandidate); ///////////////////////////////////////////////////////////////// // // BrGeantToTpcTrackCandidate // // Module constructing a set of BrTpcTrackCandidates from // BrGeantTracks and BrGeantHits in a TPC // // December 2001 // Mads Mikelsen // // // Gave this class a new name; BrGeantToTpcTrackCandidate, this class // was formerly known as BrGeantToDetectorTracks. This was done to // mach the new track and tracking classes. // // Using module: // // t1_g2tpc_trackmaker = new BrGeantToTpcTrackCandidate("T1","T1 g2track"); // t1_g2tpc_trackmaker->SetDebugLevel(0); // t1_g2tpc_trackmaker->SetMaxChisq(10.); // t1_g2tpc_trackmaker->SetMinDedx(0.001); ... other settable parameters // ..... // t1_g2tpc_trackmaker->Event(geant_node,geant_track_node); //////////////////////////////////////////////////////////////////////// #ifndef ROOT_TDirectory #include "TDirectory.h" #endif //______________________________________________________________________ BrGeantToTpcTrackCandidate::BrGeantToTpcTrackCandidate() : BrModule() { fGeantHits = new TObjArray(); fGeantHits->SetOwner(kTRUE); fTpcTrackCandidates = new TObjArray(); fMinDedx = 0.000001; fEdgeCut = 0.5; fMaxSlopeX = 1.0; fMaxSlopeY = 0.5; fMaxChisq = 1.0; fAllowedMiss = 2; fCutTime = 0; fInputModule = 0; fDetectorBit = 0; } //___________________________________________________________________ BrGeantToTpcTrackCandidate::BrGeantToTpcTrackCandidate(const Char_t* Name, const Char_t* Title) : BrModule(Name,Title) { SetState(kSetup); fGeantHits = new TObjArray(); fGeantHits->SetOwner(kTRUE); fTpcTrackCandidates = new TObjArray(); fMinDedx = 0.000001; fEdgeCut = 0.5; fMaxSlopeX = 1.0; fMaxSlopeY = 0.5; fMaxChisq = 1.0; fAllowedMiss = 2; fCutTime = 0; fInputModule = 0; fDetectorBit = 0; } //__________________________________________________________________ BrGeantToTpcTrackCandidate::~BrGeantToTpcTrackCandidate() { delete fGeantHits; delete fTpcTrackCandidates; } //_____________________________________________________________________ void BrGeantToTpcTrackCandidate::Init() { //This funktion must be called by user before event SetState(kInit); if (DebugLevel()>0 ) cout<< " Initializing "<<GetName()<<endl; BrParameterDbManager* gParamDb = BrParameterDbManager::Instance(); if (DebugLevel()>0 && (gParamDb)) cout << " Found BrParameterDbManager "<< endl; fParams_tpc = (BrDetectorParamsTPC*) gParamDb->GetDetectorParameters("BrDetectorParamsTPC", GetName()); if (DebugLevel()>0 ){ fParams_tpc->ListParameters(); } } //_________________________________________________________________ void BrGeantToTpcTrackCandidate::DefineHistograms() { // Define histograms. if (GetState() != kInit) { Stop("DefineHistograms", "Must be called after Init"); return; } TDirectory* saveDir = gDirectory; TDirectory* histDir = gDirectory->mkdir(Form("GeantToTpc_%s",GetName())); if (!histDir) { Warning("DefineHistograms","Could not create histogram subdirectory"); return; } histDir->cd(); Char_t HistName[48]; Char_t HistTitle[48]; 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, 40, 0.0, 10.0); hChi2->SetXTitle( "#chi^{2}/DOF" ); hChi2->SetYTitle( "dN/d(#chi^{2}/DOF)" ); const Int_t nRows = fParams_tpc->GetNumberOfRows(); const Int_t padsPrRow = fParams_tpc->GetPadsprow(); const Int_t timeBuckets = fParams_tpc->GetNoBuckets(); const Float_t xmin = fParams_tpc->PadToIntrinsic(padsPrRow - 1); const Float_t xmax = fParams_tpc->PadToIntrinsic(0); const Float_t ymin = fParams_tpc->TimeToIntrinsic(timeBuckets - 1); const Float_t ymax = fParams_tpc->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" ); sprintf(HistName,"hPadTime_%s",GetName()); sprintf(HistTitle,"Pad vs Time %s",GetName()); hPadTime = new TH2F(HistName, HistTitle, padsPrRow, 0, padsPrRow, 100, ymin, ymax); hPadTime->SetXTitle( "pad no." ); hPadTime->SetYTitle( "time pos [cm]" ); sprintf(HistName,"hPadRow_%s",GetName()); sprintf(HistTitle,"Pad vs Row %s",GetName()); hPadRow = new TH2F(HistName, HistTitle, padsPrRow, 0, padsPrRow, nRows, 1, nRows+1); hPadRow->SetXTitle( "pad no." ); hPadRow->SetYTitle( "row no." ); gDirectory = saveDir; } //____________________________________________________________________ void BrGeantToTpcTrackCandidate::Clear() { fGeantHits->Clear(); fTpcTrackCandidates->Delete(); } //__________________________________________________________________ void BrGeantToTpcTrackCandidate::Print(Option_t* option) const { // Standard information printout. // // Options: See BrModule::Print // TString opt(option); opt.ToLower(); BrModule::Print(option); if (opt.Contains("d")) { cout << endl << " Original author: Kris Hagel" << endl << " Revisted by: $Author: cholm $" << endl << " Revision date: $Date: 2002/03/20 00:41:22 $" << endl << " Revision number: $Revision: 1.3 $ " << endl << endl << "*************************************************" << endl; if (fInputModule) { cout << " Info for Geant input module: " << endl; fInputModule->Print(option); } cout << " Number of required hits in this TPC: " << fRequiredNHits << endl; cout << " Detector bits = " << fDetectorBit << endl; } } //_________________________________________________________________ Bool_t BrGeantToTpcTrackCandidate::CheckHit(BrGeantHit* hit) { // Check that geanmthit is inside dimensions found in // DetectorParameters.txt // Also check on energy? const Float_t xdim = TMath::Abs(fParams_tpc->PadToIntrinsic(0.0)) - fEdgeCut; const Float_t ymax = fParams_tpc->TimeToIntrinsic((Float_t)fCutTime)+5*fEdgeCut; const Float_t ymin = fParams_tpc-> TimeToIntrinsic((Float_t)fParams_tpc->GetNoBuckets())-5*fEdgeCut; const Float_t x = 0.5*(hit->LocalXPosIn() + hit->LocalXPosOut()); const Float_t xin = TMath::Abs(hit->LocalXPosIn()); const Float_t xout = TMath::Abs(hit->LocalXPosOut()); const Float_t yin = TMath::Abs(hit->LocalYPosIn()); const Float_t yout = TMath::Abs(hit->LocalYPosOut()); const Float_t dedx = hit->Dedx(); const Int_t pad = Int_t(fParams_tpc->IntrinsicToPad(x)); const Int_t row = hit->Isub(); if ((xin<xdim) && (xout<xdim) && (yin>ymin) && (yin<ymax) && (yout>ymin) && (yout<ymax) && (dedx>fMinDedx) && fParams_tpc->IsPadActive(row, pad)){ if (DebugLevel()>2) cout << " Good hit: " << hit << endl; // fill histograms if(HistOn()) { const Float_t time = 0.5*(hit->LocalYPosIn() + hit->LocalYPosOut()); hPadTime->Fill(pad, time); hPadRow->Fill(pad, row); } return kTRUE; } else { if (DebugLevel()>2) cout << " Bad hit: " << hit << endl; return kFALSE; } } //___________________________________________________________________ BrTpcHit* BrGeantToTpcTrackCandidate::GeantToTpcHit(BrGeantHit* hit) { //Defines a BrTpcHit containing information held by a BrGeantHit Float_t dhit_dpos[2]; BrTpcHit* dhit = new BrTpcHit(); dhit->SetX(0.5*(hit->LocalPosIn()[0]+hit->LocalPosOut()[0])); dhit->SetY(0.5*(hit->LocalPosIn()[1]+hit->LocalPosOut()[1])); dhit->SetZ(0.5*(hit->LocalPosIn()[2]+hit->LocalPosOut()[2])); dhit->SetXError(0.02); dhit->SetYError(0.02); return dhit; } //______________________________________________________________________ void BrGeantToTpcTrackCandidate::Event(BrEventNode* InputNode, BrEventNode* OutputNode) { // Runs through all geant traks and findis the hits associated with // each geant track. Then makes a Fit using these hits. // Defines a new TpcTrackCandidate, with the geant track's UniqueID // and the slopes obtained from tpctrackcandidate->Fit() if (GetState()==kSetup) { Warning("Event", "Init() must be called before event"); return; } SetState(kEvent); Int_t NActiveRows = fParams_tpc->GetNumberOfActiveRows(); fRequiredNHits = NActiveRows - fAllowedMiss; if (DebugLevel()>2 ) { cout << "Active rows: " << NActiveRows << endl; cout << " Required no of hits: " << fRequiredNHits << endl; } Char_t TableName[32]; sprintf(TableName,"TpcTrack %s",GetName()); if (DebugLevel()>0 ){ cout << " In Event of: " << GetName() << "," << GetTitle() << endl; } BrDataTable* gtracks = InputNode->GetDataTable("GeantTracks"); if (!gtracks) return; BrGeantTrack* gtrack = 0; BrTpcTrackCandidate* tpc_track = 0; BrGeantHit* gHit = 0; BrTpcHit* dHit = 0; BrVector3D slope; //slope of TpcTrack Char_t hitTable[100]; sprintf(hitTable, "GeantHits %s", GetName()); BrDataTable *ghits = InputNode->GetDataTable(hitTable); if(!ghits) { if(DebugLevel()>0) cout << " BrGeantToTpcTrackCandidates : No hits in " << GetName() << " for this event" << endl; return; } Int_t nTracks = gtracks->GetEntries(); if(DebugLevel()>1) cout << "Max no. of tracks : " << nTracks << endl; //loop over all geant tracks for (Int_t itr=0; itr<nTracks; itr++) { Bool_t hasHits = false; gtrack = (BrGeantTrack*)gtracks->At(itr); Int_t trackID = gtrack->TrackNo(); // Added some code to look though the Geant Hits for the detector Int_t nHits = ghits->GetEntries(); for (Int_t ihit=0; ihit<nHits; ihit++) { gHit = (BrGeantHit*)ghits->At(ihit); if (gHit->TrackNo() == trackID) { if (DebugLevel() > 4) cout << "Hit : " << endl << gHit << endl << "Track : " << endl << gtrack << endl; if (CheckHit(gHit)) { hasHits = true; dHit = GeantToTpcHit(gHit); fGeantHits->Add(dHit); } // if good Geant hit } // if Geant hit on current track } // loop over Geant hits if (DebugLevel()>2 ) { cout << " Number of hits in " << GetName() << ": " << fGeantHits->GetEntries() << " Required number is: " << fRequiredNHits << endl; } // see if track has enough hits if (!hasHits || fGeantHits->GetEntries()<fRequiredNHits) { // delete geant hits in array fGeantHits->Delete(); continue; } // make a tpc track candidate tpc_track = new BrTpcTrackCandidate(); for (Int_t ihit=0; ihit<fGeantHits->GetEntries(); ihit++){ dHit = (BrTpcHit*)fGeantHits->At(ihit); tpc_track->AddHit(dHit); } tpc_track->Fit(); slope=tpc_track->GetSlope(); Float_t slopeX = (Float_t)slope.GetX(); Float_t slopeY = (Float_t)slope.GetY(); Float_t chisq = tpc_track->GetQuality(); tpc_track->SetUniqueID(gtrack->TrackNo()); if (TMath::Abs(slopeX)<fMaxSlopeX && TMath::Abs(slopeY)<fMaxSlopeY && (chisq<fMaxChisq)) fTpcTrackCandidates->Add(tpc_track); else delete tpc_track; // delete geant hits in array fGeantHits->Delete(); } // loop over Geant tracks if (DebugLevel()>0 ){ cout << " Number of TpcTracksCandidates: " << fTpcTrackCandidates->GetEntries() << endl; } // fill number of tracks if(HistOn()) { hNtracks->Fill(fTpcTrackCandidates->GetEntries()); } // add datatable to outnode with tpctracks if there are any if (fTpcTrackCandidates->GetEntries()){ BrDataTable* track_table = new BrDataTable(TableName); OutputNode->AddObject(track_table); TIter nexttrack(fTpcTrackCandidates); BrTpcTrackCandidate* track; while ( (track = (BrTpcTrackCandidate*)nexttrack()) ){ fTpcTrackCandidates->Remove(track); track_table->Add(track); // fill histograms if(HistOn()) { hNhits->Fill(track->GetNhit()); hChi2->Fill(track->GetQuality()); hX->Fill(track->GetPos()[0]); hY->Fill(track->GetPos()[1]); } } } Clear(); } |
||||||
|
This page automatically generated by script docBrat by Christian Holm |
Copyright ; 2002 BRAHMS Collaboration
<brahmlib@rcf.rhic.bnl.gov>
|