CMS 3D CMS Logo

CaloTrkProcessing.cc

Go to the documentation of this file.
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 //#define ddebug
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   //Initialise the parameter set
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   //Get the names 
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   // define if you are at the surface of CALO  
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   //Return number of levels
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   //Get name and copy numbers
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 }

Generated on Tue Jun 9 17:46:49 2009 for CMSSW by  doxygen 1.5.4