Reference Manual

Go to the documentation of this file.
00001 
00002 // File: CastorSD.cc
00003 // Date: 02.04
00004 // UpDate: 07.04 - C3TF & C4TF semi-trapezoid added
00005 // Description: Sensitive Detector class for Castor
00007 
00008 //  Added by WC 
00009 #include "SimG4Core/Notification/interface/TrackInformation.h"
00010 
00011 #include "SimG4CMS/Forward/interface/CastorSD.h"
00012 //#include "SimDataFormats/CaloHit/interface/CastorShowerEvent.h"
00013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00014 
00015 #include "G4SDManager.hh"
00016 #include "G4Step.hh"
00017 #include "G4Track.hh"
00018 #include "G4VProcess.hh"
00019 
00020 #include "G4ios.hh"
00021 #include "G4Cerenkov.hh"
00022 #include "G4LogicalVolumeStore.hh"
00023 
00024 #include "CLHEP/Units/GlobalSystemOfUnits.h"
00025 #include "Randomize.hh"
00026 #include "G4Poisson.hh"
00027 
00028 //#define debugLog
00029 
00030 CastorSD::CastorSD(G4String name, const DDCompactView & cpv,
00031                    SensitiveDetectorCatalog & clg, 
00032                    edm::ParameterSet const & p, 
00033                    const SimTrackManager* manager) : 
00034   CaloSD(name, cpv, clg, p, manager), numberingScheme(0), lvC3EF(0),
00035   lvC3HF(0), lvC4EF(0), lvC4HF(0), lvCAST(0) {
00036   
00037   edm::ParameterSet m_CastorSD = p.getParameter<edm::ParameterSet>("CastorSD");
00038   useShowerLibrary  = m_CastorSD.getParameter<bool>("useShowerLibrary");
00039   energyThresholdSL = m_CastorSD.getParameter<double>("minEnergyInGeVforUsingSLibrary");
00040   energyThresholdSL = energyThresholdSL*GeV;   //  Convert GeV => MeV 
00041           
00042   non_compensation_factor = m_CastorSD.getParameter<double>("nonCompensationFactor");
00043   
00044   if (useShowerLibrary) showerLibrary = new CastorShowerLibrary(name, p);
00045   
00046   setNumberingScheme(new CastorNumberingScheme());
00047   
00048   edm::LogInfo("ForwardSim") 
00049     << "***************************************************\n"
00050     << "*                                                 *\n" 
00051     << "* Constructing a CastorSD  with name " << GetName() << "\n"
00052     << "*                                                 *\n"
00053     << "***************************************************";
00054     
00055   const G4LogicalVolumeStore * lvs = G4LogicalVolumeStore::GetInstance();
00056   std::vector<G4LogicalVolume*>::const_iterator lvcite;
00057   for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
00058     if (strcmp(((*lvcite)->GetName()).c_str(),"C3EF") == 0) lvC3EF = (*lvcite);
00059     if (strcmp(((*lvcite)->GetName()).c_str(),"C3HF") == 0) lvC3HF = (*lvcite);
00060     if (strcmp(((*lvcite)->GetName()).c_str(),"C4EF") == 0) lvC4EF = (*lvcite);
00061     if (strcmp(((*lvcite)->GetName()).c_str(),"C4HF") == 0) lvC4HF = (*lvcite);
00062     if (strcmp(((*lvcite)->GetName()).c_str(),"CAST") == 0) lvCAST = (*lvcite);
00063     if (lvC3EF != 0 && lvC3HF != 0 && lvC4EF != 0 && lvC4HF != 0 && lvCAST != 0) break;
00064   }
00065   edm::LogInfo("ForwardSim") << "CastorSD:: LogicalVolume pointers\n"
00066                              << lvC3EF << " for C3EF; " << lvC3HF 
00067                              << " for C3HF; " << lvC4EF << " for C4EF; " 
00068                              << lvC4HF << " for C4HF; " 
00069                              << lvCAST << " for CAST. " << std::endl;
00070 
00071   //  if(useShowerLibrary) edm::LogInfo("ForwardSim") << "\n Using Castor Shower Library \n";
00072 
00073 }
00074 
00075 //=============================================================================================
00076 
00077 CastorSD::~CastorSD() {
00078   if (useShowerLibrary) delete showerLibrary;
00079 }
00080 
00081 //=============================================================================================
00082 
00083 void CastorSD::initRun(){
00084   if (useShowerLibrary) {
00085     // showerLibrary = new CastorShowerLibrary(name, cpv, p);
00086     G4ParticleTable * theParticleTable = G4ParticleTable::GetParticleTable();
00087     showerLibrary->initParticleTable(theParticleTable);
00088     edm::LogInfo("ForwardSim") << "CastorSD::initRun: Using Castor Shower Library \n";
00089   }
00090 }
00091 
00092 //=============================================================================================
00093 
00094 double CastorSD::getEnergyDeposit(G4Step * aStep) {
00095   
00096   float NCherPhot = 0.;
00097 
00098   if (aStep == NULL) {
00099     return 0;
00100   } else {
00101     // preStepPoint information *********************************************
00102   
00103     G4StepPoint*       preStepPoint = aStep->GetPreStepPoint();
00104     G4VPhysicalVolume* currentPV    = preStepPoint->GetPhysicalVolume();
00105     G4LogicalVolume*   currentLV    = currentPV->GetLogicalVolume();
00106     G4String           name   = currentPV->GetName();
00107     std::string        nameVolume;
00108     nameVolume.assign(name,0,4);
00109     
00110 #ifdef debugLog
00111     G4SteppingControl  stepControlFlag = aStep->GetControlFlag();
00112     if (aStep->IsFirstStepInVolume()) 
00113       LogDebug("ForwardSim") << "CastorSD::getEnergyDeposit:"
00114                              << "\n IsFirstStepInVolume " ; 
00115 #endif
00116     
00117     // Get theTrack 
00118     G4Track*        theTrack = aStep->GetTrack();
00119     
00120 #ifdef debugLog
00121     if (useShowerLibrary && currentLV==lvCAST) {
00122           LogDebug("ForwardSim") << "CastorSD::getEnergyDeposit:"
00123                                  << "\n TrackID , ParentID , ParticleName ,"
00124                                  << " eta , phi , z , time ,"
00125                                  << " K , E , Mom " 
00126                                  << "\n  TRACKINFO: " 
00127                                  << theTrack->GetTrackID() 
00128                                  << " , " 
00129                                  << theTrack->GetParentID() 
00130                                  << " , "
00131                                  << theTrack->GetDefinition()->GetParticleName() 
00132                                  << " , "
00133                                  << theTrack->GetPosition().eta() 
00134                                  << " , "
00135                                  << theTrack->GetPosition().phi() 
00136                                  << " , "
00137                                  << theTrack->GetPosition().z() 
00138                                  << " , "
00139                                  << theTrack->GetGlobalTime() 
00140                                  << " , "
00141                                  << theTrack->GetKineticEnergy() 
00142                                  << " , "
00143                                  << theTrack->GetTotalEnergy() 
00144                                  << " , "
00145                                  << theTrack->GetMomentum().mag() ;
00146       if(theTrack->GetTrackID() != 1) 
00147             LogDebug("ForwardSim") << "CastorSD::getEnergyDeposit:"
00148                                    << "\n CurrentStepNumber , TrackID , Particle , VertexPosition ,"
00149                                    << " LogicalVolumeAtVertex , CreatorProcess"
00150                                    << "\n  TRACKINFO2: " 
00151                                    << theTrack->GetCurrentStepNumber() 
00152                                    << " , " 
00153                                    << theTrack->GetTrackID() 
00154                                    << " , "
00155                                    << theTrack->GetDefinition()->GetParticleName() 
00156                                    << " , "
00157                                    << theTrack->GetVertexPosition() 
00158                                    << " , "
00159                                    << theTrack->GetLogicalVolumeAtVertex()->GetName() 
00160                                    << " , " 
00161                                    << theTrack->GetCreatorProcess()->GetProcessName() ;
00162     } // end of if(useShowerLibrary)
00163 #endif
00164     
00165     // if particle moves from interaction point or "backwards (halo)
00166     bool backward = false;
00167     G4ThreeVector  hitPoint = preStepPoint->GetPosition();      
00168     G4ThreeVector  hit_mom  = preStepPoint->GetMomentumDirection();
00169     double zint = hitPoint.z();
00170     double pz   = hit_mom.z();
00171 
00172     // Check if theTrack moves backward 
00173     if (pz * zint < 0.) backward = true;
00174     
00175     // Check that theTrack is above the energy threshold to use Shower Library 
00176     bool aboveThreshold = false;
00177     if(theTrack->GetKineticEnergy() > energyThresholdSL) aboveThreshold = true;
00178     
00179     // Check if theTrack is a muon (if so, DO NOT use Shower Library) 
00180     bool notaMuon = true;
00181     G4int mumPDG  =  13;
00182     G4int mupPDG  = -13;
00183     G4int parCode = theTrack->GetDefinition()->GetPDGEncoding();
00184     if (parCode == mupPDG || parCode == mumPDG ) notaMuon = false;
00185     
00186     // angle condition
00187     double theta_max = M_PI - 3.1305; // angle in radians corresponding to -5.2 eta
00188     double R_mom=sqrt(hit_mom.x()*hit_mom.x() + hit_mom.y()*hit_mom.y());
00189     double theta = atan2(R_mom,std::abs(pz));
00190     bool angleok = false;
00191     if ( theta < theta_max) angleok = true;
00192     
00193     // OkToUse
00194     double R = sqrt(hitPoint.x()*hitPoint.x() + hitPoint.y()*hitPoint.y());
00195     bool dot = false;
00196     if ( zint < -14450. && R < 45.) dot = true;
00197     bool inRange = true;
00198     if ( zint < -14700. || R > 193.) inRange = false;
00199     bool OkToUse = false;
00200     if ( inRange && !dot) OkToUse = true;
00201     
00202     if (useShowerLibrary && aboveThreshold && notaMuon && (!backward) && OkToUse && angleok && currentLV == lvCAST ) {
00203       // Use Castor shower library if energy is above threshold, is not a muon 
00204       // and is not moving backward 
00205       getFromLibrary(aStep);
00206     
00207 #ifdef debugLog
00208       LogDebug("ForwardSim") << " Current logical volume is " << nameVolume ;
00209 #endif
00210       return NCherPhot;
00211     } else {
00212       // Usual calculations
00213       // G4ThreeVector      hitPoint = preStepPoint->GetPosition();     
00214       // G4ThreeVector      hit_mom = preStepPoint->GetMomentumDirection();
00215       G4double           stepl    = aStep->GetStepLength()/cm;
00216       G4double           beta     = preStepPoint->GetBeta();
00217       G4double           charge   = preStepPoint->GetCharge();
00218       //        G4VProcess*        curprocess   = preStepPoint->GetProcessDefinedStep();
00219       //        G4String           namePr   = preStepPoint->GetProcessDefinedStep()->GetProcessName();
00220       //        std::string nameProcess;
00221       //        nameProcess.assign(namePr,0,4);
00222 
00223       //        G4LogicalVolume*   lv    = currentPV->GetLogicalVolume();
00224       //        G4Material*        mat   = lv->GetMaterial();
00225       //        G4double           rad   = mat->GetRadlen();
00226 
00227 
00228       // postStepPoint information *********************************************
00229       G4StepPoint* postStepPoint= aStep->GetPostStepPoint();   
00230       G4VPhysicalVolume* postPV    = postStepPoint->GetPhysicalVolume();
00231       G4String           postname   = postPV->GetName();
00232       std::string        postnameVolume;
00233       postnameVolume.assign(postname,0,4);
00234 
00235       // theTrack information  *************************************************
00236       // G4Track*        theTrack = aStep->GetTrack();   
00237       //G4double        entot    = theTrack->GetTotalEnergy();
00238       G4ThreeVector   vert_mom = theTrack->GetVertexMomentumDirection();
00239 
00240       G4ThreeVector  localPoint = theTrack->GetTouchable()->GetHistory()->
00241         GetTopTransform().TransformPoint(hitPoint);
00242 
00243       // calculations...       *************************************************
00244       float phi = -100.;
00245       if (vert_mom.x() != 0) phi = atan2(vert_mom.y(),vert_mom.x()); 
00246       if (phi < 0.) phi += twopi;
00247       G4String       particleType = theTrack->GetDefinition()->GetParticleName();
00248 #ifdef debugLog
00249       float costheta =vert_mom.z()/sqrt(vert_mom.x()*vert_mom.x()+
00250                                         vert_mom.y()*vert_mom.y()+
00251                                         vert_mom.z()*vert_mom.z());
00252       float theta = acos(std::min(std::max(costheta,float(-1.)),float(1.)));
00253       float eta = -log(tan(theta/2));
00254       G4int          primaryID    = theTrack->GetTrackID();
00255       // *************************************************
00256 
00257 
00258       // *************************************************
00259       double edep   = aStep->GetTotalEnergyDeposit();
00260 #endif
00261 
00262       // *************************************************
00263 
00264 
00265       // *************************************************
00266       // take into account light collection curve for plate
00267       //      double weight = curve_Castor(nameVolume, preStepPoint);
00268       //      double edep   = aStep->GetTotalEnergyDeposit() * weight;
00269       // *************************************************
00270 
00271 
00272       // *************************************************
00273       /*    comments for sensitive volumes:      
00274             C001 ...-... CP06,CBU1 ...-...CALM --- > fibres and bundle 
00275             for first release of CASTOR
00276             CASF  --- > quartz plate  for first and second releases of CASTOR  
00277             GF2Q, GFNQ, GR2Q, GRNQ 
00278             for tests with my own test geometry of HF (on ask of Gavrilov)
00279             C3TF, C4TF - for third release of CASTOR
00280       */
00281       double meanNCherPhot=0;
00282 
00283       if (currentLV == lvC3EF || currentLV == lvC4EF || currentLV == lvC3HF ||
00284           currentLV == lvC4HF) {
00285         //      if(nameVolume == "C3EF" || nameVolume == "C4EF" || nameVolume == "C3HF" || nameVolume == "C4HF") {
00286 
00287         float bThreshold = 0.67;
00288         float nMedium = 1.4925;
00289         //     float photEnSpectrDL = (1./400.nm-1./700.nm)*10000000.cm/nm; /* cm-1  */
00290         //     float photEnSpectrDL = 10714.285714;
00291 
00292         float photEnSpectrDE = 1.24;    /* see below   */
00293         /*     E = 2pi*(1./137.)*(eV*cm/370.)/lambda  =     */
00294         /*       = 12.389184*(eV*cm)/lambda                 */
00295         /*     Emax = 12.389184*(eV*cm)/400nm*10-7cm/nm  = 3.01 eV     */
00296         /*     Emin = 12.389184*(eV*cm)/700nm*10-7cm/nm  = 1.77 eV     */
00297         /*     delE = Emax - Emin = 1.24 eV  --> */
00298         /*   */
00299         /* default for Castor nameVolume  == "CASF" or (C3TF & C4TF)  */
00300 
00301         float thFullRefl = 23.;  /* 23.dergee */
00302         float thFullReflRad = thFullRefl*pi/180.;
00303 
00304         /* default for Castor nameVolume  == "CASF" or (C3TF & C4TF)  */
00305         float thFibDir = 45.;  /* .dergee */
00306         /* for test HF geometry volumes:   
00307            if(nameVolume == "GF2Q" || nameVolume == "GFNQ" ||
00308            nameVolume == "GR2Q" || nameVolume == "GRNQ")
00309            thFibDir = 0.0; // .dergee
00310         */
00311         float thFibDirRad = thFibDir*pi/180.;
00312         /*   */
00313         /*   */
00314 
00315         // at which theta the point is located:
00316         //     float th1    = hitPoint.theta();
00317 
00318         // theta of charged particle in LabRF(hit momentum direction):
00319         float costh =hit_mom.z()/sqrt(hit_mom.x()*hit_mom.x()+
00320                                       hit_mom.y()*hit_mom.y()+
00321                                       hit_mom.z()*hit_mom.z());
00322         if (zint < 0) costh = -costh;
00323         float th = acos(std::min(std::max(costh,float(-1.)),float(1.)));
00324 
00325         // just in case (can do bot use):
00326         if (th < 0.) th += twopi;
00327 
00328 
00329 
00330         // theta of cone with Cherenkov photons w.r.t.direction of charged part.:
00331         float costhcher =1./(nMedium*beta);
00332         float thcher = acos(std::min(std::max(costhcher,float(-1.)),float(1.)));
00333 
00334         // diff thetas of charged part. and quartz direction in LabRF:
00335         float DelFibPart = fabs(th - thFibDirRad);
00336 
00337         // define real distances:
00338         float d = fabs(tan(th)-tan(thFibDirRad));   
00339 
00340         //       float a = fabs(tan(thFibDirRad)-tan(thFibDirRad+thFullReflRad));   
00341         //       float r = fabs(tan(th)-tan(th+thcher));   
00342         
00343         float a = tan(thFibDirRad)+tan(fabs(thFibDirRad-thFullReflRad));   
00344         float r = tan(th)+tan(fabs(th-thcher));   
00345 
00346 
00347         // define losses d_qz in cone of full reflection inside quartz direction
00348         float d_qz;
00349 #ifdef debugLog
00350         float variant;
00351 #endif
00352         if(DelFibPart > (thFullReflRad + thcher) ) {
00353           d_qz = 0.; 
00354 #ifdef debugLog
00355           variant=0.;
00356 #endif
00357         }
00358         // if(d > (r+a) ) {d_qz = 0.; variant=0.;}
00359         else {
00360           if((th + thcher) < (thFibDirRad+thFullReflRad) && 
00361              (th - thcher) > (thFibDirRad-thFullReflRad)) {
00362             d_qz = 1.; 
00363 #ifdef debugLog
00364             variant=1.;
00365 #endif
00366           }
00367           // if((DelFibPart + thcher) < thFullReflRad ) {d_qz = 1.; variant=1.;}
00368           // if((d+r) < a ) {d_qz = 1.; variant=1.;}
00369           else {
00370             if((thFibDirRad + thFullReflRad) < (th + thcher) && 
00371                (thFibDirRad - thFullReflRad) > (th - thcher) ) {
00372               // if((thcher - DelFibPart ) > thFullReflRad ) 
00373               // if((r-d ) > a ) 
00374               d_qz = 0.; 
00375 #ifdef debugLog
00376               variant=2.;
00377 #endif
00378             } else {
00379               // if((thcher + DelFibPart ) > thFullReflRad && 
00380               //    thcher < (DelFibPart+thFullReflRad) ) 
00381               //      {
00382               d_qz = 0.; 
00383 #ifdef debugLog
00384               variant=3.;
00385 #endif
00386 
00387               // use crossed length of circles(cone projection)
00388               // dC1/dC2 : 
00389               float arg_arcos = 0.;
00390               float tan_arcos = 2.*a*d;
00391               if(tan_arcos != 0.) arg_arcos =(r*r-a*a-d*d)/tan_arcos; 
00392               arg_arcos = fabs(arg_arcos);
00393               float th_arcos = acos(std::min(std::max(arg_arcos,float(-1.)),float(1.)));
00394               d_qz = th_arcos/pi/2.;
00395               d_qz = fabs(d_qz);
00396 
00397 
00398 
00399               //            }
00400               //             else
00401               //             {
00402               //                d_qz = 0.; variant=4.;
00403               //#ifdef debugLog
00404               // std::cout <<" ===============>variant 4 information: <===== " <<std::endl;
00405               // std::cout <<" !!!!!!!!!!!!!!!!!!!!!!  variant = " << variant  <<std::endl;
00406               //#endif 
00407               //
00408               //             }
00409             }
00410           }
00411         }
00412 
00413 
00414         // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00415                   
00416         if(charge != 0. && beta > bThreshold )  {
00417                           
00418                           meanNCherPhot = 370.*charge*charge*
00419                           ( 1. - 1./(nMedium*nMedium*beta*beta) )*
00420                           photEnSpectrDE*stepl;
00421                           
00422                           G4int poissNCherPhot = (G4int) G4Poisson(meanNCherPhot);
00423                           
00424                           if(poissNCherPhot < 0) poissNCherPhot = 0;
00425                           
00426                           float effPMTandTransport = 0.19;
00427                           double ReflPower = 0.1;
00428                           double proba = d_qz + (1-d_qz)*ReflPower;
00429                           NCherPhot = poissNCherPhot*effPMTandTransport*proba*0.307;
00430 
00431 
00432 #ifdef debugLog
00433           float thgrad = th*180./pi;
00434           float thchergrad = thcher*180./pi;
00435           float DelFibPartgrad = DelFibPart*180./pi;
00436           LogDebug("ForwardSim") << " ==============================> start all "
00437                                  << "information:<========= \n" << " =====> for "
00438                                  << "test:<===  \n" << " variant = " << variant  
00439                                  << "\n thgrad = " << thgrad  << "\n thchergrad "
00440                                  << "= " << thchergrad  << "\n DelFibPartgrad = "
00441                                  << DelFibPartgrad << "\n d_qz = " << d_qz  
00442                                  << "\n =====> Start Step Information <===  \n"
00443                                  << " ===> calo preStepPoint info <===  \n" 
00444                                  << " hitPoint = " << hitPoint  << "\n"
00445                                  << " hitMom = " << hit_mom  << "\n"
00446                                  << " stepControlFlag = " << stepControlFlag 
00447             // << "\n curprocess = " << curprocess << "\n"
00448             // << " nameProcess = " << nameProcess 
00449                                  << "\n charge = " << charge << "\n"
00450                                  << " beta = " << beta << "\n"
00451                                  << " bThreshold = " << bThreshold << "\n"
00452                                  << " thgrad =" << thgrad << "\n"
00453                                  << " effPMTandTransport=" << effPMTandTransport 
00454             // << "\n volume = " << name 
00455                                  << "\n nameVolume = " << nameVolume << "\n"
00456                                  << " nMedium = " << nMedium << "\n"
00457             //  << " rad length = " << rad << "\n"
00458             //  << " material = " << mat << "\n"
00459                                  << " stepl = " << stepl << "\n"
00460                                  << " photEnSpectrDE = " << photEnSpectrDE <<"\n"
00461                                  << " edep = " << edep << "\n"
00462                                  << " ===> calo theTrack info <=== " << "\n"
00463                                  << " particleType = " << particleType << "\n"
00464                                  << " primaryID = " << primaryID << "\n"
00465                                  << " entot= " << theTrack->GetTotalEnergy() << "\n"
00466                                  << " vert_eta= " << eta  << "\n"
00467                                  << " vert_phi= " << phi << "\n"
00468                                  << " vert_mom= " << vert_mom  << "\n"
00469                                  << " ===> calo hit preStepPointinfo <=== "<<"\n"
00470                                  << " local point = " << localPoint << "\n"
00471                                  << " ==============================> final info"
00472                                  << ":  <=== \n" 
00473                                  << " meanNCherPhot = " << meanNCherPhot << "\n"
00474                                  << " poissNCherPhot = " << poissNCherPhot <<"\n"
00475                                  << " NCherPhot = " << NCherPhot;
00476 #endif 
00477         
00478           // Included by WC
00479           //         std::cout << "\n volume = "         << name 
00480           //              << "\n nameVolume = "     << nameVolume << "\n"
00481           //              << "\n postvolume = "     << postname 
00482           //              << "\n postnameVolume = " << postnameVolume << "\n"
00483           //              << "\n particleType = "   << particleType 
00484           //              << "\n primaryID = "      << primaryID << "\n";
00485 
00486         }
00487       }
00488 
00489 
00490 #ifdef debugLog
00491       LogDebug("ForwardSim") << "CastorSD:: " << nameVolume 
00492         //      << " Light Collection Efficiency " << weight
00493                              << " Weighted Energy Deposit " << edep/MeV 
00494                              << " MeV\n";
00495 #endif
00496       // Temporary member for testing purpose only...
00497       // unit_id = setDetUnitId(aStep);
00498       // if(NCherPhot != 0) std::cout << "\n  UnitID = " << unit_id << "  ;  NCherPhot = " << NCherPhot ;
00499        
00500       return NCherPhot;
00501     }
00502   } 
00503   return 0;
00504 }
00505 
00506 //=======================================================================================
00507 
00508 uint32_t CastorSD::setDetUnitId(G4Step* aStep) {
00509   return (numberingScheme == 0 ? 0 : numberingScheme->getUnitID(aStep));
00510 }
00511 
00512 //=======================================================================================
00513 
00514 void CastorSD::setNumberingScheme(CastorNumberingScheme* scheme) {
00515 
00516   if (scheme != 0) {
00517     edm::LogInfo("ForwardSim") << "CastorSD: updates numbering scheme for " 
00518                                << GetName();
00519     if (numberingScheme) delete numberingScheme;
00520     numberingScheme = scheme;
00521   }
00522 }
00523 
00524 //=======================================================================================
00525 
00526 int CastorSD::setTrackID (G4Step* aStep) {
00527 
00528   theTrack     = aStep->GetTrack();
00529 
00530   double etrack = preStepPoint->GetKineticEnergy();
00531   TrackInformation * trkInfo = (TrackInformation *)(theTrack->GetUserInformation());
00532   int      primaryID = trkInfo->getIDonCaloSurface();
00533   if (primaryID == 0) {
00534 #ifdef debugLog
00535     edm::LogWarning("ForwardSim") << "CastorSD: Problem with primaryID **** set by force "
00536                                   << "to TkID **** " << theTrack->GetTrackID();
00537 #endif
00538     primaryID = theTrack->GetTrackID();
00539   }
00540 
00541   if (primaryID != previousID.trackID())
00542     resetForNewPrimary(preStepPoint->GetPosition(), etrack);
00543 
00544   return primaryID;
00545 }
00546 
00547 //=======================================================================================
00548 
00549 uint32_t CastorSD::rotateUnitID(uint32_t unitID, G4Track* track, CastorShowerEvent shower) {
00550 // ==============================================================
00551 //
00552 //   o   Exploit Castor phi symmetry to return newUnitID for  
00553 //       shower hits based on track phi coordinate 
00554 //
00555 // ==============================================================
00556   
00557   // Get 'track' phi:
00558   float   trackPhi = track->GetPosition().phi(); 
00559   if(trackPhi<0) trackPhi += 2*M_PI ;
00560   // Get phi from primary that gave rise to SL 'shower':
00561   float  showerPhi = shower.getPrimPhi(); 
00562   if(showerPhi<0) showerPhi += 2*M_PI ;
00563   // Delta phi:
00564   
00565   //  Find the OctSector for which 'track' and 'shower' belong
00566   int  trackOctSector = (int) (  trackPhi / (M_PI/4) ) ;
00567   int showerOctSector = (int) ( showerPhi / (M_PI/4) ) ;
00568   
00569   uint32_t  newUnitID;
00570   uint32_t         sec = ( ( unitID>>4 ) & 0xF ) ;
00571   uint32_t  complement = ( unitID & 0xFFFFFF0F ) ;
00572   
00573   // edm::LogInfo("ForwardSim") << "\n CastorSD::rotateUnitID:  " 
00574   //                            << "\n      unitID = " << unitID 
00575   //                            << "\n         sec = " << sec 
00576   //                            << "\n  complement = " << complement ; 
00577   
00578   // Get 'track' z:
00579   float   trackZ = track->GetPosition().z();
00580   
00581   int aux ;
00582   int dSec = 2*(trackOctSector - showerOctSector) ;
00583   // if(trackZ<0)  // Good for revision 1.8 of CastorNumberingScheme
00584   if(trackZ>0)  // Good for revision 1.9 of CastorNumberingScheme
00585   {
00586     int sec1 = sec-dSec;
00587     //    sec -= dSec ;
00588     if(sec1<0) sec1  += 16;
00589     if(sec1>15) sec1 -= 16;
00590     sec = (uint32_t)(sec1);
00591   } else {
00592     if( dSec<0 ) sec += 16 ;
00593     sec += dSec ;
00594     aux  = (int) (sec/16) ;
00595     sec -= aux*16 ;
00596   }
00597   sec  = sec<<4 ;
00598   newUnitID = complement | sec ;
00599   
00600 #ifdef debugLog
00601   if(newUnitID != unitID) {
00602     LogDebug("ForwardSim") << "\n CastorSD::rotateUnitID:  " 
00603                            << "\n     unitID = " << unitID 
00604                            << "\n  newUnitID = " << newUnitID ; 
00605   }
00606 #endif
00607   
00608   return newUnitID ;
00609 
00610 }
00611 
00612 //=======================================================================================
00613 
00614 void CastorSD::getFromLibrary (G4Step* aStep) {
00615 
00617 //
00618 //   Method to get hits from the Shower Library
00619 //
00620 //   CastorShowerEvent hits returned by getShowerHits are used to  
00621 //   replace the full simulation of the shower from theTrack
00622 //    
00623 //   "updateHit" save the Hits to a CaloG4Hit container
00624 //
00626 
00627   preStepPoint  = aStep->GetPreStepPoint(); 
00628   theTrack      = aStep->GetTrack();   
00629   bool ok;
00630   
00631   // ****    Call method to retrieve hits from the ShowerLibrary   ****
00632   CastorShowerEvent hits = showerLibrary->getShowerHits(aStep, ok);
00633 
00634   double etrack    = preStepPoint->GetKineticEnergy();
00635   int    primaryID = setTrackID(aStep);
00636   // int    primaryID = theTrack->GetTrackID();
00637 
00638   // Reset entry point for new primary
00639   posGlobal = preStepPoint->GetPosition();
00640   resetForNewPrimary(posGlobal, etrack);
00641 
00642   // Check whether track is EM or HAD
00643   G4int particleCode = theTrack->GetDefinition()->GetPDGEncoding();
00644   bool isEM , isHAD ;
00645   if (particleCode==emPDG || particleCode==epPDG || particleCode==gammaPDG) {
00646     isEM = true ; isHAD = false;
00647   } else {
00648     isEM = false; isHAD = true ;
00649   }
00650 
00651 #ifdef debugLog
00652   //  edm::LogInfo("ForwardSim") << "\n CastorSD::getFromLibrary:  " 
00653   LogDebug("ForwardSim") << "\n CastorSD::getFromLibrary:  " 
00654                          << hits.getNhit() << " hits for " << GetName() << " from " 
00655                          << theTrack->GetDefinition()->GetParticleName() << " of "
00656                          << preStepPoint->GetKineticEnergy()/GeV << " GeV and trackID " 
00657                          << theTrack->GetTrackID()  ;
00658 #endif
00659 
00660   // Scale to correct energy
00661   double E_track = preStepPoint->GetTotalEnergy() ;
00662   double E_SLhit = hits.getPrimE() * GeV ;
00663   double scale = E_track/E_SLhit ;
00664         
00665         //Non compensation 
00666         if (isHAD){
00667                 scale=scale*non_compensation_factor; // if hadronic extend the scale with the non-compensation factor
00668         } else {
00669                 scale=scale; // if electromagnetic, don't do anything
00670         }
00671         
00672   
00673 /*    double theTrackEnergy = theTrack->GetTotalEnergy() ; 
00674   
00675   if(fabs(theTrackEnergy-E_track)>10.) {
00676     edm::LogInfo("ForwardSim") << "\n            TrackID = " << theTrack->GetTrackID()
00677                                << "\n     theTrackEnergy = " << theTrackEnergy
00678                                << "\n preStepPointEnergy = " << E_track ;
00679     G4TrackVector tsec = *(aStep->GetSecondary());
00680     for (unsigned int kk=0; kk<tsec.size(); kk++) {
00681         edm::LogInfo("ForwardSim") << "CastorSD::getFromLibrary:"
00682                                << "\n tsec[" << kk << "]->GetTrackID() = " 
00683                                << tsec[kk]->GetTrackID() 
00684                                << " with energy " 
00685                                << tsec[kk]->GetTotalEnergy() ;
00686     }
00687   }
00688 */  
00689   //  Loop over hits retrieved from the library
00690   for (unsigned int i=0; i<hits.getNhit(); i++) {
00691     
00692     // Get nPhotoElectrons and set edepositEM / edepositHAD accordingly
00693     double nPhotoElectrons    = hits.getNphotons(i);
00694     // Apply scaling
00695       nPhotoElectrons *= scale ;
00696     if(isEM)  {
00697        // edepositEM  = nPhotoElectrons*GeV; 
00698        edepositEM  = nPhotoElectrons; 
00699        edepositHAD = 0.;                 
00700     } else if(isHAD) {
00701        edepositEM  = 0.;                  
00702        edepositHAD = nPhotoElectrons;
00703        // edepositHAD = nPhotoElectrons*GeV;
00704     }
00705     
00706     // Get hit position and time
00707     double                time = hits.getTime(i);
00708     //    math::XYZPoint    position = hits.getHitPosition(i);
00709     
00710     // Get hit detID
00711     unsigned int        unitID = hits.getDetID(i);
00712     
00713     // Make the detID "rotation" from one sector to another taking into account the 
00714     // sectors of the impinging particle (theTrack) and of the particle that produced 
00715     // the 'hits' retrieved from shower library   
00716     unsigned int rotatedUnitID = rotateUnitID(unitID , theTrack , hits);
00717     currentID.setID(rotatedUnitID, time, primaryID, 0);
00718     // currentID.setID(unitID, time, primaryID, 0);
00719    
00720     // check if it is in the same unit and timeslice as the previous one
00721     if (currentID == previousID) {
00722       updateHit(currentHit);
00723     } else {
00724       if (!checkHit()) currentHit = createNewHit();
00725     }
00726   }  //  End of loop over hits
00727 
00728   //Now kill the current track
00729   if (ok) {
00730     theTrack->SetTrackStatus(fStopAndKill);
00731 #ifdef debugLog
00732     LogDebug("ForwardSim") << "CastorSD::getFromLibrary:"
00733                            << "\n \"theTrack\" with TrackID() = " 
00734                            << theTrack->GetTrackID() 
00735                            << " and with energy " 
00736                            << theTrack->GetTotalEnergy()
00737                            << " has been set to be killed" ;
00738 #endif
00739     G4TrackVector tv = *(aStep->GetSecondary());
00740     for (unsigned int kk=0; kk<tv.size(); kk++) {
00741       if (tv[kk]->GetVolume() == preStepPoint->GetPhysicalVolume()) {
00742         tv[kk]->SetTrackStatus(fStopAndKill);
00743 #ifdef debugLog
00744         LogDebug("ForwardSim") << "CastorSD::getFromLibrary:"
00745                                << "\n tv[" << kk << "]->GetTrackID() = " 
00746                                << tv[kk]->GetTrackID() 
00747                                << " with energy " 
00748                                << tv[kk]->GetTotalEnergy()
00749                                << " has been set to be killed" ;
00750 #endif
00751       }
00752     }
00753   }
00754 }
00755
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/SimG4CMS/Forward/src/CastorSD.cc
