#include <SimG4CMS/Forward/interface/CastorSD.h>

Inheritance diagram for CastorSD:

Public Member Functions
	CastorSD (G4String, const DDCompactView &, SensitiveDetectorCatalog &clg, edm::ParameterSet const &, const SimTrackManager *)

virtual double	getEnergyDeposit (G4Step *)

virtual uint32_t	setDetUnitId (G4Step *step)

void	setNumberingScheme (CastorNumberingScheme *scheme)

virtual	~CastorSD ()

Public Member Functions inherited from CaloSD
	CaloSD (G4String aSDname, const DDCompactView &cpv, SensitiveDetectorCatalog &clg, edm::ParameterSet const &p, const SimTrackManager *, int tSlice=1, bool ignoreTkID=false)

virtual void	clear ()

virtual void	DrawAll ()

virtual void	EndOfEvent (G4HCofThisEvent *eventHC)

void	fillHits (edm::PCaloHitContainer &, std::string n)

virtual void	Initialize (G4HCofThisEvent *HCE)

virtual void	PrintAll ()

virtual bool	ProcessHits (G4Step step, G4TouchableHistory tHistory)

virtual bool	ProcessHits (G4GFlashSpot aSpot, G4TouchableHistory )

virtual	~CaloSD ()

Public Member Functions inherited from SensitiveCaloDetector
	SensitiveCaloDetector (std::string &iname, const DDCompactView &cpv, SensitiveDetectorCatalog &clg, edm::ParameterSet const &p)

Public Member Functions inherited from SensitiveDetector
virtual void	AssignSD (std::string &vname)

Local3DPoint	ConvertToLocal3DPoint (G4ThreeVector point)

Local3DPoint	FinalStepPosition (G4Step *s, coordinates)

virtual std::vector< std::string >	getNames ()

Local3DPoint	InitialStepPosition (G4Step *s, coordinates)

std::string	nameOfSD ()

void	NaNTrap (G4Step *step)

void	Register ()

	SensitiveDetector (std::string &iname, const DDCompactView &cpv, SensitiveDetectorCatalog &, edm::ParameterSet const &p)

virtual	~SensitiveDetector ()

Public Member Functions inherited from Observer< const BeginOfRun * >
	Observer ()

void	slotForUpdate (const BeginOfRun *iT)

virtual	~Observer ()

Public Member Functions inherited from Observer< const BeginOfEvent * >
	Observer ()

void	slotForUpdate (const BeginOfEvent *iT)

virtual	~Observer ()

Public Member Functions inherited from Observer< const BeginOfTrack * >
	Observer ()

void	slotForUpdate (const BeginOfTrack *iT)

virtual	~Observer ()

Public Member Functions inherited from Observer< const EndOfTrack * >
	Observer ()

void	slotForUpdate (const EndOfTrack *iT)

virtual	~Observer ()

Public Member Functions inherited from Observer< const EndOfEvent * >
	Observer ()

void	slotForUpdate (const EndOfEvent *iT)

virtual	~Observer ()

Protected Member Functions
virtual void	initRun ()

Protected Member Functions inherited from CaloSD
G4bool	checkHit ()

virtual void	clearHits ()

CaloG4Hit *	createNewHit ()

virtual bool	filterHit (CaloG4Hit *, double)

double	getAttenuation (G4Step *aStep, double birk1, double birk2, double birk3)

virtual uint16_t	getDepth (G4Step *)

double	getResponseWt (G4Track *)

virtual G4bool	getStepInfo (G4Step *aStep)

virtual int	getTrackID (G4Track *)

G4bool	hitExists ()

void	resetForNewPrimary (G4ThreeVector, double)

G4ThreeVector	setToLocal (G4ThreeVector, const G4VTouchable *)

virtual void	update (const BeginOfRun *)
	This routine will be called when the appropriate signal arrives. More...

virtual void	update (const BeginOfEvent *)
	This routine will be called when the appropriate signal arrives. More...

virtual void	update (const BeginOfTrack *trk)
	This routine will be called when the appropriate signal arrives. More...

virtual void	update (const EndOfTrack *trk)
	This routine will be called when the appropriate signal arrives. More...

virtual void	update (const ::EndOfEvent *)

void	updateHit (CaloG4Hit *)

Protected Member Functions inherited from Observer< const EndOfEvent * >
virtual void	update (const EndOfEvent *)=0
	This routine will be called when the appropriate signal arrives. More...

Private Member Functions
void	getFromLibrary (G4Step *)

uint32_t	rotateUnitID (uint32_t, G4Track *, CastorShowerEvent)

int	setTrackID (G4Step *)

Private Attributes
double	energyThresholdSL

G4LogicalVolume *	lvC3EF

G4LogicalVolume *	lvC3HF

G4LogicalVolume *	lvC4EF

G4LogicalVolume *	lvC4HF

G4LogicalVolume *	lvCAST

CastorNumberingScheme *	numberingScheme

CastorShowerLibrary *	showerLibrary

bool	useShowerLibrary

