diff --git a/SimCalorimetry/EcalSelectiveReadoutAlgos/src/EcalSelectiveReadoutSuppressor.cc b/SimCalorimetry/EcalSelectiveReadoutAlgos/src/EcalSelectiveReadoutSuppressor.cc index 8d086541accd4..4fe5de1211280 100644 --- a/SimCalorimetry/EcalSelectiveReadoutAlgos/src/EcalSelectiveReadoutSuppressor.cc +++ b/SimCalorimetry/EcalSelectiveReadoutAlgos/src/EcalSelectiveReadoutSuppressor.cc @@ -330,10 +330,10 @@ EcalSelectiveReadoutSuppressor::run(const edm::EventSetup& eventSetup, } if(ievt_ <= 10){ + int neb = (selectedBarrelDigis?selectedBarrelDigis->size():0); if(selectedEndcapDigis) LogDebug("EcalSelectiveReadout") // << __FILE__ << ":" << __LINE__ << ": " - << "Number of EB digis passing the SR: " - << selectedBarrelDigis->size() + << "Number of EB digis passing the SR: " << neb << " / " << barrelDigis.size() << "\n"; if(selectedEndcapDigis) LogDebug("EcalSelectiveReadout") // << __FILE__ << ":" << __LINE__ << ": " diff --git a/SimG4CMS/Calo/src/HFFibreFiducial.cc b/SimG4CMS/Calo/src/HFFibreFiducial.cc index cd08a7f7d66da..1d20fe99c0fc6 100644 --- a/SimG4CMS/Calo/src/HFFibreFiducial.cc +++ b/SimG4CMS/Calo/src/HFFibreFiducial.cc @@ -6,30 +6,18 @@ #include -//#define DebugLog -//#define mkodebug - int HFFibreFiducial::PMTNumber(const G4ThreeVector& pe_effect) { -#ifdef mkodebug - static double mX=0.; - static double mY=0.; - static double mZ=0.; -#endif double xv = pe_effect.x(); // X in global system double yv = pe_effect.y(); // Y in global system -#ifdef DebugLog double zv = pe_effect.z(); // Z in global system -#endif double phi = atan2(yv, xv); // In global system if (phi < 0.) phi+=CLHEP::pi; // Just for security double dph = CLHEP::pi/18; // 10 deg = a half sector width double sph = dph+dph; // 20 deg = a sector width int nphi = phi/dph; // 10 deg sector # -#ifdef DebugLog edm::LogInfo("HFShower") <<"HFFibreFiducial:***> P = " << pe_effect << ", phi = " << phi/CLHEP::deg; -#endif if (nphi > 35) nphi=35; // Just for security double xl=0.; // local sector coordinates (left/right) double yl=0.; // local sector coordinates (down/up) @@ -58,16 +46,12 @@ int HFFibreFiducial::PMTNumber(const G4ThreeVector& pe_effect) double sinr= sin(phir); yl= xv*cosr+yv*sinr; xl= yv*cosr-xv*sinr; -#ifdef DebugLog edm::LogInfo("HFShower") << "HFFibreFiducial: nr " << nr << " phi " << phir/CLHEP::deg; -#endif } if (yl < 0) yl =-yl; -#ifdef DebugLog edm::LogInfo("HFShower") << "HFFibreFiducial: Global Point " << pe_effect << " nphi " << nphi << " Local Sector Coordinates (" << xl << ", " << yl << "), widget # " << nwid; -#endif // Provides a PMT # for the (x,y) hit in the widget # nwid (M. Kosov, 11.2010) // Send comments/questions to Mikhail.Kossov@cern.ch // nwid = 1-18 for Forward HF, 19-36 for Backward HF (all equal now) @@ -76,24 +60,20 @@ int HFFibreFiducial::PMTNumber(const G4ThreeVector& pe_effect) static const int nWidM=36; if (nwid > nWidM || nwid <= 0) { -#ifdef DebugLog edm::LogInfo("HFShower") << "-Warning-HFFibreFiducial::PMTNumber: " << nwid << " == wrong widget number"; -#endif return 0; } static const double yMin= 13.1*CLHEP::cm; // start of the active area (Conv to mm?) static const double yMax=129.6*CLHEP::cm; // finish of the active area (Conv to mm?) if( yl < yMin || yl >= yMax ) { -#ifdef DebugLog edm::LogInfo("HFShower") << "-Warning-HFFibreFiducial::PMTNumber: Point " << "with y = " << yl << " outside acceptance [" << yMin << ":" << yMax << "], X = " << xv << ", Y = " << yv << ", x = " << xl << ", nW = " << nwid << ", phi = " << phi/CLHEP::deg << ", phir = " << phir/CLHEP::deg; -#endif return 0; // ===> out of the acceptance } bool left=true; // flag of the left part of the widget @@ -106,11 +86,9 @@ int HFFibreFiducial::PMTNumber(const G4ThreeVector& pe_effect) static const double tg10=.17632698070847; // phi-angular acceptance of the widget if (r > tg10) { -#ifdef DebugLog edm::LogInfo("HFShower") << "-Warning-HFFibreFiducial::PMTNumber: (x = " << xl << ", y = " << yl << ", tg = " << r << ") out of the widget acceptance tg(10) " << tg10; -#endif return 0; } @@ -1534,27 +1512,11 @@ int HFFibreFiducial::PMTNumber(const G4ThreeVector& pe_effect) int ny=static_cast((yl-yMin)/cellSize); // Layer number (starting from 0) if (ny < 0 || ny >= nLay) // Sould never happen as was checked beforehand { -#ifdef DebugLog edm::LogInfo("HFShower") << "-Warning-HFFibreFiducial::PMTNumber: " << "check limits y = " << yl << ", nL=" << nLay; -#endif return 0; } int nx=static_cast(fx); // Cell number (starting from 0) -#ifdef mkodebug - double phist=atan2(xl, yl)/CLHEP::deg; - if(sqr(mX-xv)+sqr(mY-yv)+sqr(mZ-zv) > 9.) - { - std::cout<<"HFFibreFiducial::PMTNumber:X="<nX="< out of the acceptance } int code=0; // a prototype @@ -1579,21 +1538,16 @@ int HFFibreFiducial::PMTNumber(const G4ThreeVector& pe_effect) int flag= code%10; int npmt= code/10; bool src= false; // by default: not a source-tube -#ifdef DebugLog edm::LogInfo("HFShower") << "HFFibreFiducial::nx/ny (" << nx << "," << ny << ") code/flag/npmt " << code << "/" << flag << "/" << npmt; -#endif if (!flag) return 0; // ===> no fiber in the cell else if (flag==1) npmt += 24; else if (flag==3 || flag==4) { - flag-=2; src=true; } -#ifdef DebugLog edm::LogInfo("HFShower") << "HFFibreFiducial::PMTNumber: src = " << src << ", npmt =" << npmt; -#endif - if (src) return -npmt; // return the negative number for the source + if (src) return -npmt; // return the negative number for the source return npmt; } // End of PMTNumber diff --git a/SimG4CMS/FP420/plugins/FP420Test.cc b/SimG4CMS/FP420/plugins/FP420Test.cc index fcfa235f605f6..d9ee5c5ca37ef 100644 --- a/SimG4CMS/FP420/plugins/FP420Test.cc +++ b/SimG4CMS/FP420/plugins/FP420Test.cc @@ -461,8 +461,8 @@ void FP420Test::update(const EndOfTrack * trk) { G4ThreeVector vert_pos = (*trk)()->GetVertexPosition(); // vertex ,where this track was created // float eta = 0.5 * log( (1.+vert_mom.z()) / (1.-vert_mom.z()) ); - float phi = atan2(vert_mom.y(),vert_mom.x()); - if (phi < 0.) phi += twopi; +// float phi = atan2(vert_mom.y(),vert_mom.x()); +// if (phi < 0.) phi += twopi; // float phigrad = phi*180./pi; // float XV = vert_pos.x(); // mm @@ -1082,7 +1082,6 @@ void FP420Test::update(const EndOfEvent * evt) { // Silicon Hit collection start //0) if particle goes into flat beam pipe below detector: int varia ; // = 0 -all; =1 - MI; =2 - noMI - varia = 0; // Select MI or noMI over all 3 stations // 1)MI: // if particle goes through window into detector: @@ -1154,8 +1153,8 @@ void FP420Test::update(const EndOfEvent * evt) { // double th_hit = hitPoint.theta(); // double eta_hit = -log(tan(th_hit/2)); - double phi_hit = hitPoint.phi(); - if (phi_hit < 0.) phi_hit += twopi; +// double phi_hit = hitPoint.phi(); +// if (phi_hit < 0.) phi_hit += twopi; // double phigrad_hit = phi_hit*180./pi; //UserNtuples->fillg60(eta_hit,losenergy); //UserNtuples->fillg61(eta_hit,1.); diff --git a/SimG4CMS/Forward/src/BscTest.cc b/SimG4CMS/Forward/src/BscTest.cc index cd2c116e7ee72..99866bbf3092f 100644 --- a/SimG4CMS/Forward/src/BscTest.cc +++ b/SimG4CMS/Forward/src/BscTest.cc @@ -335,12 +335,13 @@ void BscTest::update(const EndOfTrack * trk) { G4ThreeVector vert_pos = (*trk)()->GetVertexPosition(); // vertex ,where this track was created // float eta = 0.5 * log( (1.+vert_mom.z()) / (1.-vert_mom.z()) ); + /* float phi = atan2(vert_mom.y(),vert_mom.x()); if (phi < 0.) phi += twopi; if(tracklength < z4) { } - + */ // last step information const G4Step* aStep = (*trk)()->GetStep(); G4StepPoint* preStepPoint = aStep->GetPreStepPoint(); @@ -715,14 +716,13 @@ void BscTest::update(const EndOfEvent * evt) { std::cout << "BscTest: theCAFI->entries = " << theCAFI->entries() << std::endl; } int varia ; // = 0 -all; =1 - MI; =2 - noMI - varia = 0; + //varia = 0; if( lastpo.z()< z4) { varia = 1; } else{ varia = 2; } // no MI end: - for (int j=0; jentries(); j++) { BscG4Hit* aHit = (*theCAFI)[j]; CLHEP::Hep3Vector hitPoint = aHit->getEntry(); @@ -746,8 +746,8 @@ void BscTest::update(const EndOfEvent * evt) { int trackIDhit = aHit->getTrackID(); unsigned int unitID = aHit->getUnitID(); double losenergy = aHit->getEnergyLoss(); - double phi_hit = hitPoint.phi(); - if (phi_hit < 0.) phi_hit += twopi; + //double phi_hit = hitPoint.phi(); + //if (phi_hit < 0.) phi_hit += twopi; double zz = hitPoint.z(); @@ -937,11 +937,7 @@ void BscTest::update(const EndOfEvent * evt) { } else{ //UserNtuples->fillp212(vy,float(0.),1.); } - } // MI or no MI or all - end - - - } // primary end if (verbosity > 0) { diff --git a/SimG4CMS/Forward/src/CastorSD.cc b/SimG4CMS/Forward/src/CastorSD.cc index 3d5b7dfc6a974..3983bd8cedad8 100644 --- a/SimG4CMS/Forward/src/CastorSD.cc +++ b/SimG4CMS/Forward/src/CastorSD.cc @@ -410,7 +410,6 @@ double CastorSD::getEnergyDeposit(G4Step * aStep) { // if((thcher + DelFibPart ) > thFullReflRad && // thcher < (DelFibPart+thFullReflRad) ) // { - d_qz = 0.; #ifdef debugLog variant=3.; #endif @@ -422,10 +421,7 @@ double CastorSD::getEnergyDeposit(G4Step * aStep) { if(tan_arcos != 0.) arg_arcos =(r*r-a*a-d*d)/tan_arcos; arg_arcos = fabs(arg_arcos); double th_arcos = acos(std::min(std::max(arg_arcos,double(-1.)),double(1.))); - d_qz = th_arcos/pi/2.; - d_qz = fabs(d_qz); - - + d_qz = fabs(th_arcos/pi/2.); // } // else diff --git a/SimG4CMS/Forward/src/CastorShowerLibrary.cc b/SimG4CMS/Forward/src/CastorShowerLibrary.cc index b45cde50ef9f3..b2d41a62b68c7 100644 --- a/SimG4CMS/Forward/src/CastorShowerLibrary.cc +++ b/SimG4CMS/Forward/src/CastorShowerLibrary.cc @@ -15,6 +15,7 @@ #include "G4Track.hh" #include "G4ParticleTable.hh" #include "Randomize.hh" +#include "G4PhysicalConstants.hh" #include "CLHEP/Units/GlobalSystemOfUnits.h" //#define DebugLog @@ -233,14 +234,14 @@ CastorShowerEvent CastorShowerLibrary::getShowerHits(G4Step * aStep, bool & ok) ok = true; double pin = preStepPoint->GetTotalEnergy(); - double etain = momDir.getEta(); - double phiin = momDir.getPhi(); + // double etain = momDir.getEta(); + //double phiin = momDir.getPhi(); double zint = hitPoint.z(); double R=sqrt(hitPoint.x()*hitPoint.x() + hitPoint.y()*hitPoint.y()); double theta = atan2(R,std::abs(zint)); - phiin = atan2(hitPoint.y(),hitPoint.x()); - etain = -1*(std::log(std::tan((M_PI-theta)/2))); + double phiin = atan2(hitPoint.y(),hitPoint.x()); + double etain = -std::log(std::tan((pi-theta)*0.5)); // Replace "interpolation/extrapolation" by new method "select" that just randomly // selects a record from the appropriate energy bin and fills its content to diff --git a/SimG4CMS/Forward/src/DoCastorAnalysis.cc b/SimG4CMS/Forward/src/DoCastorAnalysis.cc index f14522e8797ee..aa54675992eb0 100644 --- a/SimG4CMS/Forward/src/DoCastorAnalysis.cc +++ b/SimG4CMS/Forward/src/DoCastorAnalysis.cc @@ -214,6 +214,7 @@ void DoCastorAnalysis::update(const EndOfEvent * evt) { CastorTree->Fill(); } // nentries > 0 + delete theCastorNumScheme; } void DoCastorAnalysis::update(const EndOfRun * run) {;} diff --git a/SimG4CMS/Forward/src/ZdcTestAnalysis.cc b/SimG4CMS/Forward/src/ZdcTestAnalysis.cc index 55e8f1a9a9d89..7abd99fbc8bad 100644 --- a/SimG4CMS/Forward/src/ZdcTestAnalysis.cc +++ b/SimG4CMS/Forward/src/ZdcTestAnalysis.cc @@ -342,14 +342,16 @@ void ZdcTestAnalysis::update(const EndOfEvent * evt) { for (int i = 0 ; iGetPrimaryVertex(i); - if (avertex == 0) + if (avertex == 0) { std::cout << "ZdcTest End Of Event ERR: pointer to vertex = 0" << std::endl; - std::cout << "Vertex number :" <GetNumberOfParticle(); - if (npart ==0) - std::cout << "ZdcTest End Of Event ERR: no primary!" << std::endl; - if (thePrim==0) thePrim=avertex->GetPrimary(trackID); + } else { + std::cout << "Vertex number :" <GetNumberOfParticle(); + if (npart ==0) + std::cout << "ZdcTest End Of Event ERR: no primary!" << std::endl; + if (thePrim==0) thePrim=avertex->GetPrimary(trackID); + } } double px=0.,py=0.,pz=0.; diff --git a/SimG4CMS/HcalTestBeam/plugins/HcalTB02Analysis.cc b/SimG4CMS/HcalTestBeam/plugins/HcalTB02Analysis.cc index e8e10b11596fb..0b96c63c247e4 100644 --- a/SimG4CMS/HcalTestBeam/plugins/HcalTB02Analysis.cc +++ b/SimG4CMS/HcalTestBeam/plugins/HcalTB02Analysis.cc @@ -376,6 +376,8 @@ void HcalTB02Analysis::update(const EndOfEvent * evt) { std::cout << " Event " << iEvt << std::endl; else if ((iEvt < 10000) && (iEvt%1000 == 0)) std::cout << " Event " << iEvt << std::endl; + + delete org; } void HcalTB02Analysis::fillEvent(HcalTB02HistoClass& product) { diff --git a/SimG4Core/Application/src/RunManagerMTWorker.cc b/SimG4Core/Application/src/RunManagerMTWorker.cc index 2282100af5069..e1517235a1770 100644 --- a/SimG4Core/Application/src/RunManagerMTWorker.cc +++ b/SimG4Core/Application/src/RunManagerMTWorker.cc @@ -362,13 +362,13 @@ void RunManagerMTWorker::initializeRun() { } void RunManagerMTWorker::terminateRun() { - if(m_tls->userRunAction) { + if(m_tls && m_tls->userRunAction) { m_tls->userRunAction->EndOfRunAction(m_tls->currentRun.get()); m_tls->userRunAction.reset(); } G4RunManagerKernel *kernel = G4WorkerRunManagerKernel::GetRunManagerKernel(); - if(!kernel && !m_tls->runTerminated) { + if(kernel && m_tls && !m_tls->runTerminated) { m_tls->currentEvent.reset(); m_simEvent.reset(); kernel->RunTermination(); diff --git a/SimG4Core/CustomPhysics/src/CustomParticleFactory.cc b/SimG4Core/CustomPhysics/src/CustomParticleFactory.cc index dd1b02fe9dbc7..326cf9b7c7804 100644 --- a/SimG4Core/CustomPhysics/src/CustomParticleFactory.cc +++ b/SimG4Core/CustomPhysics/src/CustomParticleFactory.cc @@ -85,8 +85,8 @@ void CustomParticleFactory::addCustomParticle(int pdgCode, double mass, const st /////////////////////// Check!!!!!!!!!!!!! G4String pType="custom"; G4String pSubType=""; - G4double spectatormass; - G4ParticleDefinition* spectator; + G4double spectatormass = 0.0; + G4ParticleDefinition* spectator = nullptr; ////////////////////// if(CustomPDGParser::s_isRHadron(pdgCode)) pType = "rhadron"; if(CustomPDGParser::s_isSLepton(pdgCode)) pType = "sLepton"; @@ -127,7 +127,7 @@ void CustomParticleFactory::addCustomParticle(int pdgCode, double mass, const st bool stable = true; double lifetime = -1; - G4DecayTable *decaytable = NULL; + G4DecayTable *decaytable = nullptr; G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable(); CustomParticle *particle = new CustomParticle(name, massGeV, width, charge, spin, @@ -170,11 +170,11 @@ void CustomParticleFactory::addCustomParticle(int pdgCode, double mass, const st if (CustomPDGParser::s_issbottomHadron(pdgCode)) { spectator = theParticleTable->FindParticle(1000005*sign); } else { - spectator = 0; + spectator = nullptr; edm::LogError("SimG4CoreCustomPhysics")<< "CustomParticleFactory: Cannot find spectator parton"; } } - spectatormass = spectator->GetPDGMass(); + if(spectator) spectatormass = spectator->GetPDGMass(); G4double cloudmass = mass-spectatormass/GeV; CustomParticle *tmpParticle = new CustomParticle( cloudname, cloudmass * GeV , 0.0*MeV, 0 , diff --git a/SimG4Core/CustomPhysics/src/DummyChargeFlipProcess.cc b/SimG4Core/CustomPhysics/src/DummyChargeFlipProcess.cc index 3d60f8d1eb0f1..81e04ddfedb80 100644 --- a/SimG4Core/CustomPhysics/src/DummyChargeFlipProcess.cc +++ b/SimG4Core/CustomPhysics/src/DummyChargeFlipProcess.cc @@ -62,7 +62,7 @@ G4VParticleChange *DummyChargeFlipProcess::PostStepDoIt( G4ParticleChange * pc = new G4ParticleChange(); pc->Initialize(aTrack); const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle(); - G4ParticleDefinition* aParticleDef = aParticle->GetDefinition(); + //G4ParticleDefinition* aParticleDef = aParticle->GetDefinition(); G4double ParentEnergy = aParticle->GetTotalEnergy(); G4ThreeVector ParentDirection(aParticle->GetMomentumDirection()); @@ -74,31 +74,31 @@ G4VParticleChange *DummyChargeFlipProcess::PostStepDoIt( pc->SetNumberOfSecondaries(numberOfSecondaries); const G4TouchableHandle thand = aTrack.GetTouchableHandle(); - // get current position of the track - aTrack.GetPosition(); - // create a new track object - G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable(); - float randomParticle = G4UniformRand(); - G4ParticleDefinition * newType = aParticleDef; - if(randomParticle < 0.333) - newType=particleTable->FindParticle(1009213); - else if(randomParticle > 0.667) - newType=particleTable->FindParticle(-1009213); - else - newType=particleTable->FindParticle(1009113); + // get current position of the track + aTrack.GetPosition(); + // create a new track object + G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable(); + float randomParticle = G4UniformRand(); + G4ParticleDefinition * newType; + if(randomParticle < 0.333) + newType=particleTable->FindParticle(1009213); + else if(randomParticle > 0.667) + newType=particleTable->FindParticle(-1009213); + else + newType=particleTable->FindParticle(1009113); - std::cout << "RHADRON: New charge = " << newType->GetPDGCharge() << std::endl; + G4cout << "RHADRON: New charge = " << newType->GetPDGCharge() << G4endl; - G4DynamicParticle * newP = new G4DynamicParticle(newType,ParentDirection,ParentEnergy); - G4Track* secondary = new G4Track( newP , - finalGlobalTime , - aTrack.GetPosition() - ); - // switch on good for tracking flag - secondary->SetGoodForTrackingFlag(); - secondary->SetTouchableHandle(thand); - // add the secondary track in the List - pc->AddSecondary(secondary); + G4DynamicParticle * newP = new G4DynamicParticle(newType,ParentDirection,ParentEnergy); + G4Track* secondary = new G4Track( newP , + finalGlobalTime , + aTrack.GetPosition() + ); + // switch on good for tracking flag + secondary->SetGoodForTrackingFlag(); + secondary->SetTouchableHandle(thand); + // add the secondary track in the List + pc->AddSecondary(secondary); // Kill the parent particle pc->ProposeTrackStatus( fStopAndKill ) ; @@ -107,9 +107,6 @@ G4VParticleChange *DummyChargeFlipProcess::PostStepDoIt( // Clear NumberOfInteractionLengthLeft ClearNumberOfInteractionLengthLeft(); - - - return pc; - + return pc; } diff --git a/SimG4Core/CustomPhysics/src/FullModelHadronicProcess.cc b/SimG4Core/CustomPhysics/src/FullModelHadronicProcess.cc index 289fcc4b5063b..5e1400353596d 100644 --- a/SimG4Core/CustomPhysics/src/FullModelHadronicProcess.cc +++ b/SimG4Core/CustomPhysics/src/FullModelHadronicProcess.cc @@ -370,7 +370,9 @@ G4VParticleChange* FullModelHadronicProcess::PostStepDoIt(const G4Track& aTrack, G4LorentzVector p_g_cms = gluinoMomentum; //gluino in CMS BEFORE collision p_g_cms.boost(trafo_full_cms); - G4LorentzVector p4_new( cloud_p4_new.v() + gluinoMomentum.v(), outgoingRhadron->GetPDGMass() ); + double e = cloud_p4_new.e() + gluinoMomentum.e(); + if(outgoingRhadron) e += outgoingRhadron->GetPDGMass(); + G4LorentzVector p4_new( cloud_p4_new.v() + gluinoMomentum.v(), e ); // G4cout<<"P4-diff: "<<(p4_new-cloud_p4_new-gluinoMomentum)/GeV<<", magnitude: " // <<(p4_new-cloud_p4_new-gluinoMomentum).m()/MeV<<" MeV" < 19 ) // commented out to duplicate ?bug? in GENXPT - x = 0.999; + //if( extraNucleonCount > 19 ) // commented out to duplicate ?bug? in GENXPT + // x = 0.999; G4double xxx = 0.95+0.05*extraNucleonCount/20.0; if( (backwardKinetic+kineticEnergy) < xxx*backwardEnergy ) { @@ -697,11 +697,13 @@ using namespace CLHEP; if( --vecLen == 0 )return false; // all the secondaries have been eliminated pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5]; pseudoParticle[6].SetMomentum( 0.0 ); // set z-momentum + /* phi = pseudoParticle[6].Angle( pseudoParticle[8] ); if( pseudoParticle[6].GetMomentum().y()/MeV < 0.0 )phi = twopi - phi; phi += pi + normal() * pi / 12.0; if( phi > twopi )phi -= twopi; if( phi < 0.0 )phi = twopi - phi; + */ } } // closes main for loop @@ -792,7 +794,7 @@ using namespace CLHEP; et = pseudoParticle[1].GetTotalEnergy()/GeV; dndl[0] = 0.0; outerCounter = 0; - eliminateThisParticle = true; // should never eliminate the target particle + //eliminateThisParticle = true; // should never eliminate the target particle resetEnergies = true; while( ++outerCounter < 3 ) // start of outer iteration loop { @@ -820,7 +822,7 @@ using namespace CLHEP; targetParticle.SetTotalEnergy( totalEnergy*GeV ); if( targetParticle.GetSide() < 0 ) { - if( extraNucleonCount > 19 )x=0.999; + //if( extraNucleonCount > 19 )x=0.999; G4double xxx = 0.95+0.05*extraNucleonCount/20.0; if( (backwardKinetic+totalEnergy-vecMass) < xxx*backwardEnergy ) { @@ -828,13 +830,15 @@ using namespace CLHEP; backwardKinetic += totalEnergy - vecMass; pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5]; pseudoParticle[6].SetMomentum( 0.0 ); // set z-momentum + /* phi = pseudoParticle[6].Angle( pseudoParticle[8] ); if( pseudoParticle[6].GetMomentum().y()/MeV < 0.0 )phi = twopi - phi; phi += pi + normal() * pi / 12.0; if( phi > twopi )phi -= twopi; if( phi < 0.0 )phi = twopi - phi; + */ outerCounter = 2; // leave outer loop - eliminateThisParticle = false; // don't eliminate this particle + //eliminateThisParticle = false; // don't eliminate this particle resetEnergies = false; break; // leave inner loop } @@ -874,7 +878,7 @@ using namespace CLHEP; pseudoParticle[4] = pseudoParticle[4] + targetParticle; outerCounter = 2; // leave outer loop - eliminateThisParticle = false; // don't eliminate this particle + //eliminateThisParticle = false; // don't eliminate this particle resetEnergies = false; break; // leave inner loop } @@ -1034,7 +1038,7 @@ using namespace CLHEP; // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen); if( tempLen >= 2 ) { - wgt = GenerateNBodyEvent( + GenerateNBodyEvent( pseudoParticle[6].GetMass(), constantCrossSection, tempV, tempLen ); // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen); if( targetParticle.GetSide() == -3 ) @@ -1516,7 +1520,7 @@ using namespace CLHEP; delete vec[vecLen-2]; vecLen -= 2; currentMass = currentParticle.GetMass()/GeV; - targetMass = targetParticle.GetMass()/GeV; + //targetMass = targetParticle.GetMass()/GeV; incidentHasChanged = true; targetHasChanged = true; currentParticle.SetKineticEnergy( ek ); @@ -1911,7 +1915,7 @@ using namespace CLHEP; } if( tempLen >= 2 ) { - wgt = GenerateNBodyEvent( pseudoParticle[3].GetMass()/MeV, + GenerateNBodyEvent( pseudoParticle[3].GetMass()/MeV, constantCrossSection, tempV, tempLen ); if( currentParticle.GetSide() == 1 ) currentParticle.Lorentz( currentParticle, pseudoParticle[5] ); @@ -1951,7 +1955,7 @@ using namespace CLHEP; } if( tempLen >= 2 ) { - wgt = GenerateNBodyEvent( pseudoParticle[4].GetMass()/MeV, + GenerateNBodyEvent( pseudoParticle[4].GetMass()/MeV, constantCrossSection, tempV, tempLen ); if( currentParticle.GetSide() == -1 ) currentParticle.Lorentz( currentParticle, pseudoParticle[6] ); @@ -2067,8 +2071,8 @@ using namespace CLHEP; { rthnve = pi*G4UniformRand(); phinve = twopi*G4UniformRand(); - targetParticle.SetMomentum( pp*std::sin(rthnve)*std::cos(rthnve)*MeV, - pp*std::sin(rthnve)*std::sin(rthnve)*MeV, + targetParticle.SetMomentum( pp*std::sin(rthnve)*std::cos(phinve)*MeV, + pp*std::sin(rthnve)*std::sin(phinve)*MeV, pp*std::cos(rthnve)*MeV ); } else @@ -2087,8 +2091,8 @@ using namespace CLHEP; { rthnve = pi*G4UniformRand(); phinve = twopi*G4UniformRand(); - currentParticle.SetMomentum( pp*std::sin(rthnve)*std::cos(rthnve)*MeV, - pp*std::sin(rthnve)*std::sin(rthnve)*MeV, + currentParticle.SetMomentum( pp*std::sin(rthnve)*std::cos(phinve)*MeV, + pp*std::sin(rthnve)*std::sin(phinve)*MeV, pp*std::cos(rthnve)*MeV ); } else @@ -2139,7 +2143,7 @@ using namespace CLHEP; if( tempLen >= 2 ) { // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen); - wgt = GenerateNBodyEvent( + GenerateNBodyEvent( pseudoParticle[4].GetTotalEnergy()/MeV+pseudoParticle[5].GetTotalEnergy()/MeV, constantCrossSection, tempV, tempLen ); theoreticalKinetic = 0.0; @@ -2320,9 +2324,7 @@ using namespace CLHEP; // const G4double mOriginal = modifiedOriginal.GetMass()/GeV; const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/GeV; G4double currentMass = currentParticle.GetMass()/GeV; - G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV; - - targetMass = targetParticle.GetMass()/GeV; + G4double targetMass = targetParticle.GetMass()/GeV; const G4double atomicWeight = targetNucleus.GetN_asInt(); // G4cout<<"Atomic weight is found to be: "< total energy (" // << totalEnergy << "MeV)" << G4endl; - totalE = totalMass; + //totalE = totalMass; //delete [] mass; //delete [] energy; //for( i=0; i<3; ++i )delete [] pcm[i]; diff --git a/SimG4Core/PhysicsLists/src/UrbanMscModel93.cc b/SimG4Core/PhysicsLists/src/UrbanMscModel93.cc index e66d55b0adcd2..cbbac7ac251f3 100644 --- a/SimG4Core/PhysicsLists/src/UrbanMscModel93.cc +++ b/SimG4Core/PhysicsLists/src/UrbanMscModel93.cc @@ -668,7 +668,7 @@ G4double UrbanMscModel93::ComputeGeomPathLength(G4double) return zPathLength; } - G4double zmean = tPathLength; + G4double zmean; if (tPathLength < currentRange*dtrl) { if(tau < taulim) zmean = tPathLength*(1.-0.5*tau) ; else zmean = lambda0*(1.-G4Exp(-tau)); diff --git a/SimMuon/RPCDigitizer/src/RPCSimAsymmetricCls.cc b/SimMuon/RPCDigitizer/src/RPCSimAsymmetricCls.cc index 9f5f35c7afb4d..fe47d5d846450 100644 --- a/SimMuon/RPCDigitizer/src/RPCSimAsymmetricCls.cc +++ b/SimMuon/RPCDigitizer/src/RPCSimAsymmetricCls.cc @@ -190,9 +190,10 @@ RPCSimAsymmetricCls::simulate(const RPCRoll* roll, int lstrip=centralStrip; // Compute the cluster size - double w = CLHEP::RandFlat::shoot(engine); - LogDebug ("RPCSimAsymmetricCls")<<"[RPCSimAsymmetricCls::simulate] Fired RandFlat :: "<getClSize(rpcId.rawId(),posX, engine); // This is for cluster size chamber by chamber LogDebug ("RPCSimAsymmetricCls")<<"Clustersize = "<getClSize(posX, engine); std::vector cls; diff --git a/SimMuon/RPCDigitizer/src/RPCSimAverageNoise.cc b/SimMuon/RPCDigitizer/src/RPCSimAverageNoise.cc index 974293491b27b..3d11b358ce41b 100644 --- a/SimMuon/RPCDigitizer/src/RPCSimAverageNoise.cc +++ b/SimMuon/RPCDigitizer/src/RPCSimAverageNoise.cc @@ -147,8 +147,8 @@ RPCSimAverageNoise::simulate(const RPCRoll* roll, int fstrip=centralStrip; int lstrip=centralStrip; // Compute the cluster size - double w = CLHEP::RandFlat::shoot(engine); - if (w < 1.e-10) w=1.e-10; + //double w = CLHEP::RandFlat::shoot(engine); + //if (w < 1.e-10) w=1.e-10; int clsize = this->getClSize(posX, engine); std::vector cls; diff --git a/SimMuon/RPCDigitizer/src/RPCSimAverageNoiseEff.cc b/SimMuon/RPCDigitizer/src/RPCSimAverageNoiseEff.cc index a14354fb40038..d2a7649410daa 100644 --- a/SimMuon/RPCDigitizer/src/RPCSimAverageNoiseEff.cc +++ b/SimMuon/RPCDigitizer/src/RPCSimAverageNoiseEff.cc @@ -162,8 +162,8 @@ RPCSimAverageNoiseEff::simulate(const RPCRoll* roll, int lstrip=centralStrip; // Compute the cluster size - double w = CLHEP::RandFlat::shoot(engine); - if (w < 1.e-10) w=1.e-10; + //double w = CLHEP::RandFlat::shoot(engine); + //if (w < 1.e-10) w=1.e-10; int clsize = this->getClSize(posX, engine); std::vector cls; diff --git a/SimMuon/RPCDigitizer/src/RPCSimAverageNoiseEffCls.cc b/SimMuon/RPCDigitizer/src/RPCSimAverageNoiseEffCls.cc index 6d1c38b0f1f82..1284d083f03c8 100644 --- a/SimMuon/RPCDigitizer/src/RPCSimAverageNoiseEffCls.cc +++ b/SimMuon/RPCDigitizer/src/RPCSimAverageNoiseEffCls.cc @@ -212,8 +212,8 @@ RPCSimAverageNoiseEffCls::simulate(const RPCRoll* roll, int lstrip=centralStrip; // Compute the cluster size - double w = CLHEP::RandFlat::shoot(engine); - if (w < 1.e-10) w=1.e-10; + // double w = CLHEP::RandFlat::shoot(engine); + //if (w < 1.e-10) w=1.e-10; // int clsize = this->getClSize(posX, engine); // This is for one and the same cls for all the chambers int clsize = this->getClSize(rpcId.rawId(),posX, engine); // This is for cluster size chamber by chamber std::vector cls; diff --git a/SimMuon/RPCDigitizer/src/RPCSimSetUp.cc b/SimMuon/RPCDigitizer/src/RPCSimSetUp.cc index fd6ffcfd21d59..4083cb23e2af2 100644 --- a/SimMuon/RPCDigitizer/src/RPCSimSetUp.cc +++ b/SimMuon/RPCDigitizer/src/RPCSimSetUp.cc @@ -223,7 +223,7 @@ void RPCSimSetUp::setRPCSetUp(const std::vector& vnoi } sslognoiseitem <<"Start Position :: current_detId = "<nstrips()):(0))<<" strips"<nstrips():0)<<" strips"<::const_iterator it = vnoise.begin(); it != vnoise.end(); ++it) { diff --git a/SimTracker/Common/src/SiG4UniversalFluctuation.cc b/SimTracker/Common/src/SiG4UniversalFluctuation.cc index 765748414d092..6fd3ebcac7f07 100644 --- a/SimTracker/Common/src/SiG4UniversalFluctuation.cc +++ b/SimTracker/Common/src/SiG4UniversalFluctuation.cc @@ -173,23 +173,8 @@ double SiG4UniversalFluctuation::SampleFluctuations(const double momentum, } } - // Glandz regime : initialisation - // -// if (material != lastMaterial) { -// f1Fluct = material->GetIonisation()->GetF1fluct(); -// f2Fluct = material->GetIonisation()->GetF2fluct(); -// e1Fluct = material->GetIonisation()->GetEnergy1fluct(); -// e2Fluct = material->GetIonisation()->GetEnergy2fluct(); -// e1LogFluct = material->GetIonisation()->GetLogEnergy1fluct(); -// e2LogFluct = material->GetIonisation()->GetLogEnergy2fluct(); -// rateFluct = material->GetIonisation()->GetRateionexcfluct(); -// ipotFluct = material->GetIonisation()->GetMeanExcitationEnergy(); -// ipotLogFluct = material->GetIonisation()->GetLogMeanExcEnergy(); -// lastMaterial = material; -// } - - double a1 = 0. , a2 = 0., a3 = 0. ; - double p1,p2,p3; + double a1 = 0., a2 = 0., a3 = 0.; + double p3; double rate = rateFluct ; double w1 = tmax/ipotFluct; @@ -222,7 +207,7 @@ double SiG4UniversalFluctuation::SampleFluctuations(const double momentum, // if (suma > sumalim) { - p1 = 0., p2 = 0 ; + double p1, p2; if((a1+a2) > 0.) { // excitation type 1 diff --git a/SimTracker/SiStripDigitizer/plugins/SiTrivialInduceChargeOnStrips.cc b/SimTracker/SiStripDigitizer/plugins/SiTrivialInduceChargeOnStrips.cc index 1b23e5112462f..5af217722e4f9 100644 --- a/SimTracker/SiStripDigitizer/plugins/SiTrivialInduceChargeOnStrips.cc +++ b/SimTracker/SiStripDigitizer/plugins/SiTrivialInduceChargeOnStrips.cc @@ -191,8 +191,13 @@ induceVector(const SiChargeCollectionDrifter::collection_type& collection_points for (int i=0; i!=N;++i) { auto delta = 1.f/(std::sqrt(2.f)*chargeSpread[i]); auto pos = delta*(float(fromStrip[i])-chargePosition[i]); - for (int j=0;j<=nStrip[i]; ++j) /// include last strip - value[kk++] = pos+float(j)*delta; + + // VI: before value[0] was not defined and value[tot] was filled + // to fix this the loop below was changed + for (int j=0;j<=nStrip[i]; ++j) { /// include last strip + value[kk] = pos+float(j)*delta; + ++kk; + } } assert(kk==tot);