CMS 3D CMS Logo

/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/SimG4CMS/Forward/src/CastorSD.cc

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