Additional Inherited Members
Public Types inherited from SensitiveDetector
enum	coordinates { WorldCoordinates, LocalCoordinates }

Protected Attributes inherited from CaloSD
int	checkHits

double	correctT

bool	corrTOFBeam

CaloG4Hit *	currentHit

CaloHitID	currentID

float	edepositEM

float	edepositHAD

double	eminHit

double	eminHitD

G4int	emPDG

double	energyCut

G4ThreeVector	entranceLocal

G4ThreeVector	entrancePoint

G4int	epPDG

bool	forceSave

G4int	gammaPDG

float	incidentEnergy

double	kmaxIon

double	kmaxNeutron

double	kmaxProton

const SimTrackManager *	m_trackManager

G4ThreeVector	posGlobal

G4StepPoint *	preStepPoint

CaloHitID	previousID

int	primIDSaved

bool	runInit

bool	suppressHeavy

G4Track *	theTrack

double	tmaxHit

bool	useMap

Detailed Description

Description: Stores hits of Castor in appropriate container

Usage: Used in sensitive detector builder

Definition at line 30 of file CastorSD.h.

Constructor & Destructor Documentation

CastorSD::CastorSD	(	G4String	name,
		const DDCompactView &	cpv,
		SensitiveDetectorCatalog &	clg,
		edm::ParameterSet const &	p,
		const SimTrackManager *	manager
	)

Definition at line 30 of file CastorSD.cc.

References energyThresholdSL, edm::ParameterSet::getParameter(), lvC3EF, lvC3HF, lvC4EF, lvC4HF, lvCAST, setNumberingScheme(), showerLibrary, and useShowerLibrary.

                                            : 
   CaloSD(name, cpv, clg, p, manager), numberingScheme(0), lvC3EF(0),
   lvC3HF(0), lvC4EF(0), lvC4HF(0), lvCAST(0) {
   
   edm::ParameterSet m_CastorSD = p.getParameter<edm::ParameterSet>("CastorSD");
   useShowerLibrary  = m_CastorSD.getParameter<bool>("useShowerLibrary");
   energyThresholdSL = m_CastorSD.getParameter<double>("minEnergyInGeVforUsingSLibrary");
   energyThresholdSL = energyThresholdSL*GeV;   //  Convert GeV => MeV 
   
   if (useShowerLibrary) showerLibrary = new CastorShowerLibrary(name, p);
   
   setNumberingScheme(new CastorNumberingScheme());
   
   edm::LogInfo("ForwardSim") 
     << "***************************************************\n"
     << "*                                                 *\n" 
     << "* Constructing a CastorSD  with name " << GetName() << "\n"
     << "*                                                 *\n"
     << "***************************************************";
     
   const G4LogicalVolumeStore * lvs = G4LogicalVolumeStore::GetInstance();
   std::vector<G4LogicalVolume*>::const_iterator lvcite;
   for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
     if (strcmp(((*lvcite)->GetName()).c_str(),"C3EF") == 0) lvC3EF = (*lvcite);
     if (strcmp(((*lvcite)->GetName()).c_str(),"C3HF") == 0) lvC3HF = (*lvcite);
     if (strcmp(((*lvcite)->GetName()).c_str(),"C4EF") == 0) lvC4EF = (*lvcite);
     if (strcmp(((*lvcite)->GetName()).c_str(),"C4HF") == 0) lvC4HF = (*lvcite);
     if (strcmp(((*lvcite)->GetName()).c_str(),"CAST") == 0) lvCAST = (*lvcite);
     if (lvC3EF != 0 && lvC3HF != 0 && lvC4EF != 0 && lvC4HF != 0 && lvCAST != 0) break;
   }
   edm::LogInfo("ForwardSim") << "CastorSD:: LogicalVolume pointers\n"
                  << lvC3EF << " for C3EF; " << lvC3HF 
                  << " for C3HF; " << lvC4EF << " for C4EF; " 
                  << lvC4HF << " for C4HF; " 
                  << lvCAST << " for CAST. " << std::endl;
 
   //  if(useShowerLibrary) edm::LogInfo("ForwardSim") << "\n Using Castor Shower Library \n";
 
 }

CastorSD::~CastorSD ( )

virtual

Definition at line 75 of file CastorSD.cc.

References showerLibrary, and useShowerLibrary.

                     {
   if (useShowerLibrary) delete showerLibrary;
 }

Member Function Documentation

double CastorSD::getEnergyDeposit ( G4Step * aStep )

virtual

Reimplemented from CaloSD.

Definition at line 92 of file CastorSD.cc.

