|
// $Id: BrGeantToDetectorTracks.cxx,v 1.2 2002/02/25 11:26:14 trulsml Exp $ // // $Log: BrGeantToDetectorTracks.cxx,v $ // Revision 1.2 2002/02/25 11:26:14 trulsml // Removed 'SetDebugLevel(2)' in both CTOR methods. // // Revision 1.1.1.1 2001/06/21 14:55:07 hagel // Initial revision of brat2 // // Revision 1.5 2001/06/19 19:08:33 trine // Default Dedx limit changed to 0.000001 // Possibility to take into account "preprocessor" cut in low timebins // (for more realistic definition of "reconstructible" tracks.) // // Revision 1.4 2001/03/07 12:19:15 cholm // * Made the method BrModule::Info() const obsolete in favour of // BrModule::Print(Option_t* option="B") const. // // Revision 1.3 2000/09/04 14:07:23 videbaek // Updated for LinkDef modifications // // Revision 1.2 2000/06/05 13:23:18 trine // Added function BrDetectorHit* GeantToDetectorHit(BrGeantHit* ) // converting BrGeantHits to BrDetectorHits. // Track fits now using BrLocalTrack::Fit() instead of previous BrDetectorTrack* // FitTrackToGHits() // Output track structure now BrDetectorTrack->BrLocalTrack->BrDetectorHits // (same as for reconstruction code.) // Added settable parameters fEdgeCut, fMaxChisq, fMaxSlopeX, fMaxSlopeY // // #include "BrGeantToDetectorTracks.h" ClassImp(BrGeantToDetectorTracks) /////////////////////////////////////////////////////////////////////////////////// // // // BrGeantToDetectorTracks // // Module constructing a set of BrDetectorTracks from BrGeantTracks // // in one detector (TPC or DC). // // // // March 2000 // // trine // // // // Removed need for info about DetectorIDs (BrGeantInput), so one // // can use the geantdig and then this class. Also added Set function // // for the number of missed rows // // Peter, May 2000 // // // // Introduced new function converting BrGeantHits to BrDetectorHits: // // BrDetectorHit* BrGeantToDetectorTracks::GeantToDetectorHit(BrGeantHit* hit) // // after checking hit quality: // // Bool_t BrGeantToDetectorTracks::CheckHit(BrGeantHit* hit) // // For the moment: Using dpos[0] = dpos[1] = 0.02 for the BrDetectorHits. // // Using BrLocalTrack::Fit() instead of the now obsolete // // BrDetectorTrack* BrGeantToDetectorTracks::FitTrackToGHits(). // // Cuts on track chisquare and track slopes (settable thresholds.) // // Outputs BrDetectorTracks with pointers to BrLocalTracks and associated // // hits. // // Trine, May 2000 // // // // Using module: // // t1_g2dtrackmaker = new BrGeantToDetectorTracks("TPM1","TPM1 g2d maker"); // // t1_g2dtrackmaker->SetDebugLevel(0); // // t1_g2dtrackmaker->SetMaxChisq(10.); // // t1_g2dtrackmaker->SetMinDedx(0.001); ... other settable parameters // // ..... // // t1_g2dtrackmaker->Event(geant_node,geant_track_node); // // // /////////////////////////////////////////////////////////////////////////////////// BrGeantToDetectorTracks::BrGeantToDetectorTracks() : BrModule() { fGeantHits = new TObjArray(); fDetectorTracks = 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; } BrGeantToDetectorTracks::BrGeantToDetectorTracks(Char_t* Name,Char_t* Title) : BrModule(Name,Title) { fGeantHits = new TObjArray(); fDetectorTracks = 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; if ((!strcmp(GetName(),"T3"))||(!strcmp(GetName(),"T4"))||(!strcmp(GetName(),"T5"))){ fIsTPC = kFALSE; if (DebugLevel()>0 ){ cout << GetName() << " This is a DC! " << endl; } } else { fIsTPC = kTRUE; if (DebugLevel()>0 ){ cout << GetName() << " This is a TPC! " << endl; } } Init(); } BrGeantToDetectorTracks::~BrGeantToDetectorTracks() { delete fGeantHits; delete fDetectorTracks; } void BrGeantToDetectorTracks::Init() { if (DebugLevel()>0 ){ cout << " Initializing " << GetName() << endl; } BrParameterDbManager* gParamDb = BrParameterDbManager::Instance(); if (DebugLevel()>0 && (gParamDb)){ cout << " Found BrParameterDbManager " << endl; } if (fIsTPC) { fParams_tpc = (BrDetectorParamsTPC*) gParamDb->GetDetectorParameters("BrDetectorParamsTPC", GetName()); if (DebugLevel()>0 ){ fParams_tpc->ListParameters(); } } else { fParams_dc = (BrDetectorParamsDC*) gParamDb->GetDetectorParameters("BrDetectorParamsDC", GetName()); } } void BrGeantToDetectorTracks::SetInputModule(BrGeantInput *input) { fInputModule = input; if (fIsTPC) { if (DebugLevel()>0 ){ fParams_tpc->ListParameters(); } fDetectorBit |= fInputModule->GetDetectorBit(GetName()); if (DebugLevel()>0 ){ cout << " fDetectorBit for " << GetName() << ": " << fDetectorBit << endl; } } else { Char_t* SubName = new Char_t[4]; for (Int_t i=0; i<3; i++){ sprintf(SubName,"%sA%d",GetName(),i+1); fDetectorBit |= fInputModule->GetDetectorBit((Text_t*)SubName); if (DebugLevel()>0 ){ cout << " fDetectorBit for " << SubName << ": " << fDetectorBit << endl; } } } } void BrGeantToDetectorTracks::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: trulsml $" << endl << " Revision date: $Date: 2002/02/25 11:26:14 $" << endl << " Revision number: $Revision: 1.2 $ " << 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; } } void BrGeantToDetectorTracks::Clear() { fGeantHits->Clear(); fDetectorTracks->Delete(); } void BrGeantToDetectorTracks::Event(BrEventNode* InputNode, BrEventNode* OutputNode) { if (fIsTPC) { Int_t NActiveRows = fParams_tpc->GetNumberOfActiveRows(); fRequiredNHits = NActiveRows - fAllowedMiss; if (DebugLevel()>2 ){ cout << "Active rows: " << NActiveRows << endl; cout << " Required no of hits: " << fRequiredNHits << endl; } } else { Int_t NActiveRows = fParams_dc->GetNplane(); fRequiredNHits = NActiveRows - fAllowedMiss; } Char_t TableName[32]; sprintf(TableName,"DetectorTrack %s",GetName()); if (DebugLevel()>0 ){ cout << " In Event of: " << GetName() << "," << GetTitle() << endl; } BrDataTable* gtracks = InputNode->GetDataTable("GeantTracks"); if (!gtracks) return; BrGeantTrack* gtrack = 0; BrDetectorTrack* det_track = 0; BrLocalTrack* loc_track = 0; BrGeantHit* gHit = 0; BrDetectorHit* dHit = 0; Char_t hitTable[100]; sprintf(hitTable, "GeantHits %s", GetName()); BrDataTable *ghits = InputNode->GetDataTable(hitTable); if(!hitTable) { if(DebugLevel()>0) cout << " BrGeantToDetectorTracks : No hits in " << GetName() << " for this event" << endl; return; } Int_t nTracks = gtracks->GetEntries(); if(DebugLevel()>1) cout << "Max no. of tracks : " << nTracks << endl; 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 = GeantToDetectorHit(gHit); fGeantHits->Add(dHit); } // if good Geant hit } // if Geant hit on current track } // loop over Geant hits if (hasHits){ if (DebugLevel()>2 ){ cout << " Number of hits in " << GetName() << ": " << fGeantHits->GetEntries() << " Required number is: " << fRequiredNHits << endl; } if (fGeantHits->GetEntries()>=fRequiredNHits){ loc_track = new BrLocalTrack(); for (Int_t ihit=0; ihit<fGeantHits->GetEntries(); ihit++){ dHit = (BrDetectorHit*)fGeantHits->At(ihit); loc_track->AddDetectorHit(dHit); } loc_track->Fit(); Float_t slopeX = (Float_t)loc_track->GetVec()[0]; Float_t slopeY = (Float_t)loc_track->GetVec()[1]; Float_t chisq = loc_track->GetChisq(); det_track = new BrDetectorTrack(); det_track->SetLocalTrack(loc_track); det_track->SetPos(loc_track->GetPos()); det_track->SetAlpha(loc_track->GetVec()); det_track->SetID(gtrack->TrackNo()); if (DebugLevel()> 2){ cout << " Geant Track properties: " << det_track->GetAlpha() << " " << det_track->GetPos() << " " << chisq << endl; } if (TMath::Abs(slopeX)<fMaxSlopeX && TMath::Abs(slopeY)<fMaxSlopeY && (chisq<fMaxChisq)){ fDetectorTracks->Add(det_track); if (DebugLevel()> 2){ cout << " Good Geant Track properties: " << det_track->GetAlpha() << " " << det_track->GetPos() << " " << chisq << endl; } } // if good track properties } // if sufficient number of good hits } // if track has hits in this detector fGeantHits->Clear(); } // loop over Geant tracks if (DebugLevel()>0 ){ cout << " Number of DetectorTracks: " << fDetectorTracks->GetEntries() << endl; } if (fDetectorTracks->GetEntries()){ BrDataTable* track_table = new BrDataTable(TableName); OutputNode->AddObject(track_table); TIter nexttrack(fDetectorTracks); BrDetectorTrack* track; while ( (track = (BrDetectorTrack*)nexttrack()) ){ fDetectorTracks->Remove(track); track_table->Add(track); } } Clear(); } Bool_t BrGeantToDetectorTracks::CheckHit(BrGeantHit* hit) { // Check that hit is inside dimensions found in DetectorParameters.txt // Also check on energy? Float_t xdim = TMath::Abs(fParams_tpc->PadToIntrinsic(0.0)) - fEdgeCut; // Float_t ydim = TMath::Abs(fParams_tpc->TimeToIntrinsic(0.0)) - fEdgeCut; Float_t ymax = fParams_tpc->TimeToIntrinsic((Float_t)fCutTime)+5*fEdgeCut; Float_t ymin = fParams_tpc->TimeToIntrinsic((Float_t)fParams_tpc->GetNoBuckets())-5*fEdgeCut; Float_t xin = TMath::Abs(hit->LocalXPosIn()); Float_t xout = TMath::Abs(hit->LocalXPosOut()); Float_t yin = TMath::Abs(hit->LocalYPosIn()); Float_t yout = TMath::Abs(hit->LocalYPosOut()); Float_t dedx = hit->Dedx(); if ((xin<xdim) && (xout<xdim) && (yin>ymin) && (yin<ymax) && (yout>ymin) && (yout<ymax) && (dedx>fMinDedx)){ if (DebugLevel()>2) cout << " Good hit: " << hit << endl; return kTRUE; } else { if (DebugLevel()>2) cout << " Bad hit: " << hit << endl; return kFALSE; } } BrDetectorHit* BrGeantToDetectorTracks::GeantToDetectorHit(BrGeantHit* hit) { Float_t dhit_pos[3]; Float_t dhit_dpos[2]; BrDetectorHit* dhit = new BrDetectorHit(); dhit_pos[0] = 0.5*(hit->LocalPosIn()[0]+hit->LocalPosOut()[0]); dhit_pos[1] = 0.5*(hit->LocalPosIn()[1]+hit->LocalPosOut()[1]); dhit_pos[2] = 0.5*(hit->LocalPosIn()[2]+hit->LocalPosOut()[2]); dhit->SetPos(dhit_pos); dhit_dpos[0] = 0.02; dhit_dpos[1] = 0.02; dhit->SetDpos(dhit_dpos); Int_t Idedx = (Int_t)(100000*hit->Dedx()); dhit->SetStat(Idedx); if (DebugLevel()>2) cout << " Detector hit: " << dhit << endl; return dhit; } BrDetectorTrack* BrGeantToDetectorTracks::FitTrackToGHits() { Double_t s_z2 = 0.0; Double_t s_z = 0.0; Double_t s_s = 0.0; Double_t s_zx = 0.0; Double_t s_x = 0.0; Double_t s_zy = 0.0; Double_t s_y = 0.0; Double_t ax00 = 0.0; Double_t x00 = 0.0; Double_t ay00 = 0.0; Double_t y00 = 0.0; BrGeantHit* hit = 0; BrVector3D Posin(0.0,0.0,0.0); BrVector3D Posout(0.0,0.0,0.0); for (Int_t ihit=0; ihit<fGeantHits->GetEntries(); ihit++){ hit =(BrGeantHit*)fGeantHits->At(ihit); Posin = hit->LocalPosIn(); Posout = hit->LocalPosOut(); s_z2 += Posin[2]*Posin[2]; s_z2 += Posout[2]*Posout[2]; s_z += Posin[2]; s_z += Posout[2]; s_s += 2; s_zx += Posin[0]*Posin[2]; s_zx += Posout[0]*Posout[2]; s_zy += Posin[1]*Posin[2]; s_zy += Posout[1]*Posout[2]; s_x += Posin[0]; s_x += Posout[0]; s_y += Posin[1]; s_y += Posout[1]; } if ((s_z2*s_s - s_z*s_z) != 0.0){ ax00 = (s_zx*s_s - s_z*s_x)/(s_z2*s_s - s_z*s_z); x00 = (s_z2*s_x - s_z*s_zx)/(s_z2*s_s - s_z*s_z); ay00 = (s_zy*s_s - s_z*s_y)/(s_z2*s_s - s_z*s_z); y00 = (s_z2*s_y - s_z*s_zy)/(s_z2*s_s - s_z*s_z); } BrVector3D testvec(ax00,ay00,(Double_t)1.0); BrVector3D testpt(x00,y00,(Double_t)0.0); BrDetectorTrack* track = new BrDetectorTrack(); track->SetPos(testpt); track->SetAlpha(testvec); Float_t Quality = 0.0; // Calculate Chi**2 or rather just the average offset for (Int_t jhit=0; jhit<fGeantHits->GetEntries(); jhit++){ hit =(BrGeantHit*)fGeantHits->At(jhit); Posin = hit->LocalPosIn(); Posout = hit->LocalPosOut(); Float_t x = ( Posin[ 0 ] + Posout[ 0 ] ) / (Float_t) 2.0; Float_t y = ( Posin[ 1 ] + Posout[ 1 ] ) / (Float_t) 2.0; Float_t z = ( Posin[ 2 ] + Posout[ 2 ] ) / (Float_t) 2.0; Float_t xfit = testvec[ 0 ] * z + testpt[ 0 ]; Float_t yfit = testvec[ 1 ] * z + testpt[ 1 ]; Quality += sqrt( pow( x - xfit, 2.0 ) + pow( y - yfit, 2.0 ) ); } Quality /= (Float_t) fGeantHits->GetEntries(); track->SetQuality( Quality ); if (DebugLevel()>2){ cout << " Dir: " << track->GetAlpha() << " Pos: " << track->GetPos() << " Hits: " << fGeantHits->GetEntries() << endl; } return track; } |
||||||
|
This page automatically generated by script docBrat by Christian Holm |
Copyright ; 2002 BRAHMS Collaboration
<brahmlib@rcf.rhic.bnl.gov>
|