CMSSW Reference Manual

00001 // -*- C++ -*-
00002 //
00003 // Package:     Forward
00004 // Class  :     TotemSD
00005 //
00006 // Implementation:
00007 //     <Notes on implementation>
00008 //
00009 // Original Author: 
00010 //         Created:  Tue May 16 10:14:34 CEST 2006
00011 // $Id: TotemSD.cc,v 1.4 2007/11/20 12:37:21 fabiocos Exp $
00012 //
00013 
00014 // system include files
00015 
00016 // user include files
00017 #include "FWCore/Framework/interface/ESHandle.h"
00018 #include "FWCore/Framework/interface/EventSetup.h"
00019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00020 
00021 #include "SimG4Core/Notification/interface/TrackInformation.h"
00022 #include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
00023 #include "SimG4Core/Physics/interface/G4ProcessTypeEnumerator.h"
00024  
00025 #include "SimDataFormats/TrackingHit/interface/UpdatablePSimHit.h"
00026 #include "SimDataFormats/SimHitMaker/interface/TrackingSlaveSD.h"
00027 
00028 #include "SimG4CMS/Forward/interface/TotemSD.h"
00029 #include "SimG4CMS/Forward/interface/TotemNumberMerger.h"
00030 #include "SimG4CMS/Forward/interface/TotemT1NumberingScheme.h"
00031 #include "SimG4CMS/Forward/interface/TotemT2NumberingSchemeGem.h"
00032 #include "SimG4CMS/Forward/interface/TotemRPNumberingScheme.h"
00033 
00034 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
00035 #include "DataFormats/GeometryVector/interface/LocalVector.h"
00036 
00037 #include "G4SDManager.hh"
00038 #include "G4Step.hh"
00039 #include "G4Track.hh"
00040 #include "G4VProcess.hh"
00041 
00042 //
00043 // constructors and destructor
00044 //
00045 TotemSD::TotemSD(std::string name, const DDCompactView & cpv,
00046                  SensitiveDetectorCatalog & clg, 
00047                  edm::ParameterSet const & p, const SimTrackManager* manager) :
00048   SensitiveTkDetector(name, cpv, clg, p), numberingScheme(0), name(name),
00049   hcID(-1), theHC(0), theManager(manager), currentHit(0), theTrack(0), 
00050   currentPV(0), unitID(0),  previousUnitID(0), preStepPoint(0), 
00051   postStepPoint(0), eventno(0){
00052 
00053   //Add Totem Sentitive Detector Names
00054   collectionName.insert(name);
00055 
00056   //Parameters
00057   edm::ParameterSet m_p = p.getParameter<edm::ParameterSet>("TotemSD");
00058   int verbn = m_p.getUntrackedParameter<int>("Verbosity");
00059  
00060   SetVerboseLevel(verbn);
00061   LogDebug("ForwardSim") 
00062     << "*******************************************************\n"
00063     << "*                                                     *\n"
00064     << "* Constructing a TotemSD  with name " << name << "\n"
00065     << "*                                                     *\n"
00066     << "*******************************************************";
00067 
00068   slave  = new TrackingSlaveSD(name);
00069 
00070   //
00071   // Now attach the right detectors (LogicalVolumes) to me
00072   //
00073   std::vector<std::string> lvNames = clg.logicalNames(name);
00074   this->Register();
00075   for (std::vector<std::string>::iterator it=lvNames.begin();  
00076        it !=lvNames.end(); it++) {
00077     this->AssignSD(*it);
00078     edm::LogInfo("ForwardSim") << "TotemSD : Assigns SD to LV " << (*it);
00079   }
00080 
00081   if      (name == "TotemHitsT1") {
00082     numberingScheme = dynamic_cast<TotemVDetectorOrganization*>(new TotemT1NumberingScheme(1));
00083   } else if (name == "TotemHitsT2Si") {
00084     numberingScheme = dynamic_cast<TotemVDetectorOrganization*>(new TotemT2NumberingSchemeGem(3));
00085   } else if (name == "TotemHitsT2Gem") {
00086     numberingScheme = dynamic_cast<TotemVDetectorOrganization*>(new TotemT2NumberingSchemeGem(4));
00087   } else if (name == "TotemHitsRP") {
00088     numberingScheme = dynamic_cast<TotemVDetectorOrganization*>(new TotemRPNumberingScheme(3));
00089   } else {
00090     edm::LogWarning("ForwardSim") << "TotemSD: ReadoutName not supported\n";
00091   }
00092   
00093   edm::LogInfo("ForwardSim") << "TotemSD: Instantiation completed";
00094 } 
00095 
00096 TotemSD::~TotemSD() { 
00097   if (slave)           delete slave; 
00098   if (numberingScheme) delete numberingScheme;
00099 }
00100 
00101 bool TotemSD::ProcessHits(G4Step * aStep, G4TouchableHistory * ) {
00102 
00103   if (aStep == NULL) {
00104     return true;
00105   } else {
00106     GetStepInfo(aStep);
00107     if (HitExists() == false && edeposit>0.) { 
00108       CreateNewHit();
00109       return true;
00110     }
00111     if (HitExists() == false && (((unitID==1111 || unitID==2222) && 
00112                                   ParentId==0 && ParticleType==2212))) { 
00113       CreateNewHitEvo();
00114       return true;
00115     }
00116   }
00117   return true;
00118 }
00119 
00120 uint32_t TotemSD::setDetUnitId(G4Step * aStep) { 
00121 
00122   return (numberingScheme == 0 ? 0 : numberingScheme->GetUnitID(aStep));
00123 }
00124 
00125 void TotemSD::Initialize(G4HCofThisEvent * HCE) { 
00126 
00127   LogDebug("ForwardSim") << "TotemSD : Initialize called for " << name;
00128 
00129   theHC = new TotemG4HitCollection(name, collectionName[0]);
00130   if (hcID<0) 
00131     hcID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
00132   HCE->AddHitsCollection(hcID, theHC);
00133 
00134   tsID   = -2;
00135   primID = -2;
00136 }
00137 
00138 void TotemSD::EndOfEvent(G4HCofThisEvent* ) {
00139 
00140   // here we loop over transient hits and make them persistent
00141   for (int j=0; j<theHC->entries() && j<15000; j++) {
00142     TotemG4Hit* aHit = (*theHC)[j];
00143 #ifdef ddebug
00144     LogDebug("ForwardSim") << "HIT NUMERO " << j << "unit ID = "
00145                            << aHit->getUnitID() << "\n"
00146                            << "               " << "enrty z " 
00147                            << aHit->getEntry().z() << "\n"
00148                            << "               " << "theta   " 
00149                            << aHit->getThetaAtEntry() << "\n";
00150 #endif
00151     Local3DPoint theExitPoint(0,0,0);
00152     Local3DPoint Entrata(aHit->getEntry().x(),
00153                          aHit->getEntry().y(),
00154                          aHit->getEntry().z());
00155     slave->processHits(PSimHit(Entrata,theExitPoint,
00156                                aHit->getPabs(), aHit->getTof(),
00157                                aHit->getEnergyLoss(), aHit->getParticleType(),
00158                                aHit->getUnitID(), aHit->getTrackID(),
00159                                aHit->getThetaAtEntry(),aHit->getPhiAtEntry()));
00160 
00161   }
00162   Summarize();
00163 }
00164 
00165 void TotemSD::clear() {
00166 } 
00167 
00168 void TotemSD::DrawAll() {
00169 } 
00170 
00171 void TotemSD::PrintAll() {
00172   LogDebug("ForwardSim") << "TotemSD: Collection " << theHC->GetName();
00173   theHC->PrintAllHits();
00174 } 
00175 
00176 void TotemSD::fillHits(edm::PSimHitContainer& c, std::string n) {
00177   if (slave->name() == n) c=slave->hits();
00178 }
00179 
00180 void TotemSD::update (const BeginOfEvent * i) {
00181   LogDebug("ForwardSim") << " Dispatched BeginOfEvent for " << GetName()
00182                        << " !" ;
00183    clearHits();
00184    eventno = (*i)()->GetEventID();
00185 }
00186 
00187 void TotemSD::update (const ::EndOfEvent*) {
00188 }
00189 
00190 void TotemSD::clearHits(){
00191   slave->Initialize();
00192 }
00193 
00194 G4ThreeVector TotemSD::SetToLocal(G4ThreeVector global) {
00195 
00196   G4ThreeVector       localPoint;
00197   const G4VTouchable* touch= preStepPoint->GetTouchable();
00198   localPoint = touch->GetHistory()->GetTopTransform().TransformPoint(global);
00199   return localPoint;  
00200 }
00201 
00202 void TotemSD::GetStepInfo(G4Step* aStep) {
00203   
00204   preStepPoint = aStep->GetPreStepPoint(); 
00205   postStepPoint= aStep->GetPostStepPoint(); 
00206   theTrack     = aStep->GetTrack();   
00207   //Local3DPoint theEntryPoint = SensitiveDetector::InitialStepPosition(aStep,LocalCoordinates);  
00208   //Local3DPoint theExitPoint  = SensitiveDetector::FinalStepPosition(aStep,LocalCoordinates);
00209   hitPoint     = preStepPoint->GetPosition();   
00210   currentPV    = preStepPoint->GetPhysicalVolume();
00211 
00212   // double weight = 1; 
00213   G4String name = currentPV->GetName();
00214   name.assign(name,0,4);
00215   G4String particleType = theTrack->GetDefinition()->GetParticleName();
00216   edeposit = aStep->GetTotalEnergyDeposit();
00217   
00218   tSlice    = (postStepPoint->GetGlobalTime() )/nanosecond;
00219   tSliceID  = (int) tSlice;
00220   unitID    = setDetUnitId(aStep);
00221 #ifdef debug
00222   LogDebug("ForwardSim") << "UNITa " << unitID;
00223 #endif
00224   primaryID = theTrack->GetTrackID();
00225 
00226 
00227   Posizio = hitPoint;
00228   Pabs    = aStep->GetPreStepPoint()->GetMomentum().mag()/GeV;
00229   Tof     = aStep->GetPostStepPoint()->GetGlobalTime()/nanosecond;
00230    
00231   Eloss   = aStep->GetTotalEnergyDeposit()/GeV;
00232   ParticleType = theTrack->GetDefinition()->GetPDGEncoding();      
00233 
00234   ThetaAtEntry = aStep->GetPreStepPoint()->GetPosition().theta()/deg;
00235   PhiAtEntry   = aStep->GetPreStepPoint()->GetPosition().phi()/deg;
00236 
00237   ParentId = theTrack->GetParentID();
00238   Vx = theTrack->GetVertexPosition().x();
00239   Vy = theTrack->GetVertexPosition().y();
00240   Vz = theTrack->GetVertexPosition().z();
00241 }
00242 
00243 bool TotemSD::HitExists() {
00244    
00245   if (primaryID<1) {
00246     edm::LogWarning("ForwardSim") << "***** TotemSD error: primaryID = " 
00247                                   << primaryID
00248                                   << " maybe detector name changed";
00249   }
00250    
00251   // Update if in the same detector, time-slice and for same track   
00252   //  if (primaryID == primID && tSliceID == tsID && unitID==previousUnitID) {
00253   if (tSliceID == tsID && unitID==previousUnitID) {
00254     UpdateHit();
00255     return true;
00256   }
00257    
00258   // Reset entry point for new primary
00259   if (primaryID != primID)
00260     ResetForNewPrimary();
00261    
00262   //look in the HitContainer whether a hit with the same primID, unitID,
00263   //tSliceID already exists:
00264    
00265   bool found = false;
00266 
00267   for (int j=0; j<theHC->entries()&&!found; j++) {
00268     TotemG4Hit* aPreviousHit = (*theHC)[j];
00269     if (aPreviousHit->getTrackID()     == primaryID &&
00270         aPreviousHit->getTimeSliceID() == tSliceID  &&
00271         aPreviousHit->getUnitID()      == unitID       ) {
00272       currentHit = aPreviousHit;
00273       found      = true;
00274     }
00275   }          
00276 
00277   if (found) {
00278     UpdateHit();
00279     return true;
00280   } else {
00281     return false;
00282   }    
00283 }
00284 
00285 void TotemSD::CreateNewHit() {
00286 
00287 #ifdef debug
00288   LogDebug("ForwardSim") << "TotemSD CreateNewHit for"
00289                          << " PV "     << currentPV->GetName()
00290                          << " PVid = " << currentPV->GetCopyNo()
00291                          << " MVid = " << currentPV->GetMother()->GetCopyNo()
00292                          << " Unit "   << unitID << "\n"
00293                          << " primary "    << primaryID
00294                          << " time slice " << tSliceID 
00295                          << " For Track  " << theTrack->GetTrackID()
00296                          << " which is a " 
00297                          << theTrack->GetDefinition()->GetParticleName();
00298            
00299   if (theTrack->GetTrackID()==1) {
00300     LogDebug("ForwardSim") << " of energy "     << theTrack->GetTotalEnergy();
00301   } else {
00302     LogDebug("ForwardSim") << " daughter of part. " << theTrack->GetParentID();
00303   }
00304 
00305   cout  << " and created by " ;
00306   if (theTrack->GetCreatorProcess()!=NULL)
00307     LogDebug("ForwardSim") << theTrack->GetCreatorProcess()->GetProcessName() ;
00308   else 
00309     LogDebug("ForwardSim") << "NO process";
00310 #endif          
00311     
00312 
00313   currentHit = new TotemG4Hit;
00314   currentHit->setTrackID(primaryID);
00315   currentHit->setTimeSlice(tSlice);
00316   currentHit->setUnitID(unitID);
00317   currentHit->setIncidentEnergy(incidentEnergy);
00318 
00319   currentHit->setPabs(Pabs);
00320   currentHit->setTof(Tof);
00321   currentHit->setEnergyLoss(Eloss);
00322   currentHit->setParticleType(ParticleType);
00323   currentHit->setThetaAtEntry(ThetaAtEntry);
00324   currentHit->setPhiAtEntry(PhiAtEntry);
00325 
00326   currentHit->setEntry(Posizio.x(),Posizio.y(),Posizio.z());
00327 
00328   currentHit->setParentId(ParentId);
00329   currentHit->setVx(Vx);
00330   currentHit->setVy(Vy);
00331   currentHit->setVz(Vz);
00332 
00333   UpdateHit();
00334   
00335   StoreHit(currentHit);
00336 }        
00337 
00338 void TotemSD::CreateNewHitEvo() {
00339 
00340 // LogDebug("ForwardSim") << "INSIDE CREATE NEW HIT EVO ";
00341 
00342   currentHit = new TotemG4Hit;
00343   currentHit->setTrackID(primaryID);
00344   currentHit->setTimeSlice(tSlice);
00345   currentHit->setUnitID(unitID);
00346   currentHit->setIncidentEnergy(incidentEnergy);
00347 
00348   currentHit->setPabs(Pabs);
00349   currentHit->setTof(Tof);
00350   currentHit->setEnergyLoss(Eloss);
00351   currentHit->setParticleType(ParticleType);
00352   currentHit->setThetaAtEntry(ThetaAtEntry);
00353   currentHit->setPhiAtEntry(PhiAtEntry);
00354 
00355   //  LogDebug("ForwardSim") << Posizio.x() << " " << Posizio.y() << " " << Posizio.z();
00356 
00357   currentHit->setParentId(ParentId);
00358   currentHit->setVx(Vx);
00359   currentHit->setVy(Vy);
00360   currentHit->setVz(Vz);
00361 
00362   G4ThreeVector _PosizioEvo;
00363   int flagAcc=0;
00364   _PosizioEvo=PosizioEvo(Posizio,Vx,Vy,Vz,Pabs,flagAcc);
00365 
00366   if(flagAcc==1){
00367     currentHit->setEntry(_PosizioEvo.x(),_PosizioEvo.y(),_PosizioEvo.z());
00368 
00369     // if(flagAcc==1)
00370     UpdateHit();
00371   
00372     StoreHit(currentHit);
00373   }
00374   // LogDebug("ForwardSim") << "STORED HIT IN: " << unitID;
00375 }        
00376  
00377 G4ThreeVector TotemSD::PosizioEvo(G4ThreeVector Pos, double vx, double vy,
00378                                   double vz, double pabs, int& accettanza) {
00379   accettanza=0;
00380   //Pos.xyz() in mm
00381   G4ThreeVector PosEvo; 
00382   double ThetaX=atan((Pos.x()-vx)/(Pos.z()-vz));                 
00383   double ThetaY=atan((Pos.y()-vy)/(Pos.z()-vz));                
00384   double X_at_0 =(vx-((Pos.x()-vx)/(Pos.z()-vz))*vz)/1000.;   
00385   double Y_at_0 =(vy-((Pos.y()-vy)/(Pos.z()-vz))*vz)/1000.;  
00386   
00387   //  double temp_evo_X;
00388   //  double temp_evo_Y;
00389   //  double temp_evo_Z;
00390   //  temp_evo_Z = fabs(Pos.z())/Pos.z()*220000.; 
00391  
00392   //csi=-dp/d
00393   double csi = fabs((7000.-pabs)/7000.);
00394 
00395   // all in m 
00396   const int no_rp=4;
00397   double x_par[no_rp+1];
00398   double y_par[no_rp+1];
00399   //rp z position
00400   double rp[no_rp]={141.,149.,198.,220.};
00401   //{lx0,mlx} for each rp; Lx=lx0+mlx*csi
00402   double leffx[][2]={{122.5429,-46.9312},{125.4194,-49.1849},{152.6,-81.157},{98.8914,-131.8390}};
00403   //{ly0,mly} for each rp; Ly=ly0+mly*csi
00404   double leffy[][2]={{124.2314,-55.4852},{127.7825,-57.4503},{179.455,-76.274},{273.0931,-40.4626}};
00405   //{vx0,mvx0} for each rp; vx=vx0+mvx*csi
00406   double avx[][2]={{0.515483,-1.0123},{0.494122,-1.0534},{0.2217,-1.483},{0.004633,-1.0719}};
00407   //{vy0,mvy0} for each rp; vy=vy0+mvy*csi
00408   double avy[][2]={{0.371418,-1.6327},{0.349035,-1.6955},{0.0815,-2.59},{0.007592,-4.0841}};                
00409   //{D0,md,a,b} for each rp; D=D0+(md+a*thetax)*csi+b*thetax
00410   double ddx[][4]= {{-0.082336,-0.092513,112.3436,-82.5029},{-0.086927,-0.097670,114.9513,-82.9835},
00411                     {-0.092117,-0.0915,180.6236,-82.443},{-0.050470,0.058837,208.1106,20.8198}};
00412   // {10sigma_x+0.5mm,10sigma_y+0.5mm}
00413   double detlim[][2]={{0,0},{0.0028,0.0021},{0,0},{0.0008,0.0013}};   
00414   //{rmax,dmax}
00415   double pipelim[][2]={{0.026,0.026},{0.04,0.04},{0.0226,0.0177},{0.04,0.04}};
00416   
00417   
00418   for(int j=0; j<no_rp ; j++)  { 
00419     //y=Ly*thetay+vy*y0
00420     //x=Lx*thetax+vx*x0-csi*D   
00421     y_par[j]=ThetaY*(leffy[j][0]+leffy[j][1]*csi)+(avy[j][0]+avy[j][1]*csi)*Y_at_0;
00422     x_par[j]=ThetaX*(leffx[j][0]+leffx[j][1]*csi)+(avx[j][0]+avx[j][1]*csi)*X_at_0-
00423       csi*(ddx[j][0]+(ddx[j][1]+ddx[j][2]*ThetaX)*csi+ddx[j][3]*ThetaX);
00424   }
00425    
00426    
00427   //pass TAN@141
00428   if (fabs(y_par[0])<pipelim[0][1] && sqrt((y_par[0]*y_par[0])+(x_par[0]*x_par[0]))<pipelim[0][0])  {
00429     //pass 149
00430     if ((sqrt((y_par[1]*y_par[1])+(x_par[1]*x_par[1]))<pipelim[1][0]) &&
00431         (fabs(y_par[1])>detlim[1][1] || x_par[1]>detlim[1][0]))  {
00432       accettanza = 1;
00433     }
00434   }
00435 
00436       
00437   //pass TAN@141
00438   if (fabs(y_par[0])<pipelim[0][1] && sqrt((y_par[0])*(y_par[0])+(x_par[0])*(x_par[0]))<pipelim[0][0]) {
00439     //pass Q5@198
00440     if (fabs(y_par[2])<pipelim[2][1] && sqrt((y_par[2]*y_par[2])+(x_par[2]*x_par[2]))<pipelim[2][0]) {
00441       //pass 220
00442       if ((sqrt((y_par[3]*y_par[3])+(x_par[3]*x_par[3]))<pipelim[3][0]) &&
00443           (fabs(y_par[3])>detlim[3][1] || x_par[3]>detlim[3][0])) {
00444         accettanza = 1;
00445         
00446         PosEvo.setX(1000*x_par[3]);
00447         PosEvo.setY(1000*y_par[3]);
00448         PosEvo.setZ(1000*rp[3]);          
00449         if(Pos.z()<vz)PosEvo.setZ(-1000*rp[3]);
00450       }
00451     }
00452     
00453   }
00454 /*
00455   LogDebug("ForwardSim") << "\n"
00456                          << "ACCETTANZA: "<<accettanza << "\n" 
00457                          << "CSI: "<< csi << "\n"
00458                          << "Theta_X: " << ThetaX << "\n"
00459                          << "Theta_Y: " << ThetaY << "\n"
00460                          << "X_at_0: "<< X_at_0 << "\n"
00461                          << "Y_at_0: "<< Y_at_0 << "\n" 
00462                          << "x_par[3]: "<< x_par[3] << "\n"
00463                          << "y_par[3]: "<< y_par[3] << "\n"
00464                          << "pos " << Pos.x() << " " << Pos.y() << " " 
00465                          << Pos.z() << "\n" << "V "<< vx << " " << vy << " "
00466                          << vz << "\n"
00467 */
00468 // --------------
00469   return PosEvo;
00470 }
00471  
00472 
00473 void TotemSD::UpdateHit() {
00474   //
00475   if (Eloss > 0.) {
00476     //  currentHit->addEnergyDeposit(edepositEM,edepositHAD);
00477 
00478 #ifdef debug
00479     LogDebug("ForwardSim") << "G4TotemT1SD updateHit: add eloss " << Eloss 
00480                            << "\nCurrentHit=" << currentHit
00481                            << ", PostStepPoint=" 
00482                            << postStepPoint->GetPosition();
00483 #endif
00484 
00485     currentHit->setEnergyLoss(Eloss);
00486   }  
00487   //  if(PostStepPoint->GetPhysicalVolume() != CurrentPV){
00488   //  currentHit->setExitPoint(SetToLocal(postStepPoint->GetPosition()));
00489   // Local3DPoint exit=currentHit->exitPoint();
00490 /*
00491 #ifdef debug
00492   LogDebug("ForwardSim") << "G4TotemT1SD updateHit: exit point " 
00493                          << exit.x() << " " << exit.y() << " " << exit.z();
00494 //  LogDebug("ForwardSim") << "Energy deposit in Unit " << unitID << " em " << edepositEM/MeV
00495 // << " hadronic " << edepositHAD/MeV << " MeV";
00496 #endif
00497 */
00498 
00499   // buffer for next steps:
00500   tsID           = tSliceID;
00501   primID         = primaryID;
00502   previousUnitID = unitID;
00503 }
00504 
00505 void TotemSD::StoreHit(TotemG4Hit* hit) {
00506 
00507   if (primID<0) return;
00508   if (hit == 0 ) {
00509     edm::LogWarning("ForwardSim") << "TotemSD: hit to be stored is NULL !!";
00510     return;
00511   }
00512 
00513   theHC->insert( hit );
00514 }
00515 
00516 void TotemSD::ResetForNewPrimary() {
00517   
00518   entrancePoint  = SetToLocal(hitPoint);
00519   incidentEnergy = preStepPoint->GetKineticEnergy();
00520 }
00521 
00522 void TotemSD::Summarize() {
00523 }