References a, abs, beta, DeDxDiscriminatorTools::charge(), dot(), energyThresholdSL, eta(), getFromLibrary(), funct::log(), LogDebug, lvC3EF, lvC3HF, lvC4EF, lvC4HF, lvCAST, M_PI, max(), min, SensitiveDetector::name, NULL, ExpressReco_HICollisions_FallBack::particleType, phi, pi, CaloSD::preStepPoint, csvReporter::r, dttmaxenums::R, mathSSE::sqrt(), funct::tan(), theta(), CaloSD::theTrack, and useShowerLibrary.

                                                 {
   
   float NCherPhot = 0.;
 
   if (aStep == NULL) {
     return 0;
   } else {
     // preStepPoint information *********************************************
   
     G4StepPoint*       preStepPoint = aStep->GetPreStepPoint();
     G4VPhysicalVolume* currentPV    = preStepPoint->GetPhysicalVolume();
     G4LogicalVolume*   currentLV    = currentPV->GetLogicalVolume();
     G4String           name   = currentPV->GetName();
     std::string        nameVolume;
     nameVolume.assign(name,0,4);
     
 #ifdef debugLog
     G4SteppingControl  stepControlFlag = aStep->GetControlFlag();
     if (aStep->IsFirstStepInVolume()) 
       LogDebug("ForwardSim") << "CastorSD::getEnergyDeposit:"
                  << "\n IsFirstStepInVolume " ; 
 #endif
     
     // Get theTrack 
     G4Track*        theTrack = aStep->GetTrack();
     
 #ifdef debugLog
     if (useShowerLibrary && currentLV==lvCAST) {
           LogDebug("ForwardSim") << "CastorSD::getEnergyDeposit:"
                  << "\n TrackID , ParentID , ParticleName ,"
                  << " eta , phi , z , time ,"
                  << " K , E , Mom " 
                  << "\n  TRACKINFO: " 
                  << theTrack->GetTrackID() 
                  << " , " 
                  << theTrack->GetParentID() 
                  << " , "
                  << theTrack->GetDefinition()->GetParticleName() 
                  << " , "
                  << theTrack->GetPosition().eta() 
                  << " , "
                  << theTrack->GetPosition().phi() 
                  << " , "
                  << theTrack->GetPosition().z() 
                  << " , "
                  << theTrack->GetGlobalTime() 
                  << " , "
                  << theTrack->GetKineticEnergy() 
                  << " , "
                  << theTrack->GetTotalEnergy() 
                  << " , "
                                  << theTrack->GetMomentum().mag() ;
       if(theTrack->GetTrackID() != 1) 
         LogDebug("ForwardSim") << "CastorSD::getEnergyDeposit:"
                    << "\n CurrentStepNumber , TrackID , Particle , VertexPosition ,"
                    << " LogicalVolumeAtVertex , CreatorProcess"
                    << "\n  TRACKINFO2: " 
                    << theTrack->GetCurrentStepNumber() 
                    << " , " 
                  << theTrack->GetTrackID() 
                  << " , "
                  << theTrack->GetDefinition()->GetParticleName() 
                  << " , "
                    << theTrack->GetVertexPosition() 
                    << " , "
                    << theTrack->GetLogicalVolumeAtVertex()->GetName() 
                    << " , " 
                    << theTrack->GetCreatorProcess()->GetProcessName() ;
     } // end of if(useShowerLibrary)
 #endif
     
     // if particle moves from interaction point or "backwards (halo)
     bool backward = false;
     G4ThreeVector  hitPoint = preStepPoint->GetPosition();  
     G4ThreeVector  hit_mom  = preStepPoint->GetMomentumDirection();
     double zint = hitPoint.z();
     double pz   = hit_mom.z();
 
     // Check if theTrack moves backward 
     if (pz * zint < 0.) backward = true;
     
     // Check that theTrack is above the energy threshold to use Shower Library 
     bool aboveThreshold = false;
     if(theTrack->GetKineticEnergy() > energyThresholdSL) aboveThreshold = true;
     
     // Check if theTrack is a muon (if so, DO NOT use Shower Library) 
     bool notaMuon = true;
     G4int mumPDG  =  13;
     G4int mupPDG  = -13;
     G4int parCode = theTrack->GetDefinition()->GetPDGEncoding();
     if (parCode == mupPDG || parCode == mumPDG ) notaMuon = false;
     
     // angle condition
     double theta_max = M_PI - 3.1305; // angle in radians corresponding to -5.2 eta
     double R_mom=sqrt(hit_mom.x()*hit_mom.x() + hit_mom.y()*hit_mom.y());
     double theta = atan2(R_mom,std::abs(pz));
     bool angleok = false;
     if ( theta < theta_max) angleok = true;
     
     // OkToUse
     double R = sqrt(hitPoint.x()*hitPoint.x() + hitPoint.y()*hitPoint.y());
     bool dot = false;
     if ( zint < -14450. && R < 45.) dot = true;
     bool inRange = true;
     if ( zint < -14700. || R > 193.) inRange = false;
     bool OkToUse = false;
     if ( inRange && !dot) OkToUse = true;
     
     if (useShowerLibrary && aboveThreshold && notaMuon && (!backward) && OkToUse && angleok && currentLV == lvCAST ) {
       // Use Castor shower library if energy is above threshold, is not a muon 
       // and is not moving backward 
       getFromLibrary(aStep);
     
 #ifdef debugLog
       LogDebug("ForwardSim") << " Current logical volume is " << nameVolume ;
 #endif
       return NCherPhot;
     } else {
       // Usual calculations
       // G4ThreeVector      hitPoint = preStepPoint->GetPosition(); 
       // G4ThreeVector      hit_mom = preStepPoint->GetMomentumDirection();
       G4double           stepl    = aStep->GetStepLength()/cm;
       G4double           beta     = preStepPoint->GetBeta();
       G4double           charge   = preStepPoint->GetCharge();
       //        G4VProcess*        curprocess   = preStepPoint->GetProcessDefinedStep();
       //        G4String           namePr   = preStepPoint->GetProcessDefinedStep()->GetProcessName();
       //        std::string nameProcess;
       //        nameProcess.assign(namePr,0,4);
 
       //        G4LogicalVolume*   lv    = currentPV->GetLogicalVolume();
       //        G4Material*        mat   = lv->GetMaterial();
       //        G4double           rad   = mat->GetRadlen();
 
 
       // postStepPoint information *********************************************
       G4StepPoint* postStepPoint= aStep->GetPostStepPoint();   
       G4VPhysicalVolume* postPV    = postStepPoint->GetPhysicalVolume();
       G4String           postname   = postPV->GetName();
       std::string        postnameVolume;
       postnameVolume.assign(postname,0,4);
 
       // theTrack information  *************************************************
       // G4Track*        theTrack = aStep->GetTrack();   
       //G4double        entot    = theTrack->GetTotalEnergy();
       G4ThreeVector   vert_mom = theTrack->GetVertexMomentumDirection();
 
       G4ThreeVector  localPoint = theTrack->GetTouchable()->GetHistory()->
     GetTopTransform().TransformPoint(hitPoint);
 
       // calculations...       *************************************************
       float phi = -100.;
       if (vert_mom.x() != 0) phi = atan2(vert_mom.y(),vert_mom.x()); 
       if (phi < 0.) phi += twopi;
       G4String       particleType = theTrack->GetDefinition()->GetParticleName();
 #ifdef debugLog
       float costheta =vert_mom.z()/sqrt(vert_mom.x()*vert_mom.x()+
                     vert_mom.y()*vert_mom.y()+
                     vert_mom.z()*vert_mom.z());
       float theta = acos(std::min(std::max(costheta,float(-1.)),float(1.)));
       float eta = -log(tan(theta/2));
       G4int          primaryID    = theTrack->GetTrackID();
       // *************************************************
 
 
       // *************************************************
       double edep   = aStep->GetTotalEnergyDeposit();
 #endif
 
       // *************************************************
 
 
       // *************************************************
       // take into account light collection curve for plate
       //      double weight = curve_Castor(nameVolume, preStepPoint);
       //      double edep   = aStep->GetTotalEnergyDeposit() * weight;
       // *************************************************
 
 
       // *************************************************
       /*    comments for sensitive volumes:      
         C001 ...-... CP06,CBU1 ...-...CALM --- > fibres and bundle 
         for first release of CASTOR
         CASF  --- > quartz plate  for first and second releases of CASTOR  
         GF2Q, GFNQ, GR2Q, GRNQ 
         for tests with my own test geometry of HF (on ask of Gavrilov)
         C3TF, C4TF - for third release of CASTOR
       */
       double meanNCherPhot=0;
 
       if (currentLV == lvC3EF || currentLV == lvC4EF || currentLV == lvC3HF ||
       currentLV == lvC4HF) {
     //      if(nameVolume == "C3EF" || nameVolume == "C4EF" || nameVolume == "C3HF" || nameVolume == "C4HF") {
 
     float bThreshold = 0.67;
     float nMedium = 1.4925;
     //     float photEnSpectrDL = (1./400.nm-1./700.nm)*10000000.cm/nm; /* cm-1  */
     //     float photEnSpectrDL = 10714.285714;
 
     float photEnSpectrDE = 1.24;    /* see below   */
     /*     E = 2pi*(1./137.)*(eV*cm/370.)/lambda  =     */
     /*       = 12.389184*(eV*cm)/lambda                 */
     /*     Emax = 12.389184*(eV*cm)/400nm*10-7cm/nm  = 3.01 eV     */
     /*     Emin = 12.389184*(eV*cm)/700nm*10-7cm/nm  = 1.77 eV     */
     /*     delE = Emax - Emin = 1.24 eV  --> */
     /*   */
     /* default for Castor nameVolume  == "CASF" or (C3TF & C4TF)  */
     float effPMTandTransport = 0.19;
     /* for test HF geometry volumes:   
        if(nameVolume == "GF2Q" || nameVolume == "GFNQ" || 
        nameVolume == "GR2Q" || nameVolume == "GRNQ")
        effPMTandTransport = 0.19;
     */
 
     float thFullRefl = 23.;  /* 23.dergee */
     float thFullReflRad = thFullRefl*pi/180.;
 
     /* default for Castor nameVolume  == "CASF" or (C3TF & C4TF)  */
     float thFibDir = 45.;  /* .dergee */
     /* for test HF geometry volumes:   
        if(nameVolume == "GF2Q" || nameVolume == "GFNQ" ||
        nameVolume == "GR2Q" || nameVolume == "GRNQ")
        thFibDir = 0.0; // .dergee
     */
     float thFibDirRad = thFibDir*pi/180.;
     /*   */
     /*   */
 
     // at which theta the point is located:
     //     float th1    = hitPoint.theta();
 
     // theta of charged particle in LabRF(hit momentum direction):
     float costh =hit_mom.z()/sqrt(hit_mom.x()*hit_mom.x()+
                       hit_mom.y()*hit_mom.y()+
                       hit_mom.z()*hit_mom.z());
     if (zint < 0) costh = -costh;
     float th = acos(std::min(std::max(costh,float(-1.)),float(1.)));
 
     // just in case (can do bot use):
     if (th < 0.) th += twopi;
 
 
 
     // theta of cone with Cherenkov photons w.r.t.direction of charged part.:
     float costhcher =1./(nMedium*beta);
     float thcher = acos(std::min(std::max(costhcher,float(-1.)),float(1.)));
 
     // diff thetas of charged part. and quartz direction in LabRF:
     float DelFibPart = fabs(th - thFibDirRad);
 
     // define real distances:
     float d = fabs(tan(th)-tan(thFibDirRad));   
 
     //       float a = fabs(tan(thFibDirRad)-tan(thFibDirRad+thFullReflRad));   
     //       float r = fabs(tan(th)-tan(th+thcher));   
     
     float a = tan(thFibDirRad)+tan(fabs(thFibDirRad-thFullReflRad));   
     float r = tan(th)+tan(fabs(th-thcher));   
 
 
     // define losses d_qz in cone of full reflection inside quartz direction
     float d_qz;
     float variant;
 
     if(DelFibPart > (thFullReflRad + thcher) ) {d_qz = 0.; variant=0.;}
     // if(d > (r+a) ) {d_qz = 0.; variant=0.;}
     else {
       if((th + thcher) < (thFibDirRad+thFullReflRad) && 
          (th - thcher) > (thFibDirRad-thFullReflRad) 
          ) {d_qz = 1.; variant=1.;}
       // if((DelFibPart + thcher) < thFullReflRad ) {d_qz = 1.; variant=1.;}
       // if((d+r) < a ) {d_qz = 1.; variant=1.;}
       else {
         if((thFibDirRad + thFullReflRad) < (th + thcher) && 
            (thFibDirRad - thFullReflRad) > (th - thcher) ) {
           // if((thcher - DelFibPart ) > thFullReflRad ) 
           // if((r-d ) > a ) 
           d_qz = 0.; variant=2.;
 
         } else {
           // if((thcher + DelFibPart ) > thFullReflRad && 
           //    thcher < (DelFibPart+thFullReflRad) ) 
           //      {
           d_qz = 0.; variant=3.;
 
 
           // use crossed length of circles(cone projection)
           // dC1/dC2 : 
           float arg_arcos = 0.;
           float tan_arcos = 2.*a*d;
           if(tan_arcos != 0.) arg_arcos =(r*r-a*a-d*d)/tan_arcos; 
           arg_arcos = fabs(arg_arcos);
           float th_arcos = acos(std::min(std::max(arg_arcos,float(-1.)),float(1.)));
           d_qz = th_arcos/pi/2.;
           d_qz = fabs(d_qz);
 
 
 
           //        }
           //             else
           //         {
           //                d_qz = 0.; variant=4.;
           //#ifdef debugLog
           // std::cout <<" ===============>variant 4 information: <===== " <<std::endl;
           // std::cout <<" !!!!!!!!!!!!!!!!!!!!!!  variant = " << variant  <<std::endl;
           //#endif 
           //
           //         }
         }
       }
     }
 
 
     // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
     if(charge != 0. && beta > bThreshold && d_qz != 0. ) {
 
       meanNCherPhot = 370.*charge*charge* 
         ( 1. - 1./(nMedium*nMedium*beta*beta) )*
         photEnSpectrDE*stepl;
 
       // dLamdX:
       //   meanNCherPhot = (2.*pi/137.)*charge*charge* 
       //                     ( 1. - 1./(nMedium*nMedium*beta*beta) )*
       //                     photEnSpectrDL*stepl;
 
 
       //     NCherPhot = meanNCherPhot;
       // Poisson:
       G4int poissNCherPhot = (G4int) G4Poisson(meanNCherPhot);
 
       if(poissNCherPhot < 0) poissNCherPhot = 0; 
       
       NCherPhot = poissNCherPhot * effPMTandTransport * d_qz;
 
 
 
 #ifdef debugLog
       float thgrad = th*180./pi;
       float thchergrad = thcher*180./pi;
       float DelFibPartgrad = DelFibPart*180./pi;
       LogDebug("ForwardSim") << " ==============================> start all "
                  << "information:<========= \n" << " =====> for "
                  << "test:<===  \n" << " variant = " << variant  
                  << "\n thgrad = " << thgrad  << "\n thchergrad "
                  << "= " << thchergrad  << "\n DelFibPartgrad = "
                  << DelFibPartgrad << "\n d_qz = " << d_qz  
                  << "\n =====> Start Step Information <===  \n"
                  << " ===> calo preStepPoint info <===  \n" 
                  << " hitPoint = " << hitPoint  << "\n"
                  << " hitMom = " << hit_mom  << "\n"
                  << " stepControlFlag = " << stepControlFlag 
         // << "\n curprocess = " << curprocess << "\n"
         // << " nameProcess = " << nameProcess 
                  << "\n charge = " << charge << "\n"
                  << " beta = " << beta << "\n"
                  << " bThreshold = " << bThreshold << "\n"
                  << " thgrad =" << thgrad << "\n"
                  << " effPMTandTransport=" << effPMTandTransport 
         // << "\n volume = " << name 
                  << "\n nameVolume = " << nameVolume << "\n"
                  << " nMedium = " << nMedium << "\n"
         //  << " rad length = " << rad << "\n"
         //  << " material = " << mat << "\n"
                  << " stepl = " << stepl << "\n"
                  << " photEnSpectrDE = " << photEnSpectrDE <<"\n"
                  << " edep = " << edep << "\n"
                  << " ===> calo theTrack info <=== " << "\n"
                  << " particleType = " << particleType << "\n"
                  << " primaryID = " << primaryID << "\n"
                  << " entot= " << theTrack->GetTotalEnergy() << "\n"
                  << " vert_eta= " << eta  << "\n"
                  << " vert_phi= " << phi << "\n"
                  << " vert_mom= " << vert_mom  << "\n"
                  << " ===> calo hit preStepPointinfo <=== "<<"\n"
                  << " local point = " << localPoint << "\n"
                  << " ==============================> final info"
                  << ":  <=== \n" 
                  << " meanNCherPhot = " << meanNCherPhot << "\n"
                  << " poissNCherPhot = " << poissNCherPhot <<"\n"
                  << " NCherPhot = " << NCherPhot;
 #endif 
     
       // Included by WC
       //         std::cout << "\n volume = "         << name 
       //              << "\n nameVolume = "     << nameVolume << "\n"
       //              << "\n postvolume = "     << postname 
       //              << "\n postnameVolume = " << postnameVolume << "\n"
       //              << "\n particleType = "   << particleType 
       //              << "\n primaryID = "      << primaryID << "\n";
 
     }
       }
 
 
 #ifdef debugLog
       LogDebug("ForwardSim") << "CastorSD:: " << nameVolume 
     //      << " Light Collection Efficiency " << weight
                  << " Weighted Energy Deposit " << edep/MeV 
                  << " MeV\n";
 #endif
       // Temporary member for testing purpose only...
       // unit_id = setDetUnitId(aStep);
       // if(NCherPhot != 0) std::cout << "\n  UnitID = " << unit_id << "  ;  NCherPhot = " << NCherPhot ;
        
       return NCherPhot;
     }
   } 
   return 0;
 }

