CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/SimG4CMS/Forward/src/ZdcSD.cc

Go to the documentation of this file.
00001 
00002 // File: ZdcSD.cc
00003 // Date: 03.01
00004 // Description: Sensitive Detector class for Zdc
00005 // Modifications: 
00007 #include "SimG4CMS/Forward/interface/ZdcSD.h"
00008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00009 #include "G4SDManager.hh"
00010 #include "G4Step.hh"
00011 #include "G4Track.hh"
00012 #include "G4VProcess.hh"
00013  #include "SimG4Core/Notification/interface/TrackInformation.h"
00014 #include "G4ios.hh"
00015 #include "G4Cerenkov.hh"
00016 #include "G4ParticleTable.hh"
00017 #include "CLHEP/Units/GlobalSystemOfUnits.h"
00018 #include "Randomize.hh"
00019 #include "G4Poisson.hh"
00020 
00021 ZdcSD::ZdcSD(G4String name, const DDCompactView & cpv,
00022              SensitiveDetectorCatalog & clg, 
00023              edm::ParameterSet const & p,const SimTrackManager* manager) : 
00024   CaloSD(name, cpv, clg, p, manager), numberingScheme(0) {
00025   edm::ParameterSet m_ZdcSD = p.getParameter<edm::ParameterSet>("ZdcSD");
00026   useShowerLibrary = m_ZdcSD.getParameter<bool>("UseShowerLibrary");
00027   useShowerHits    = m_ZdcSD.getParameter<bool>("UseShowerHits");
00028   zdcHitEnergyCut  = m_ZdcSD.getParameter<double>("ZdcHitEnergyCut")*GeV;
00029   verbosity  = m_ZdcSD.getParameter<int>("Verbosity");
00030   int verbn  = verbosity/10;
00031   verbosity %= 10;
00032   ZdcNumberingScheme* scheme;
00033   scheme = new ZdcNumberingScheme(verbn);
00034   setNumberingScheme(scheme);
00035   
00036   edm::LogInfo("ForwardSim")
00037     << "***************************************************\n"
00038     << "*                                                 *\n"
00039     << "* Constructing a ZdcSD  with name " << name <<"   *\n"
00040     << "*                                                 *\n"
00041     << "***************************************************";
00042 
00043   edm::LogInfo("ForwardSim")
00044      << "\nUse of shower library is set to " 
00045      << useShowerLibrary 
00046      << "\nUse of Shower hits method is set to "
00047      << useShowerHits;                  
00048  
00049   edm::LogInfo("ForwardSim")
00050      << "\nEnergy Threshold Cut set to " 
00051      << zdcHitEnergyCut/GeV
00052      <<" (GeV)";
00053   
00054   if(useShowerLibrary){
00055     showerLibrary = new ZdcShowerLibrary(name, cpv, p);
00056   }
00057 }
00058 
00059 ZdcSD::~ZdcSD() {
00060   
00061   if(numberingScheme) delete numberingScheme;
00062   if(showerLibrary)delete showerLibrary;
00063 
00064   edm::LogInfo("ForwardSim") 
00065     <<"end of ZdcSD\n";
00066 }
00067 
00068 void ZdcSD::initRun(){
00069   if(useShowerLibrary){
00070     G4ParticleTable *theParticleTable = G4ParticleTable::GetParticleTable();
00071     showerLibrary->initRun(theParticleTable); 
00072   }
00073   hits.clear();  
00074 }
00075 
00076 bool ZdcSD::ProcessHits(G4Step * aStep, G4TouchableHistory * ) {
00077 
00078   NaNTrap( aStep ) ;
00079 
00080   if (aStep == NULL) {
00081     return true;
00082   } else {
00083     if(useShowerLibrary){
00084       getFromLibrary(aStep);
00085     }
00086     if(useShowerHits){
00087       if (getStepInfo(aStep)) {
00088         if (hitExists() == false && edepositEM+edepositHAD>0.)
00089           currentHit = CaloSD::createNewHit();
00090       }
00091     }
00092   }
00093   return true;
00094 }
00095 
00096 void ZdcSD::getFromLibrary (G4Step* aStep) {
00097   bool ok = true;
00098 
00099   preStepPoint  = aStep->GetPreStepPoint(); 
00100   theTrack      = aStep->GetTrack();   
00101 
00102   double etrack    = preStepPoint->GetKineticEnergy();
00103   int primaryID = setTrackID(aStep);  
00104 
00105   hits.clear();
00106 
00107   /*
00108     if (etrack >= zdcHitEnergyCut) {
00109     primaryID    = theTrack->GetTrackID();
00110     } else {
00111     primaryID    = theTrack->GetParentID();
00112     if (primaryID == 0) primaryID = theTrack->GetTrackID();
00113     }
00114   */
00115     
00116   // Reset entry point for new primary
00117   posGlobal = preStepPoint->GetPosition();
00118   resetForNewPrimary(posGlobal, etrack);
00119 
00120   if (etrack >= zdcHitEnergyCut){
00121     // create hits only if above threshold
00122 
00123     LogDebug("ForwardSim")
00124       //std::cout
00125       <<"----------------New track------------------------------\n"
00126       <<"Incident EnergyTrack: "<<etrack<< " MeV \n"
00127       <<"Zdc Cut Energy for Hits: "<<zdcHitEnergyCut<<" MeV \n"
00128       << "ZdcSD::getFromLibrary " <<hits.size() <<" hits for "
00129       << GetName() << " of " << primaryID << " with " 
00130       << theTrack->GetDefinition()->GetParticleName() << " of " 
00131       << preStepPoint->GetKineticEnergy()<< " MeV\n"; 
00132     
00133     hits.swap(showerLibrary->getHits(aStep, ok));    
00134   }
00135  
00136   entrancePoint = preStepPoint->GetPosition();
00137   for (unsigned int i=0; i<hits.size(); i++) {
00138     posGlobal           = hits[i].position;
00139     entranceLocal       = hits[i].entryLocal;
00140     double time         = hits[i].time;
00141     unsigned int unitID = hits[i].detID;
00142     edepositHAD         = hits[i].DeHad;
00143     edepositEM          = hits[i].DeEM;
00144     currentID.setID(unitID, time, primaryID);
00145       
00146     // check if it is in the same unit and timeslice as the previous on    
00147     if (currentID == previousID) {
00148       updateHit(currentHit);    
00149     } else {
00150       currentHit = createNewHit();
00151     }
00152       
00153     //  currentHit->setPosition(hitPoint.x(),hitPoint.y(),hitPoint.z());
00154     //  currentHit->setEM(eEM);
00155     //   currentHit->setHadr(eHAD);
00156     currentHit->setIncidentEnergy(etrack);
00157     //  currentHit->setEntryLocal(hitEntry.x(),hitEntry.y(),hitEntry.z());
00158       
00159     LogDebug("ForwardSim") << "ZdcSD: Final Hit number:"<<i<<"-->"
00160                            <<"New HitID: "<<currentHit->getUnitID()
00161                            <<" New Hit trackID: "<<currentHit->getTrackID()
00162                            <<" New EM Energy: "<<currentHit->getEM()/GeV
00163                            <<" New HAD Energy: "<<currentHit->getHadr()/GeV
00164                            <<" New HitEntryPoint: "<<currentHit->getEntryLocal()
00165                            <<" New IncidentEnergy: "<<currentHit->getIncidentEnergy()/GeV
00166                            <<" New HitPosition: "<<posGlobal;
00167   }
00168   
00169   //Now kill the current track
00170   if (ok) {
00171     theTrack->SetTrackStatus(fStopAndKill);
00172     G4TrackVector tv = *(aStep->GetSecondary());
00173     for (unsigned int kk=0; kk<tv.size(); kk++) {
00174       if (tv[kk]->GetVolume() == preStepPoint->GetPhysicalVolume())
00175         tv[kk]->SetTrackStatus(fStopAndKill);
00176     }
00177   }
00178 }
00179 
00180 double ZdcSD::getEnergyDeposit(G4Step * aStep, edm::ParameterSet const & p ) {
00181 
00182   float NCherPhot = 0.;
00183   //std::cout<<"I go through here"<<std::endl;
00184 
00185   if (aStep == NULL) {
00186     LogDebug("ForwardSim") << "ZdcSD::  getEnergyDeposit: aStep is NULL!";
00187     return 0;
00188   } else {
00189     // preStepPoint information
00190     G4SteppingControl  stepControlFlag = aStep->GetControlFlag();
00191     G4StepPoint*       preStepPoint = aStep->GetPreStepPoint();
00192     G4VPhysicalVolume* currentPV    = preStepPoint->GetPhysicalVolume();
00193     G4String           nameVolume   = currentPV->GetName();
00194 
00195     G4ThreeVector      hitPoint = preStepPoint->GetPosition();  
00196     G4ThreeVector      hit_mom = preStepPoint->GetMomentumDirection();
00197     G4double           stepL = aStep->GetStepLength()/cm;
00198     G4double           beta     = preStepPoint->GetBeta();
00199     G4double           charge   = preStepPoint->GetCharge();
00200 
00201     // G4VProcess*        curprocess   = preStepPoint->GetProcessDefinedStep();
00202     // G4String           namePr   = preStepPoint->GetProcessDefinedStep()->GetProcessName();
00203     // G4LogicalVolume*   lv    = currentPV->GetLogicalVolume();
00204     // G4Material*        mat   = lv->GetMaterial();
00205     // G4double           rad   = mat->GetRadlen();
00206 
00207     // postStepPoint information
00208     G4StepPoint* postStepPoint = aStep->GetPostStepPoint();   
00209     G4VPhysicalVolume* postPV = postStepPoint->GetPhysicalVolume();
00210     G4String postnameVolume = postPV->GetName();
00211 
00212     // theTrack information
00213     G4Track* theTrack = aStep->GetTrack();   
00214     G4String particleType = theTrack->GetDefinition()->GetParticleName();
00215     G4int primaryID = theTrack->GetTrackID();
00216     G4double entot = theTrack->GetTotalEnergy();
00217     G4ThreeVector vert_mom = theTrack->GetVertexMomentumDirection();
00218     G4ThreeVector localPoint = theTrack->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
00219 
00220     // calculations
00221     float costheta = vert_mom.z()/sqrt(vert_mom.x()*vert_mom.x()+
00222                                        vert_mom.y()*vert_mom.y()+
00223                                        vert_mom.z()*vert_mom.z());
00224     float theta = acos(std::min(std::max(costheta,float(-1.)),float(1.)));
00225     float eta = -log(tan(theta/2));
00226     float phi = -100.;
00227     if (vert_mom.x() != 0) phi = atan2(vert_mom.y(),vert_mom.x()); 
00228     if (phi < 0.) phi += twopi;
00229 
00230     // Get the total energy deposit
00231     double stepE   = aStep->GetTotalEnergyDeposit();
00232     LogDebug("ForwardSim") 
00233       << "ZdcSD::  getEnergyDeposit: "
00234       <<"*****************HHHHHHHHHHHHHHHHHHHHHHHHHHLLLLLLLLLlllllllllll&&&&&&&&&&\n"
00235       << "  preStepPoint: " << nameVolume << "," << stepL << "," << stepE 
00236       << "," << beta << "," << charge << "\n"
00237       << "  postStepPoint: " << postnameVolume << "," << costheta << "," 
00238       << theta << "," << eta << "," << phi << "," << particleType << "," 
00239       << primaryID;
00240 
00241     float bThreshold = 0.67;
00242     if ((beta > bThreshold) && (charge != 0) && (nameVolume == "ZDC_EMFiber" || nameVolume == "ZDC_HadFiber")) {
00243       LogDebug("ForwardSim") << "ZdcSD::  getEnergyDeposit:  pass "; 
00244 
00245       float nMedium = 1.4925;
00246       // float photEnSpectrDL = 10714.285714;
00247       //       photEnSpectrDL = (1./400.nm-1./700.nm)*10000000.cm/nm; /* cm-1  */
00248 
00249       float photEnSpectrDE = 1.24;
00250       // E = 2pi*(1./137.)*(eV*cm/370.)/lambda = 12.389184*(eV*cm)/lambda
00251       // Emax = 12.389184*(eV*cm)/400nm*10-7cm/nm  = 3.01 eV
00252       // Emin = 12.389184*(eV*cm)/700nm*10-7cm/nm  = 1.77 eV
00253       // delE = Emax - Emin = 1.24 eV
00254 
00255       float effPMTandTransport = 0.15;
00256 
00257       // Check these values
00258       float thFullRefl = 23.;
00259       float thFullReflRad = thFullRefl*pi/180.;
00260 
00261       edm::ParameterSet m_ZdcSD = p.getParameter<edm::ParameterSet>("ZdcSD");
00262       thFibDir  = m_ZdcSD.getParameter<double>("FiberDirection");
00263       //float thFibDir = 90.;
00264       float thFibDirRad = thFibDir*pi/180.;
00265 
00266       // at which theta the point is located:
00267       //   float th1 = hitPoint.theta();
00268 
00269       // theta of charged particle in LabRF(hit momentum direction):
00270       float costh = hit_mom.z()/sqrt(hit_mom.x()*hit_mom.x()+
00271                                      hit_mom.y()*hit_mom.y()+
00272                                      hit_mom.z()*hit_mom.z());
00273       float th = acos(std::min(std::max(costh,float(-1.)),float(1.)));
00274       // just in case (can do both standard ranges of phi):
00275       if (th < 0.) th += twopi;
00276 
00277       // theta of cone with Cherenkov photons w.r.t.direction of charged part.:
00278       float costhcher =1./(nMedium*beta);
00279       float thcher = acos(std::min(std::max(costhcher,float(-1.)),float(1.)));
00280 
00281       // diff thetas of charged part. and quartz direction in LabRF:
00282       float DelFibPart = fabs(th - thFibDirRad);
00283 
00284       // define real distances:
00285       float d = fabs(tan(th)-tan(thFibDirRad));   
00286 
00287       // float a = fabs(tan(thFibDirRad)-tan(thFibDirRad+thFullReflRad));   
00288       // float r = fabs(tan(th)-tan(th+thcher));   
00289       float a = tan(thFibDirRad)+tan(fabs(thFibDirRad-thFullReflRad));   
00290       float r = tan(th)+tan(fabs(th-thcher));   
00291       
00292       // std::cout.testOut << "  d=|tan(" << th << ")-tan(" << thFibDirRad << ")| "
00293       //              << "=|" << tan(th) << "-" << tan(thFibDirRad) << "| = " << d;
00294       // std::cout.testOut << "  a=tan(" << thFibDirRad << ")=" << tan(thFibDirRad) 
00295       //              << " + tan(|" << thFibDirRad << " - " << thFullReflRad << "|)="
00296       //              << tan(fabs(thFibDirRad-thFullReflRad)) << " = " << a;
00297       // std::cout.testOut << "  r=tan(" << th << ")=" << tan(th) << " + tan(|" << th 
00298       //              << " - " << thcher << "|)=" << tan(fabs(th-thcher)) << " = " << r;
00299 
00300       // define losses d_qz in cone of full reflection inside quartz direction
00301       float d_qz = -1;
00302       float variant = -1;
00303 
00304       // if (d > (r+a))
00305       if (DelFibPart > (thFullReflRad + thcher) ) {
00306         variant = 0.; d_qz = 0.;
00307       } else {
00308         // if ((DelFibPart + thcher) < thFullReflRad )  [(d+r) < a]
00309         if ((th + thcher) < (thFibDirRad+thFullReflRad) && (th - thcher) > (thFibDirRad-thFullReflRad) ) {
00310           variant = 1.; d_qz = 1.;
00311         } else {
00312           // if ((thcher - DelFibPart ) > thFullReflRad )  [(r-d) > a]
00313           if ((thFibDirRad + thFullReflRad) < (th + thcher) && (thFibDirRad - thFullReflRad) > (th - thcher) ) {
00314             variant = 2.; d_qz = 0.;
00315           } else {
00316             // if ((thcher + DelFibPart ) > thFullReflRad && thcher < (DelFibPart+thFullReflRad) ) {  [(r+d) > a && (r-d) < a)]
00317             variant = 3.; // d_qz is calculated below
00318 
00319             // use crossed length of circles(cone projection) - dC1/dC2 : 
00320             float arg_arcos = 0.;
00321             float tan_arcos = 2.*a*d;
00322             if (tan_arcos != 0.) arg_arcos =(r*r-a*a-d*d)/tan_arcos; 
00323             // std::cout.testOut << "  d_qz: " << r << "," << a << "," << d << " " << tan_arcos << " " << arg_arcos;
00324             arg_arcos = fabs(arg_arcos);
00325             // std::cout.testOut << "," << arg_arcos;
00326             float th_arcos = acos(std::min(std::max(arg_arcos,float(-1.)),float(1.)));
00327             // std::cout.testOut << " " << th_arcos;
00328             d_qz = th_arcos/pi/2.;
00329             // std::cout.testOut << " " << d_qz;
00330             d_qz = fabs(d_qz);
00331             // std::cout.testOut << "," << d_qz;
00332           }
00333         }
00334       }
00335 
00336       //  std::cout<< std::endl;
00337 
00338       double meanNCherPhot = 0.;
00339       G4int poissNCherPhot = 0;
00340       if (d_qz > 0) {
00341         meanNCherPhot = 370.*charge*charge*( 1. - 1./(nMedium*nMedium*beta*beta) ) * photEnSpectrDE * stepL;
00342 
00343         // dLamdX:  meanNCherPhot = (2.*pi/137.)*charge*charge* 
00344         //                          ( 1. - 1./(nMedium*nMedium*beta*beta) ) * photEnSpectrDL * stepL;
00345         poissNCherPhot = (G4int) G4Poisson(meanNCherPhot);
00346 
00347         if (poissNCherPhot < 0) poissNCherPhot = 0; 
00348 
00349         // NCherPhot = meanNCherPhot;
00350         NCherPhot = poissNCherPhot * effPMTandTransport * d_qz;
00351       }
00352 
00353       LogDebug("ForwardSim") 
00354         << "ZdcSD::  getEnergyDeposit:  gED: "
00355         << stepE
00356         << "," << costh
00357         << "," << th
00358         << "," << costhcher
00359         << "," << thcher
00360         << "," << DelFibPart
00361         << "," << d
00362         << "," << a
00363         << "," << r
00364         << "," << hitPoint
00365         << "," << hit_mom
00366         << "," << stepControlFlag
00367         << "," << entot
00368         << "," << vert_mom
00369         << "," << localPoint
00370         << "," << charge
00371         << "," << beta
00372         << "," << stepL
00373         << "," << d_qz
00374         << "," << variant
00375         << "," << meanNCherPhot
00376         << "," << poissNCherPhot
00377         << "," << NCherPhot;
00378       // --constants-----------------
00379       // << "," << photEnSpectrDE
00380       // << "," << nMedium
00381       // << "," << bThreshold
00382       // << "," << thFibDirRad
00383       // << "," << thFullReflRad
00384       // << "," << effPMTandTransport
00385       // --other variables-----------
00386       // << "," << curprocess
00387       // << "," << nameProcess
00388       // << "," << name
00389       // << "," << rad
00390       // << "," << mat
00391 
00392     } else {
00393       // determine failure mode: beta, charge, and/or nameVolume
00394       if (beta <= bThreshold)
00395         LogDebug("ForwardSim") 
00396           << "ZdcSD::  getEnergyDeposit: fail beta=" << beta;
00397       if (charge == 0)
00398         LogDebug("ForwardSim") 
00399           << "ZdcSD::  getEnergyDeposit: fail charge=0";
00400       if ( !(nameVolume == "ZDC_EMFiber" || nameVolume == "ZDC_HadFiber") )
00401         LogDebug("ForwardSim") 
00402           << "ZdcSD::  getEnergyDeposit: fail nv=" << nameVolume;
00403     }
00404 
00405     return NCherPhot;
00406   } 
00407 }
00408 
00409 uint32_t ZdcSD::setDetUnitId(G4Step* aStep) {
00410   uint32_t returnNumber = 0;
00411   if(numberingScheme != 0)returnNumber = numberingScheme->getUnitID(aStep);
00412   // edm: return (numberingScheme == 0 ? 0 : numberingScheme->getUnitID(aStep));
00413   return returnNumber;
00414 }
00415 
00416 void ZdcSD::setNumberingScheme(ZdcNumberingScheme* scheme) {
00417   if (scheme != 0) {
00418     edm::LogInfo("ForwardSim") << "ZdcSD: updates numbering scheme for " 
00419                                << GetName();
00420     if (numberingScheme) delete numberingScheme;
00421     numberingScheme = scheme;
00422   }
00423 }
00424 
00425 int ZdcSD::setTrackID (G4Step* aStep) {
00426   theTrack     = aStep->GetTrack();
00427   double etrack = preStepPoint->GetKineticEnergy();
00428   TrackInformation * trkInfo = (TrackInformation *)(theTrack->GetUserInformation());
00429   int primaryID = trkInfo->getIDonCaloSurface();
00430   if (primaryID == 0) {
00431 #ifdef DebugLog
00432     LogDebug("ZdcSD") << "ZdcSD: Problem with primaryID **** set by force "
00433                         << "to TkID **** " << theTrack->GetTrackID();
00434 #endif
00435     primaryID = theTrack->GetTrackID();
00436     }
00437   if (primaryID != previousID.trackID())
00438       resetForNewPrimary(preStepPoint->GetPosition(), etrack); 
00439   return primaryID;
00440 }