00001 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
00002 #include "SimG4Core/Notification/interface/TrackWithHistory.h"
00003 #include "SimG4Core/Notification/interface/TrackInformation.h"
00004
00005 #include "SimG4CMS/Calo/interface/CaloTrkProcessing.h"
00006
00007 #include "SimG4Core/Application/interface/SimTrackManager.h"
00008 #include "DetectorDescription/Core/interface/DDFilter.h"
00009 #include "DetectorDescription/Core/interface/DDFilteredView.h"
00010 #include "DetectorDescription/Core/interface/DDSolid.h"
00011 #include "DetectorDescription/Core/interface/DDValue.h"
00012 #include "FWCore/Utilities/interface/Exception.h"
00013
00014 #include "G4EventManager.hh"
00015 #include "G4LogicalVolumeStore.hh"
00016 #include "G4LogicalVolume.hh"
00017 #include "G4Step.hh"
00018 #include "G4Track.hh"
00019
00020
00021
00022 CaloTrkProcessing::CaloTrkProcessing(G4String name,
00023 const DDCompactView & cpv,
00024 SensitiveDetectorCatalog & clg,
00025 edm::ParameterSet const & p,
00026 const SimTrackManager* manager) :
00027 SensitiveCaloDetector(name, cpv, clg, p), lastTrackID(-1),
00028 m_trackManager(manager) {
00029
00030
00031 edm::ParameterSet m_p = p.getParameter<edm::ParameterSet>("CaloTrkProcessing");
00032 testBeam = m_p.getParameter<bool>("TestBeam");
00033 eMin = m_p.getParameter<double>("EminTrack")*MeV;
00034
00035 edm::LogInfo("CaloSim") << "CaloTrkProcessing: Initailised with TestBeam = "
00036 << testBeam << " Emin = " << eMin << " MeV";
00037
00038
00039 G4String attribute = "ReadOutName";
00040 DDSpecificsFilter filter;
00041 DDValue ddv(attribute,name,0);
00042 filter.setCriteria(ddv,DDSpecificsFilter::equals);
00043 DDFilteredView fv(cpv);
00044 fv.addFilter(filter);
00045 fv.firstChild();
00046 DDsvalues_type sv(fv.mergedSpecifics());
00047
00048 G4String value = "Calorimeter";
00049 std::vector<std::string> caloNames = getNames (value, sv);
00050 #ifdef ddebug
00051 LogDebug("CaloSim") << "CaloTrkProcessing: Names for " << value << ":";
00052 for (unsigned int i=0; i<caloNames.size(); i++)
00053 LogDebug("CaloSim") << " (" << i << ") " << caloNames[i];
00054 #endif
00055
00056 value = "Levels";
00057 std::vector<double> levels = getNumbers (value, sv);
00058 #ifdef ddebug
00059 LogDebug("CaloSim") << "CaloTrkProcessing: Names for " << value << ":";
00060 for (unsigned int i=0; i<levels.size(); i++)
00061 LogDebug("CaloSim") << " (" << i << ") " << levels[i];
00062 #endif
00063
00064 value = "Neighbours";
00065 std::vector<double> neighbours = getNumbers (value, sv);
00066 #ifdef ddebug
00067 LogDebug("CaloSim") << "CaloTrkProcessing: Names for " << value << ":";
00068 for (unsigned int i=0; i<neighbours.size(); i++)
00069 LogDebug("CaloSim") << " (" << i << ") " << neighbours[i];
00070 #endif
00071
00072 value = "Inside";
00073 std::vector<std::string> insideNames = getNames (value, sv);
00074 #ifdef ddebug
00075 LogDebug("CaloSim") << "CaloTrkProcessing: Names for " << value << ":";
00076 for (unsigned int i=0; i<insideNames.size(); i++)
00077 LogDebug("CaloSim") << " (" << i << ") " << insideNames[i];
00078 #endif
00079
00080 value = "InsideLevel";
00081 std::vector<double> insideLevel = getNumbers (value, sv);
00082 #ifdef ddebug
00083 LogDebug("CaloSim") << "CaloTrkProcessing: Names for " << value << ":";
00084 for (unsigned int i=0; i<insideLevel.size(); i++)
00085 LogDebug("CaloSim") << " (" << i << ") " << insideLevel[i];
00086 #endif
00087
00088 if (caloNames.size() < neighbours.size()) {
00089 edm::LogError("CaloSim") << "CaloTrkProcessing: # of Calorimeter bins "
00090 << caloNames.size() << " does not match with "
00091 << neighbours.size() << " ==> illegal ";
00092 throw cms::Exception("Unknown", "CaloTrkProcessing")
00093 << "Calorimeter array size does not match with size of neighbours\n";
00094 }
00095
00096 const G4LogicalVolumeStore * lvs = G4LogicalVolumeStore::GetInstance();
00097 std::vector<G4LogicalVolume *>::const_iterator lvcite;
00098 int istart = 0;
00099 for (unsigned int i=0; i<caloNames.size(); i++) {
00100 G4LogicalVolume* lv = 0;
00101 G4String name = caloNames[i];
00102 for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++)
00103 if ((*lvcite)->GetName() == name) {
00104 lv = (*lvcite);
00105 break;
00106 }
00107 CaloTrkProcessing::Detector detector;
00108 detector.name = name;
00109 detector.lv = lv;
00110 detector.level = static_cast<int>(levels[i]);
00111 int number = static_cast<int>(neighbours[i]);
00112 if (istart+number > (int)(insideNames.size())) {
00113 edm::LogError("CaloSim") << "CaloTrkProcessing: # of InsideNames bins "
00114 << insideNames.size() << " too few compaerd to "
00115 << istart+number << " requested ==> illegal ";
00116 throw cms::Exception("Unknown", "CaloTrkProcessing")
00117 << "InsideNames array size does not match with list of neighbours\n";
00118 }
00119 std::vector<std::string> inside;
00120 std::vector<G4LogicalVolume*> insideLV;
00121 std::vector<int> insideLevels;
00122 for (int k = 0; k < number; k++) {
00123 lv = 0;
00124 name = insideNames[istart+k];
00125 for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++)
00126 if ((*lvcite)->GetName() == name) {
00127 lv = (*lvcite);
00128 break;
00129 }
00130 inside.push_back(name);
00131 insideLV.push_back(lv);
00132 insideLevels.push_back(static_cast<int>(insideLevel[istart+k]));
00133 }
00134 detector.fromDets = inside;
00135 detector.fromDetL = insideLV;
00136 detector.fromLevels = insideLevels;
00137 detectors.push_back(detector);
00138 istart += number;
00139 }
00140
00141 edm::LogInfo("CaloSim") << "CaloTrkProcessing: with " << detectors.size()
00142 << " calorimetric volumes";
00143 for (unsigned int i=0; i<detectors.size(); i++) {
00144 edm::LogInfo("CaloSim") << "CaloTrkProcessing: Calorimeter volume " << i
00145 << " " << detectors[i].name << " LV "
00146 << detectors[i].lv << " at level "
00147 << detectors[i].level << " with "
00148 << detectors[i].fromDets.size() << " neighbours";
00149 for (unsigned int k=0; k<detectors[i].fromDets.size(); k++)
00150 edm::LogInfo("CaloSim") << " Element " << k << " "
00151 << detectors[i].fromDets[k] << " LV "
00152 << detectors[i].fromDetL[k] << " at level "
00153 << detectors[i].fromLevels[k];
00154 }
00155 }
00156
00157 CaloTrkProcessing::~CaloTrkProcessing() {
00158 edm::LogInfo("CaloSim") << "CaloTrkProcessing: Deleted";
00159 }
00160
00161 void CaloTrkProcessing::update(const BeginOfEvent * evt) {
00162 lastTrackID = -1;
00163 }
00164
00165 void CaloTrkProcessing::update(const G4Step * aStep) {
00166
00167
00168
00169 G4Track* theTrack = aStep->GetTrack();
00170 int id = theTrack->GetTrackID();
00171
00172 TrackInformation* trkInfo = dynamic_cast<TrackInformation*>
00173 (theTrack->GetUserInformation());
00174
00175 if (trkInfo == 0) {
00176 edm::LogError("CaloSim") << "CaloTrkProcessing: No trk info !!!! abort ";
00177 throw cms::Exception("Unknown", "CaloTrkProcessing")
00178 << "cannot get trkInfo for Track " << id << "\n";
00179 }
00180
00181 if (testBeam) {
00182 if (trkInfo->getIDonCaloSurface() == 0) {
00183 #ifdef ddebug
00184 LogDebug("CaloSim") << "CaloTrkProcessing set IDonCaloSurface to " << id
00185 << " at step Number "
00186 << theTrack->GetCurrentStepNumber();
00187 #endif
00188 trkInfo->setIDonCaloSurface(id,0,0);
00189 lastTrackID = id;
00190 if (theTrack->GetKineticEnergy()/MeV > eMin)
00191 trkInfo->putInHistory();
00192 }
00193 } else {
00194 #ifdef ddebug
00195 LogDebug("CaloSim") << "CaloTrkProcessing Entered for " << id
00196 << " at stepNumber "<< theTrack->GetCurrentStepNumber()
00197 << " IDonCaloSur.. " << trkInfo->getIDonCaloSurface()
00198 << " CaloCheck " << trkInfo->caloIDChecked();
00199 #endif
00200 if (trkInfo->getIDonCaloSurface() != 0) {
00201 if (trkInfo->caloIDChecked() == false) {
00202 G4StepPoint* postStepPoint = aStep->GetPostStepPoint();
00203 const G4VTouchable* post_touch = postStepPoint->GetTouchable();
00204
00205 if (isItInside(post_touch, trkInfo->getIDCaloVolume(),
00206 trkInfo->getIDLastVolume()) > 0) {
00207 trkInfo->setIDonCaloSurface(0, -1, -1);
00208 } else {
00209 trkInfo->setCaloIDChecked(true);
00210 }
00211 }
00212 } else {
00213 G4StepPoint* postStepPoint = aStep->GetPostStepPoint();
00214 const G4VTouchable* post_touch = postStepPoint->GetTouchable();
00215 int ical = isItCalo(post_touch);
00216 if (ical >= 0) {
00217 G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
00218 const G4VTouchable* pre_touch = preStepPoint->GetTouchable();
00219 int inside = isItInside(pre_touch, ical, -1);
00220 if (inside >= 0 || (theTrack->GetCurrentStepNumber()==1)) {
00221 trkInfo->setIDonCaloSurface(id,ical,inside);
00222 trkInfo->setCaloIDChecked(true);
00223 lastTrackID = id;
00224 if (theTrack->GetKineticEnergy()/MeV > eMin) trkInfo->putInHistory();
00225 #ifdef ddebug
00226 LogDebug("CaloSim") <<"CaloTrkProcessing: set ID on Calo " << ical
00227 << " surface (Inside " << inside << ") to "
00228 << id << " of a Track with Kinetic Energy "
00229 << theTrack->GetKineticEnergy()/MeV << " MeV";
00230 #endif
00231 }
00232 }
00233 }
00234 }
00235 }
00236
00237 std::vector<std::string> CaloTrkProcessing::getNames(const G4String str,
00238 const DDsvalues_type &sv){
00239
00240 LogDebug("CaloSim") << "CaloTrkProcessing::getNames called for " << str;
00241
00242 DDValue value(str);
00243 if (DDfetch(&sv,value)) {
00244 LogDebug("CaloSim") << value;
00245 const std::vector<std::string> & fvec = value.strings();
00246 int nval = fvec.size();
00247 if (nval < 1) {
00248
00249 edm::LogError("CaloSim") << "CaloTrkProcessing: # of " << str
00250 << " bins " << nval << " < 1 ==> illegal ";
00251 throw cms::Exception("Unknown", "CaloTrkProcessing")
00252 << "nval < 2 for array " << str << "\n";
00253 }
00254
00255 return fvec;
00256 } else {
00257 edm::LogError("CaloSim") << "CaloTrkProcessing: cannot get array " << str ;
00258 throw cms::Exception("Unknown", "CaloTrkProcessing")
00259 << "cannot get array " << str << "\n";
00260 }
00261 }
00262
00263 std::vector<double> CaloTrkProcessing::getNumbers(const G4String str,
00264 const DDsvalues_type &sv) {
00265
00266 LogDebug("CaloSim") << "CaloTrkProcessing::getNumbers called for " << str;
00267
00268 DDValue value(str);
00269 if (DDfetch(&sv,value)) {
00270 LogDebug("CaloSim") << value;
00271 const std::vector<double> & fvec = value.doubles();
00272 int nval = fvec.size();
00273 if (nval < 1) {
00274 edm::LogError("CaloSim") << "CaloTrkProcessing: # of " << str
00275 << " bins " << nval << " < 1 ==> illegal ";
00276 throw cms::Exception("Unknown", "CaloTrkProcessing")
00277 << "nval < 2 for array " << str << "\n";
00278 }
00279
00280 return fvec;
00281 } else {
00282 edm::LogError("CaloSim") << "CaloTrkProcessing: cannot get array " << str ;
00283 throw cms::Exception("Unknown", "CaloTrkProcessing")
00284 << "cannot get array " << str << "\n";
00285 }
00286 }
00287
00288 int CaloTrkProcessing::isItCalo(const G4VTouchable* touch) {
00289
00290 int lastLevel = -1;
00291 G4LogicalVolume* lv=0;
00292 for (unsigned int it=0; it < detectors.size(); it++) {
00293 if (lastLevel != detectors[it].level) {
00294 lastLevel = detectors[it].level;
00295 lv = detLV(touch, lastLevel);
00296 #ifdef ddebug
00297 std::string name1 = "Unknown";
00298 if (lv != 0) name1 = lv->GetName();
00299 LogDebug("CaloSim") << "CaloTrkProcessing: volume " << name1
00300 << " at Level " << lastLevel;
00301 int levels = detLevels(touch);
00302 if (levels > 0) {
00303 G4String name2[20]; int copyno2[20];
00304 detectorLevel(touch, levels, copyno2, name2);
00305 for (int i2=0; i2<levels; i2++)
00306 LogDebug("CaloSim") << " " << i2 << " " << name2[i2] << " "
00307 << copyno2[i2];
00308 }
00309 #endif
00310 }
00311 bool ok = (lv == detectors[it].lv);
00312 if (ok) return it;
00313 }
00314 return -1;
00315 }
00316
00317 int CaloTrkProcessing::isItInside(const G4VTouchable* touch, int idcal,
00318 int idin) {
00319 int lastLevel = -1;
00320 G4LogicalVolume* lv=0;
00321 int id1, id2;
00322 if (idcal < 0) {id1 = 0; id2 = static_cast<int>(detectors.size());}
00323 else {id1 = idcal; id2 = id1+1;}
00324 for (int it1 = id1; it1 < id2; it1++) {
00325 if (idin < 0) {
00326 for (unsigned int it2 = 0; it2 < detectors[it1].fromDets.size(); it2++) {
00327 if (lastLevel != detectors[it1].fromLevels[it2]) {
00328 lastLevel = detectors[it1].fromLevels[it2];
00329 lv = detLV(touch,lastLevel);
00330 #ifdef ddebug
00331 std::string name1 = "Unknown";
00332 if (lv != 0) name1 = lv->GetName();
00333 LogDebug("CaloSim") << "CaloTrkProcessing: volume " << name1
00334 << " at Level " << lastLevel;
00335 int levels = detLevels(touch);
00336 if (levels > 0) {
00337 G4String name2[20]; int copyno2[20];
00338 detectorLevel(touch, levels, copyno2, name2);
00339 for (int i2=0; i2<levels; i2++)
00340 LogDebug("CaloSim") << " " << i2 << " " << name2[i2] << " "
00341 << copyno2[i2];
00342 }
00343 #endif
00344 }
00345 bool ok = (lv == detectors[it1].fromDetL[it2]);
00346 if (ok) return it2;
00347 }
00348 } else {
00349 lastLevel = detectors[it1].fromLevels[idin];
00350 lv = detLV(touch,lastLevel);
00351 #ifdef ddebug
00352 std::string name1 = "Unknown";
00353 if (lv != 0) name1 = lv->GetName();
00354 LogDebug("CaloSim") << "CaloTrkProcessing: volume " << name1
00355 << " at Level " << lastLevel;
00356 int levels = detLevels(touch);
00357 if (levels > 0) {
00358 G4String name2[20]; int copyno2[20];
00359 detectorLevel(touch, levels, copyno2, name2);
00360 for (int i2=0; i2<levels; i2++)
00361 LogDebug("CaloSim") << " " << i2 << " " << name2[i2] << " "
00362 << copyno2[i2];
00363 }
00364 #endif
00365 bool ok = (lv == detectors[it1].fromDetL[idin]);
00366 if (ok) return idin;
00367 }
00368 }
00369 return -1;
00370 }
00371
00372 int CaloTrkProcessing::detLevels(const G4VTouchable* touch) const {
00373
00374
00375 if (touch)
00376 return ((touch->GetHistoryDepth())+1);
00377 else
00378 return 0;
00379 }
00380
00381 G4LogicalVolume* CaloTrkProcessing::detLV(const G4VTouchable* touch,
00382 int currentlevel) const {
00383
00384 G4LogicalVolume* lv=0;
00385 if (touch) {
00386 int level = ((touch->GetHistoryDepth())+1);
00387 if (level > 0 && level >= currentlevel) {
00388 int ii = level - currentlevel;
00389 lv = touch->GetVolume(ii)->GetLogicalVolume();
00390 return lv;
00391 }
00392 }
00393 return lv;
00394 }
00395
00396 void CaloTrkProcessing::detectorLevel(const G4VTouchable* touch, int& level,
00397 int* copyno, G4String* name) const {
00398
00399 static const std::string unknown("Unknown");
00400
00401 if (level > 0) {
00402 for (int ii = 0; ii < level; ii++) {
00403 int i = level - ii - 1;
00404 G4VPhysicalVolume* pv = touch->GetVolume(i);
00405 if (pv != 0)
00406 name[ii] = pv->GetName();
00407 else
00408 name[ii] = unknown;
00409 copyno[ii] = touch->GetReplicaNumber(i);
00410 }
00411 }
00412 }