void CastorSD::getFromLibrary ( G4Step * aStep )

private

Definition at line 610 of file CastorSD.cc.

References CaloSD::checkHit(), CaloSD::createNewHit(), CaloSD::currentHit, CaloSD::currentID, CaloSD::edepositEM, CaloSD::edepositHAD, CaloSD::emPDG, CaloSD::epPDG, CaloSD::gammaPDG, CastorShowerEvent::getDetID(), CastorShowerEvent::getHitPosition(), CastorShowerEvent::getNhit(), CastorShowerEvent::getNphotons(), CastorShowerEvent::getPrimE(), CastorShowerLibrary::getShowerHits(), CastorShowerEvent::getTime(), GetVolume(), i, LogDebug, convertSQLiteXML::ok, CaloSD::posGlobal, position, CaloSD::preStepPoint, CaloSD::previousID, CaloSD::resetForNewPrimary(), rotateUnitID(), CaloHitID::setID(), setTrackID(), showerLibrary, CaloSD::theTrack, cond::rpcobgas::time, funct::true, and CaloSD::updateHit().

Referenced by getEnergyDeposit().

                                             {
 
 //
 //   Method to get hits from the Shower Library
 //
 //   CastorShowerEvent hits returned by getShowerHits are used to  
 //   replace the full simulation of the shower from theTrack
 //    
 //   "updateHit" save the Hits to a CaloG4Hit container
 //
 
   preStepPoint  = aStep->GetPreStepPoint(); 
   theTrack      = aStep->GetTrack();   
   bool ok;
   
   // ****    Call method to retrieve hits from the ShowerLibrary   ****
   CastorShowerEvent hits = showerLibrary->getShowerHits(aStep, ok);
 
   double etrack    = preStepPoint->GetKineticEnergy();
   int    primaryID = setTrackID(aStep);
   // int    primaryID = theTrack->GetTrackID();
 
   // Reset entry point for new primary
   posGlobal = preStepPoint->GetPosition();
   resetForNewPrimary(posGlobal, etrack);
 
   // Check whether track is EM or HAD
   G4int particleCode = theTrack->GetDefinition()->GetPDGEncoding();
   bool isEM , isHAD ;
   if (particleCode==emPDG || particleCode==epPDG || particleCode==gammaPDG) {
     isEM = true ; isHAD = false;
   } else {
     isEM = false; isHAD = true ;
   }
 
 #ifdef debugLog
   //  edm::LogInfo("ForwardSim") << "\n CastorSD::getFromLibrary:  " 
   LogDebug("ForwardSim") << "\n CastorSD::getFromLibrary:  " 
              << hits.getNhit() << " hits for " << GetName() << " from " 
              << theTrack->GetDefinition()->GetParticleName() << " of "
              << preStepPoint->GetKineticEnergy()/GeV << " GeV and trackID " 
              << theTrack->GetTrackID()  ;
 #endif
 
   // Scale to correct energy
   double E_track = preStepPoint->GetTotalEnergy() ;
   double E_SLhit = hits.getPrimE() * GeV ;
   double scale = E_track/E_SLhit ;
   
 /*    double theTrackEnergy = theTrack->GetTotalEnergy() ; 
   
   if(fabs(theTrackEnergy-E_track)>10.) {
     edm::LogInfo("ForwardSim") << "\n            TrackID = " << theTrack->GetTrackID()
                                << "\n     theTrackEnergy = " << theTrackEnergy
                                << "\n preStepPointEnergy = " << E_track ;
     G4TrackVector tsec = *(aStep->GetSecondary());
     for (unsigned int kk=0; kk<tsec.size(); kk++) {
     edm::LogInfo("ForwardSim") << "CastorSD::getFromLibrary:"
                    << "\n tsec[" << kk << "]->GetTrackID() = " 
                    << tsec[kk]->GetTrackID() 
                    << " with energy " 
                    << tsec[kk]->GetTotalEnergy() ;
     }
   }
 */  
   //  Loop over hits retrieved from the library
   for (unsigned int i=0; i<hits.getNhit(); i++) {
     
     // Get nPhotoElectrons and set edepositEM / edepositHAD accordingly
     double nPhotoElectrons    = hits.getNphotons(i);
     // Apply scaling
       nPhotoElectrons *= scale ;
     if(isEM)  {
        // edepositEM  = nPhotoElectrons*GeV; 
        edepositEM  = nPhotoElectrons; 
        edepositHAD = 0.;                 
     } else if(isHAD) {
        edepositEM  = 0.;                  
        edepositHAD = nPhotoElectrons;
        // edepositHAD = nPhotoElectrons*GeV;
     }
     
     // Get hit position and time
     double                time = hits.getTime(i);
     math::XYZPoint    position = hits.getHitPosition(i);
     
     // Get hit detID
     unsigned int        unitID = hits.getDetID(i);
     
     // Make the detID "rotation" from one sector to another taking into account the 
     // sectors of the impinging particle (theTrack) and of the particle that produced 
     // the 'hits' retrieved from shower library   
     unsigned int rotatedUnitID = rotateUnitID(unitID , theTrack , hits);
     currentID.setID(rotatedUnitID, time, primaryID, 0);
     // currentID.setID(unitID, time, primaryID, 0);
    
     // check if it is in the same unit and timeslice as the previous one
     if (currentID == previousID) {
       updateHit(currentHit);
     } else {
       if (!checkHit()) currentHit = createNewHit();
     }
   }  //  End of loop over hits
 
   //Now kill the current track
   if (ok) {
     theTrack->SetTrackStatus(fStopAndKill);
 #ifdef debugLog
     LogDebug("ForwardSim") << "CastorSD::getFromLibrary:"
                << "\n \"theTrack\" with TrackID() = " 
                << theTrack->GetTrackID() 
                << " and with energy " 
                << theTrack->GetTotalEnergy()
                << " has been set to be killed" ;
 #endif
     G4TrackVector tv = *(aStep->GetSecondary());
     for (unsigned int kk=0; kk<tv.size(); kk++) {
       if (tv[kk]->GetVolume() == preStepPoint->GetPhysicalVolume()) {
     tv[kk]->SetTrackStatus(fStopAndKill);
 #ifdef debugLog
     LogDebug("ForwardSim") << "CastorSD::getFromLibrary:"
                    << "\n tv[" << kk << "]->GetTrackID() = " 
                    << tv[kk]->GetTrackID() 
                    << " with energy " 
                    << tv[kk]->GetTotalEnergy()
                    << " has been set to be killed" ;
 #endif
       }
     }
   }
 }

