/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

DsRpcSensDet.cc
Go to the documentation of this file.
00001 
00007 #include "DsRpcSensDet.h"
00008 #include "DsRpcModel.h"
00009 
00010 #include "Event/SimRpcHit.h"
00011 
00012 #include "Conventions/Detectors.h"
00013 
00014 #include "GiGaCnv/GiGaVolumeUtils.h"
00015 
00016 #include "DetDesc/DetectorElement.h"
00017 
00018 #include "GaudiKernel/IRndmGenSvc.h"
00019 
00020 #include "G4DataHelpers/ITouchableToDetectorElement.h"
00021 
00022 #include "G4Step.hh"
00023 #include "G4HCofThisEvent.hh"
00024 #include "G4SDManager.hh"
00025 #include "G4StepPoint.hh"
00026 #include "G4Region.hh"
00027 #include "G4LogicalVolume.hh"
00028 #include "G4VFastSimulationModel.hh"
00029 #include "G4ProductionCuts.hh"
00030 #include "G4ThreeVector.hh"
00031 #include "G4TouchableHistory.hh"
00032 
00033 //#include "G4OpticalPhoton.hh"
00034 
00035 #include "CLHEP/Units/PhysicalConstants.h"
00036 
00037 #include <string>
00038 #include <vector>
00039 using namespace std;
00040 
00041 
00042 DetectorId::DetectorId_t detector_ids1[] = 
00043 {  DetectorId::kRPC,
00044    (DetectorId::DetectorId_t)-1 };
00045 
00046 Site::Site_t site_ids1[] = 
00047 { Site::kDayaBay, Site::kLingAo, Site::kFar, (Site::Site_t)-1 };
00048 
00049 
00050 DsRpcSensDet::DsRpcSensDet(const std::string& type,
00051                            const std::string& name, 
00052                            const IInterface*  parent)
00053     : G4VSensitiveDetector(name)
00054     , GiGaSensDetBase(type,name,parent)
00055     , m_t2de(0)
00056 {
00057     info() << "DsRpcSensDet (" << type << "/" << name << ") created" << endreq;
00058 
00059     declareProperty("StripLogicalVolume",
00060                     m_stripLogVol="/dd/Geometry/RPC/lvRPCStrip",
00061                     "RPC read-out strip logical volume to which this SD is attached.");
00062 
00063     declareProperty("TouchableToDetelem", m_t2deName = "TH2DE",
00064                     "The ITouchableToDetectorElement to use to resolve sensor.");
00065 
00066     declareProperty("SensorStructures",m_sensorStructures,
00067                     "TDS Paths in which to look for sensor detector elements"
00068                     " using this sensitive detector");
00069 
00070     declareProperty("PackedIdPropertyName",m_idParameter="RpcID",
00071                     "The name of the user property holding the RPC ID.");
00072 
00073 
00074 }
00075 
00076 DsRpcSensDet::~DsRpcSensDet()
00077 {
00078 }
00079 
00080 StatusCode DsRpcSensDet::initialize()
00081 {
00082     this->GiGaSensDetBase::initialize();
00083     info() << "DsRpcSensDet initialize" << endreq;
00084 
00085     m_t2de = tool<ITouchableToDetectorElement>(m_t2deName);
00086 
00087     // define collections.
00088 
00089     // first a catch-all
00090     collectionName.insert("unknown");
00091     for (int isite=0; site_ids1[isite] >= 0; ++isite) {
00092         Site::Site_t site = site_ids1[isite];
00093 
00094         for (int idet=0; detector_ids1[idet] >= 0; ++idet) {
00095             DetectorId::DetectorId_t detid = detector_ids1[idet];
00096 
00097             DayaBay::Detector det(site,detid);
00098 
00099             if (det.bogus()) continue;
00100 
00101             string name=det.detName();
00102             collectionName.insert(name.c_str());
00103             debug() << "Inserted collectionName " << collectionName.size()-1 
00104                     << " = \"" << name << "\"" 
00105                     << "(" << isite << "," << idet << "): "
00106                     << "(" << site << "," << detid << ")"
00107                     << endreq;
00108         }
00109     }
00110 
00111     // Set up a fast sim model for RPC.  Note, this is
00112     // kind of circular, but I find no other way to attach a fast sim
00113     // to a sens det.
00114 
00115     const G4LogicalVolume* lv = GiGaVolumeUtils::findLVolume(m_stripLogVol);
00116     if (!lv) {
00117         error() << "+++++++++++++++++++++++++++++++Failed to get " << m_stripLogVol << endreq;
00118         return StatusCode::FAILURE;
00119     }
00120 
00121     G4Region* rpc_reg = new G4Region("RPC_region");
00122     rpc_reg->AddRootLogicalVolume(const_cast<G4LogicalVolume*>(lv));
00123 
00124     vector<double> cuts; 
00125     cuts.push_back(1.0*mm);
00126     cuts.push_back(1.0*mm);
00127     cuts.push_back(1.0*mm);
00128     rpc_reg->SetProductionCuts(new G4ProductionCuts());
00129     rpc_reg->GetProductionCuts()->SetProductionCuts(cuts);
00130 
00131     DsRpcModel* rpcModel = new DsRpcModel("RpcModel",rpc_reg);
00132     rpcModel = 0;
00133 
00134     // set up random numbers
00135     IRndmGenSvc *rgs = 0;
00136     if (service("RndmGenSvc",rgs,true).isFailure()) {
00137         fatal() << "Failed to get random service" << endreq;
00138         return StatusCode::FAILURE;        
00139     }
00140     StatusCode sc = m_uni.initialize(rgs, Rndm::Flat(0,1));
00141     if (sc.isFailure()) {
00142         fatal() << "Failed to initialize uniform random numbers" << endreq;
00143         return StatusCode::FAILURE;
00144     }
00145 
00146     return StatusCode::SUCCESS;
00147 }
00148 
00149 StatusCode DsRpcSensDet::finalize()
00150 {
00151     info() << "DsRpcSensDet finalize" << endreq;
00152     return this->GiGaSensDetBase::finalize();
00153 }
00154 
00155 
00156 void DsRpcSensDet::Initialize(G4HCofThisEvent* hce)
00157 {
00158     m_hc.clear();
00159 
00160     G4DhHitCollection* hc = new G4DhHitCollection(SensitiveDetectorName,collectionName[0]);
00161     m_hc[0] = hc;
00162     int hcid = G4SDManager::GetSDMpointer()->GetCollectionID(hc);
00163     hce->AddHitsCollection(hcid,hc);
00164 
00165     for (int isite=0; site_ids1[isite] >= 0; ++isite) {
00166         for (int idet=0; detector_ids1[idet] >= 0; ++idet) {
00167             DayaBay::Detector det(site_ids1[isite],detector_ids1[idet]);
00168 
00169             if (det.bogus()) continue;
00170 
00171             string name=det.detName();
00172             debug()<<"det.detName(): "<<name<<"*********************"<<endreq;  // was info()
00173             G4DhHitCollection* hc = new G4DhHitCollection(SensitiveDetectorName,name.c_str());
00174             short int id = det.siteDetPackedData();
00175             debug()<<"det.siteDetPackedData(): "<<id<<"**************"<<endreq; // was info()
00176             m_hc[id] = hc;
00177 
00178             int hcid = G4SDManager::GetSDMpointer()->GetCollectionID(hc);
00179             hce->AddHitsCollection(hcid,hc);
00180             //info()
00181             debug() << "Add hit collection with hcid=" << hcid << ", cached ID=" 
00182                     << (void*)id 
00183                     << " name= \"" << SensitiveDetectorName << "/" << name << "\"" 
00184                     << endreq;
00185         }
00186     }
00187 
00188     debug() << "DsRpcSensDet Initialize, made "
00189            << hce->GetNumberOfCollections() << " collections"
00190            << endreq;
00191     
00192 }
00193 
00194 // Return the SensitiveDetector ID for the given touchable history.
00195 const DetectorElement* DsRpcSensDet::SensDetElem(const G4TouchableHistory& hist)
00196 {
00197     const IDetectorElement* idetelem = 0;
00198     int steps=0;
00199 
00200     if (!hist.GetHistoryDepth()) {
00201         error() << "DsRpcSensDet::SensDetElem given empty touchable history" << endreq;
00202         return 0;
00203     }
00204 
00205     StatusCode sc = 
00206         m_t2de->GetBestDetectorElement(&hist,m_sensorStructures,idetelem,steps);
00207     if (sc.isFailure()) {      // verbose warning
00208         warning() << "Failed to find detector element in:\n";
00209         for (size_t ind=0; ind<m_sensorStructures.size(); ++ind) {
00210             warning() << "\t\t" << m_sensorStructures[ind] << "\n";
00211         }
00212         warning() << "\tfor touchable history:\n";
00213         for (int ind=0; ind < hist.GetHistoryDepth(); ++ind) {
00214             warning() << "\t (" << ind << ") " 
00215                       << hist.GetVolume(ind)->GetName() << "\n";
00216         }
00217         warning() << endreq;
00218         return 0;
00219     }
00220     
00221     return dynamic_cast<const DetectorElement*>(idetelem);
00222 }
00223 
00224 int  DsRpcSensDet::SensDetId(const DetectorElement& de)
00225 {
00226     const DetectorElement* detelem = &de;
00227     //info()<<"+++++++++++++++++++DetectorElement+++++++++"<<detelem<<endreq; 
00228     while (detelem) {
00229        // info()<<"detelem = dynamic_cast<const DetectorElement*>(detelem->parentIDetectorElement())::::::::"<<detelem<<endreq;
00230         if (detelem->params()->exists(m_idParameter)) {
00231             break;
00232         }
00233         detelem = dynamic_cast<const DetectorElement*>(detelem->parentIDetectorElement());
00234     }
00235     if (!detelem) {
00236         warning() << "Could not get RPC detector element starting from " << de << endreq;
00237         return 0;
00238     }
00239 
00240     return detelem->params()->param<int>(m_idParameter);
00241 }
00242 
00243 
00244 bool DsRpcSensDet::ProcessHits(G4Step* step,
00245                                G4TouchableHistory* /*history*/)
00246 {
00247     //if (!step) return false; just crash for now if not defined
00248 
00249     // Find out what detector we are in (RPC)
00250     G4StepPoint* preStepPoint = step->GetPreStepPoint();
00251 
00252     double energyDep = step->GetTotalEnergyDeposit();
00253  
00254     // Do not process non-positive energy deposits 
00255     if (energyDep <= 0.0) {
00256         debug()<< "Hit energy too low: " << energyDep/CLHEP::eV << endreq;
00257         return true;
00258     }
00259 
00260 
00261 
00262 
00263     const G4TouchableHistory* hist = 
00264         dynamic_cast<const G4TouchableHistory*>(preStepPoint->GetTouchable());
00265     if (!hist or !hist->GetHistoryDepth()) {
00266       error() << "ProcessHits: step has no or empty touchable history." 
00267               << " Particle name for track assoc. with this step is " << step->GetTrack()->GetDefinition()->GetParticleName() << endreq;
00268         return false;
00269     }
00270 
00271     const DetectorElement* de = this->SensDetElem(*hist);
00272     if (!de) return false;
00273 
00274     int rpcid = this->SensDetId(*de);
00275     DayaBay::Detector detector(rpcid);
00276  
00277     G4Track* track = step->GetTrack();
00278     double weight = track->GetWeight();
00279     weight = 1.;
00280 
00281 
00282 #if 1
00283     if (!rpcid) {
00284         warning () << "Dumping TouchableHistory: (Particle name for track assoc. with this step is " 
00285                    << step->GetTrack()->GetDefinition()->GetParticleName() <<" .)\n";
00286         for (int ind=0; ind<hist->GetHistoryDepth(); ++ind) {
00287             warning () << "\t (" << ind << ") " << hist->GetVolume(ind)->GetName() 
00288                     << " #" << hist->GetReplicaNumber(ind)
00289                     << " physvol @ " << (void*)hist->GetVolume(ind)
00290                     << "\n";
00291         }
00292         warning() << endreq;
00293     }
00294 
00295     verbose() << "Dumping TouchableHistory:\n";
00296     for (int ind=0; ind<hist->GetHistoryDepth(); ++ind) {
00297         verbose() << "\t (" << ind << ") " << hist->GetVolume(ind)->GetName() 
00298                   << " #" << hist->GetReplicaNumber(ind)
00299                   << " physvol @ " << (void*)hist->GetVolume(ind)
00300                   << "\n";
00301     }
00302     verbose() << endreq;
00303 #endif
00304 
00305 
00306     DayaBay::SimRpcHit* srhit = new DayaBay::SimRpcHit();
00307     srhit->setEnergyDep(energyDep);
00308 
00309     // base hit
00310 
00311     // Time since event created
00312     srhit->setHitTime(preStepPoint->GetGlobalTime());
00313 
00314     //#include "G4NavigationHistory.hh"
00315 
00316     const G4AffineTransform& trans = hist->GetHistory()->GetTopTransform();
00317     const G4ThreeVector& global_pos = preStepPoint->GetPosition();
00318     G4ThreeVector pos = trans.TransformPoint(global_pos);
00319     srhit->setLocalPos(pos);
00320     srhit->setSensDetId(rpcid);
00321     
00322     // rpc hit
00323     // sprit->setDir(...);       // for now
00324     //G4ThreeVector pol = trans.TransformAxis(track->GetPolarization());
00325    // pol = pol.unit();
00326    // G4ThreeVector dir = trans.TransformAxis(track->GetMomentum());
00327     //dir = dir.unit();
00328     //srhit->setPol(pol);
00329     //srhit->setDir(dir);
00330     //srhit->setWavelength(wavelength);
00331     //srhit->setType(0);
00332     //srhit->setWeight(weight);
00333 
00334     // Dump info on the hit and touchable history
00335 #if 0
00336     debug() << "hit: id=" << (void*)pmtid << ", "
00337             << wavelength/CLHEP::nm << " nm, "
00338             << "pol=[" << pol[0] << "," << pol[1] << "," << pol[2] << "] "
00339             << "pos=[" << pos[0] << "," << pos[1] << "," << pos[2] << "] "
00340             << endreq;
00341 
00342     for (int ind=0; ind<hist->GetHistoryDepth(); ++ind) {
00343         verbose () << "\t (" << ind << ") " << hist->GetVolume(ind)->GetName() 
00344                 << " #" << hist->GetReplicaNumber(ind)
00345                 << " physvol @ " << (void*)hist->GetVolume(ind)
00346                 << "\n";
00347     }
00348     verbose() << endreq;
00349 #endif
00350 
00351 
00352     int trackid = track->GetTrackID();
00353     this->StoreHit(srhit,trackid);
00354 
00355     return true;
00356 }
00357 
00358 void DsRpcSensDet::StoreHit(DayaBay::SimRpcHit* hit, int trackid)
00359 {
00360     int did = hit->sensDetId();
00361     DayaBay::Detector det(did);
00362     short int sdid = det.siteDetPackedData();
00363 
00364     G4DhHitCollection* hc = m_hc[sdid];
00365     if (!hc) {
00366         warning() << "Got hit with no hit collection.  ID = " << (void*)did
00367                   << " which is detector: \"" << DayaBay::Detector(did).detName()
00368                   << "\". Storing to the " << collectionName[0] << " collection"
00369                   << endreq;
00370         sdid = 0;
00371         hc = m_hc[sdid];
00372     }
00373 
00374 #if 1
00375     verbose() << "Storing hit RPC: " << (void*)did 
00376               << " from " << DayaBay::Detector(did).detName()
00377               << " in hc #"<<  sdid << " = "
00378               << hit->hitTime()/CLHEP::ns << "[ns] "
00379               << hit->localPos()/CLHEP::cm << "[cm] "
00380               <<"Storing to the "<<collectionName[0]<<"collection****************"
00381               << endreq;
00382 #endif
00383 
00384     hc->insert(new G4DhHit(hit,trackid));
00385 }
00386 
00387 void DsRpcSensDet::EndOfEvent(G4HCofThisEvent* hce)
00388 {
00389     //info() 
00390     debug()<< "Cache has " << m_hc.size() << " collections*******************" << endreq;
00391 
00392     int ncols = hce->GetNumberOfCollections();
00393     debug() << "DsRpcSensDet EndOfEvent " << ncols << " collections\n";
00394     for (int ind=0; ind<ncols; ++ind) {
00395         G4VHitsCollection* hc = hce->GetHC(ind);
00396         debug() << ind << ": " 
00397                 << hc->GetSDname() << "//" << hc->GetName() << " has " 
00398                 << hc->GetSize() << " hits\n";
00399     }
00400     debug() << endreq;
00401 }
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:17:58 for DetSim by doxygen 1.7.4