void CastorSD::initRun ( )

protectedvirtual

Reimplemented from CaloSD.

Definition at line 81 of file CastorSD.cc.

References CastorShowerLibrary::initParticleTable(), showerLibrary, and useShowerLibrary.

                       {
   if (useShowerLibrary) {
     // showerLibrary = new CastorShowerLibrary(name, cpv, p);
     G4ParticleTable * theParticleTable = G4ParticleTable::GetParticleTable();
     showerLibrary->initParticleTable(theParticleTable);
     edm::LogInfo("ForwardSim") << "CastorSD::initRun: Using Castor Shower Library \n";
   }
 }

uint32_t CastorSD::rotateUnitID	(	uint32_t	unitID,
		G4Track *	track,
		CastorShowerEvent	shower
	)

private

Definition at line 545 of file CastorSD.cc.

References printConversionInfo::aux, CastorShowerEvent::getPrimPhi(), LogDebug, M_PI, and reco::btau::trackPhi.

Referenced by getFromLibrary().

                                                                                          {
 // ==============================================================
 //
 //   o   Exploit Castor phi symmetry to return newUnitID for  
 //       shower hits based on track phi coordinate 
 //
 // ==============================================================
   
   // Get 'track' phi:
   float   trackPhi = track->GetPosition().phi(); 
   if(trackPhi<0) trackPhi += 2*M_PI ;
   // Get phi from primary that gave rise to SL 'shower':
   float  showerPhi = shower.getPrimPhi(); 
   if(showerPhi<0) showerPhi += 2*M_PI ;
   // Delta phi:
   
   //  Find the OctSector for which 'track' and 'shower' belong
   int  trackOctSector = (int) (  trackPhi / (M_PI/4) ) ;
   int showerOctSector = (int) ( showerPhi / (M_PI/4) ) ;
   
   uint32_t  newUnitID;
   uint32_t         sec = ( ( unitID>>4 ) & 0xF ) ;
   uint32_t  complement = ( unitID & 0xFFFFFF0F ) ;
   
   // edm::LogInfo("ForwardSim") << "\n CastorSD::rotateUnitID:  " 
   //                            << "\n      unitID = " << unitID 
   //                            << "\n         sec = " << sec 
   //                            << "\n  complement = " << complement ; 
   
   // Get 'track' z:
   float   trackZ = track->GetPosition().z();
   
   int aux ;
   int dSec = 2*(trackOctSector - showerOctSector) ;
   // if(trackZ<0)  // Good for revision 1.8 of CastorNumberingScheme
   if(trackZ>0)  // Good for revision 1.9 of CastorNumberingScheme
   {
     int sec1 = sec-dSec;
     //    sec -= dSec ;
     if(sec1<0) sec1  += 16;
     if(sec1>15) sec1 -= 16;
     sec = (uint32_t)(sec1);
   } else {
     if( dSec<0 ) sec += 16 ;
     sec += dSec ;
     aux  = (int) (sec/16) ;
     sec -= aux*16 ;
   }
   sec  = sec<<4 ;
   newUnitID = complement | sec ;
   
 #ifdef debugLog
   if(newUnitID != unitID) {
     LogDebug("ForwardSim") << "\n CastorSD::rotateUnitID:  " 
                << "\n     unitID = " << unitID 
                << "\n  newUnitID = " << newUnitID ; 
   }
 #endif
   
   return newUnitID ;
 
 }

uint32_t CastorSD::setDetUnitId ( G4Step * step )

virtual

Implements CaloSD.

Definition at line 504 of file CastorSD.cc.

References CastorNumberingScheme::getUnitID(), and numberingScheme.

                                              {
   return (numberingScheme == 0 ? 0 : numberingScheme->getUnitID(aStep));
 }

void CastorSD::setNumberingScheme ( CastorNumberingScheme * scheme )

Definition at line 510 of file CastorSD.cc.

References numberingScheme.

Referenced by CastorSD().

                                                                {
 
   if (scheme != 0) {
     edm::LogInfo("ForwardSim") << "CastorSD: updates numbering scheme for " 
                    << GetName();
     if (numberingScheme) delete numberingScheme;
     numberingScheme = scheme;
   }
 }

int CastorSD::setTrackID ( G4Step * aStep )

private

Definition at line 522 of file CastorSD.cc.

References TrackInformation::getIDonCaloSurface(), CaloSD::preStepPoint, CaloSD::previousID, CaloSD::resetForNewPrimary(), CaloSD::theTrack, and CaloHitID::trackID().

Referenced by getFromLibrary().

                                        {
 
   theTrack     = aStep->GetTrack();
 
   double etrack = preStepPoint->GetKineticEnergy();
   TrackInformation * trkInfo = (TrackInformation *)(theTrack->GetUserInformation());
   int      primaryID = trkInfo->getIDonCaloSurface();
   if (primaryID == 0) {
 #ifdef debugLog
     edm::LogWarning("ForwardSim") << "CastorSD: Problem with primaryID **** set by force "
                   << "to TkID **** " << theTrack->GetTrackID();
 #endif
     primaryID = theTrack->GetTrackID();
   }
 
   if (primaryID != previousID.trackID())
     resetForNewPrimary(preStepPoint->GetPosition(), etrack);
 
   return primaryID;
 }