00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "BesEmcSD.hh"
00011
00012 #include "BesEmcHit.hh"
00013 #include "BesEmcConstruction.hh"
00014 #include "BesEmcGeometry.hh"
00015 #include "G4VPhysicalVolume.hh"
00016 #include "G4Step.hh"
00017 #include "G4VTouchable.hh"
00018 #include "G4TouchableHistory.hh"
00019 #include "G4SDManager.hh"
00020 #include "G4Event.hh"
00021 #include "G4ios.hh"
00022 #include "BesSensitiveManager.hh"
00023 #include "BesTruthTrack.hh"
00024 #include <sstream>
00025 #include "Identifier/EmcID.h"
00026
00027
00028
00029 BesEmcSD::BesEmcSD(G4String name,
00030 BesEmcConstruction* det,
00031 BesEmcGeometry* besEMCGeometry)
00032 :BesSensitiveDetector(name)
00033 ,Detector(det),fBesEmcGeometry(besEMCGeometry)
00034 {
00035 collectionName.insert("BesEmcHitsCollection");
00036 collectionName.insert("BesEmcHitsList");
00037 collectionName.insert("BesEmcTruthHitsList");
00038 HitID = new G4int[20000];
00039 }
00040
00041
00042
00043 BesEmcSD::~BesEmcSD()
00044 {
00045 delete [] HitID;
00046 }
00047
00048
00049
00050 void BesEmcSD::Initialize(G4HCofThisEvent*)
00051 {
00052 CalCollection = new BesEmcHitsCollection
00053 (SensitiveDetectorName,collectionName[0]);
00054 for (G4int j=0;j<20000;j++)
00055 {
00056 HitID[j] = -1;
00057 }
00058 nHit=0;
00059 }
00060
00061 void BesEmcSD::BeginOfTruthEvent(const G4Event* )
00062 {
00063 CalList = new BesEmcHitsCollection
00064 (SensitiveDetectorName,collectionName[1]);
00065 CalTruthList = new BesEmcTruthHitsCollection
00066 (SensitiveDetectorName,collectionName[2]);
00067 m_trackIndex = -99;
00068 }
00069
00070 void BesEmcSD::EndOfTruthEvent(const G4Event* evt)
00071 {
00072 static G4int HLID=-1;
00073 if(HLID<0)
00074 HLID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[1]);
00075 G4HCofThisEvent* HCE = evt->GetHCofThisEvent();
00076 HCE->AddHitsCollection(HLID,CalList);
00077
00078 static G4int HTID=-1;
00079 if(HTID<0)
00080 HTID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[2]);
00081 HCE->AddHitsCollection(HTID,CalTruthList);
00082 }
00083
00084
00085 G4bool BesEmcSD::ProcessHits(G4Step* aStep,G4TouchableHistory*)
00086 {
00087
00088 G4double edep = aStep->GetTotalEnergyDeposit();
00089 G4double stepl = aStep->GetStepLength();
00090 if ((edep==0.)&&(stepl==0.)) return false;
00091
00092 G4TouchableHistory* theTouchable
00093 = (G4TouchableHistory*)(aStep->GetPreStepPoint()->GetTouchable());
00094 G4VPhysicalVolume* physVol = theTouchable->GetVolume();
00095
00096 if(physVol->GetName().contains("physicalCrystal"))
00097 {
00098 PartId=1;
00099 CryNumberPhi=theTouchable->GetReplicaNumber(2);
00100 CryNumberTheta=theTouchable->GetReplicaNumber(1);
00101 }
00102 else if(physVol->GetName().contains("logicalCrystal"))
00103 {
00104 PartId=1;
00105 std::istringstream thetaBuf((theTouchable->GetVolume(1)->GetName()).substr(16,2));
00106 thetaBuf >> CryNumberTheta ;
00107 CryNumberPhi = 308-theTouchable->GetCopyNumber(2);
00108 }
00109 else if(physVol->GetName().contains("physicalEndCrystal"))
00110 {
00111 PartId=theTouchable->GetReplicaNumber(3);
00112 G4int endSector=theTouchable->GetReplicaNumber(2);
00113 G4int endNb=theTouchable->GetReplicaNumber(0);
00114 ComputeThetaPhi(PartId,endSector,ComputeEndCopyNb(endNb));
00115 }
00116 else if(physVol->GetName().contains("logicalEndCrystal"))
00117 {
00118 PartId=2-2*theTouchable->GetCopyNumber(3);
00119 G4int endSector=theTouchable->GetCopyNumber(2);
00120 G4int endNb,endNbGDML;
00121 std::istringstream thetaBuf((theTouchable->GetVolume(0)->GetName()).substr(20,2));
00122 thetaBuf >> endNb ;
00123 endNbGDML=ComputeEndCopyNb(endNb);
00124 G4int endSectorGDML=EndPhiNumberForGDML(endSector);
00125 ComputeThetaPhi(PartId,endSectorGDML,endNbGDML);
00126 }
00127 else if(physVol->GetName().contains("physicalPD"))
00128 {
00129 edep*=50.;
00130 PartId=1;
00131 CryNumberPhi=theTouchable->GetReplicaNumber(2);
00132 CryNumberTheta=theTouchable->GetReplicaNumber(1);
00133 }
00134 else if(physVol->GetName().contains("logicalPD"))
00135 {
00136 edep*=50.;
00137 PartId=1;
00138 G4int nb=theTouchable->GetCopyNumber(1);
00139 if(nb==216) {
00140 CryNumberTheta=0;
00141 } else {
00142 CryNumberTheta = 43-(nb-2)/5;
00143 }
00144 CryNumberPhi = 308-theTouchable->GetCopyNumber(2);
00145 }
00146
00147 if (verboseLevel>1)
00148 G4cout << "(Check ID)New EMC Hit on crystal: "
00149 << CryNumberPhi << "(phi)"
00150 << CryNumberTheta << "(theta)"
00151 << " and the volume is " << physVol->GetName()
00152 << G4endl;
00153
00154
00155 G4int trackId = aStep->GetTrack()->GetTrackID();
00156
00157 BesEmcHit* calHit = new BesEmcHit();
00158 calHit->SetTrackIndex(-99);
00159 calHit->SetG4Index(trackId);
00160 calHit->SetEdepCrystal(edep);
00161 calHit->SetTrakCrystal(stepl);
00162 calHit->SetPosCrystal(aStep->GetPreStepPoint()->GetPosition());
00163 calHit->SetTimeCrystal(aStep->GetPreStepPoint()->GetGlobalTime());
00164 calHit->SetMomentum(aStep->GetPreStepPoint()->GetMomentum());
00165 calHit->SetNumCrystal(PartId,CryNumberTheta,CryNumberPhi);
00166
00167 if(edep>0&&nHit<20000)
00168 HitID[nHit]=CalCollection->insert(calHit)-1;
00169 else
00170 {
00171 delete calHit;
00172
00173
00174 }
00175 nHit++;
00176
00177
00178 if(CalList)
00179 {
00180 G4int trackIndex, g4TrackId;
00181 GetCurrentTrackIndex(trackIndex, g4TrackId);
00182 calHit->SetTrackIndex(trackIndex);
00183
00184 if(m_trackIndex != trackIndex)
00185 {
00186 m_trackIndex = trackIndex;
00187
00188 G4int flag=1;
00189 G4int pdg = abs(aStep->GetTrack()->GetDefinition()->GetPDGEncoding());
00190 if(pdg==12 || pdg==14 || pdg==16) {
00191 flag=0;
00192 }
00193
00194 if(flag && aStep->GetTrack()->GetTrackID()==g4TrackId )
00195 {
00196 BesEmcHit* truHit = new BesEmcHit();
00197 truHit->SetTrackIndex(trackIndex);
00198 truHit->SetG4Index(trackId);
00199 truHit->SetEdepCrystal(edep);
00200 truHit->SetTrakCrystal(stepl);
00201 truHit->SetPosCrystal(aStep->GetPreStepPoint()->GetPosition());
00202 truHit->SetTimeCrystal(aStep->GetPreStepPoint()->GetGlobalTime());
00203 truHit->SetMomentum(aStep->GetPreStepPoint()->GetMomentum());
00204 truHit->SetNumCrystal(PartId,CryNumberTheta,CryNumberPhi);
00205 CalList->insert(truHit);
00206 }
00207 }
00208
00209 else if(m_trackIndex == trackIndex)
00210 {
00211 G4int length = CalList->entries();
00212 if(length>=1)
00213 {
00214 for(G4int i=0;i<length;i++)
00215 {
00216 BesEmcHit* oldHit =(*CalList)[i];
00217 if(oldHit->GetTrackIndex()==trackIndex)
00218 {
00219 oldHit->SetEdepCrystal(oldHit->GetEdepCrystal() +edep);
00220 break;
00221 }
00222 }
00223 }
00224 }
00225
00226 }
00227
00228
00229 if(CalTruthList)
00230 {
00231 G4int trackIndex, g4TrackId;
00232 GetCurrentTrackIndex(trackIndex, g4TrackId);
00233 Identifier id = EmcID::crystal_id(PartId,CryNumberTheta,CryNumberPhi);
00234
00235 G4int flag=1;
00236 if(CalTruthList->entries()>0) {
00237 for(G4int i=0;i<CalTruthList->entries();i++) {
00238 BesEmcTruthHit* oldHit=(*CalTruthList)[i];
00239 if(oldHit->GetTrackIndex()==trackIndex) {
00240 flag=0;
00241 oldHit->SetEDep(oldHit->GetEDep()+edep);
00242 if(oldHit->Find(id)!=oldHit->End()) {
00243 oldHit->AddEHit(id,edep);
00244 } else {
00245 oldHit->Insert(id,edep);
00246 }
00247 break;
00248 }
00249 }
00250 }
00251
00252 if(flag==1) {
00253 G4int pdg = abs(aStep->GetTrack()->GetDefinition()->GetPDGEncoding());
00254 if(!(pdg==12 || pdg==14 || pdg==16)) {
00255 BesEmcTruthHit* truHit=new BesEmcTruthHit;
00256 truHit->SetTrackIndex(trackIndex);
00257 truHit->SetG4TrackId(g4TrackId);
00258 truHit->SetEDep(edep);
00259 truHit->SetIdentify(id);
00260 truHit->Insert(id,edep);
00261
00262 if(aStep->GetTrack()->GetTrackID()==g4TrackId) {
00263 truHit->SetHitEmc(1);
00264 truHit->SetPDGCode(aStep->GetTrack()->GetDefinition()->GetPDGEncoding());
00265 truHit->SetPDGCharge(aStep->GetTrack()->GetDefinition()->GetPDGCharge());
00266 truHit->SetParticleName(aStep->GetTrack()->GetDefinition()->GetParticleName());
00267 truHit->SetTime(aStep->GetPreStepPoint()->GetGlobalTime());
00268 truHit->SetPosition(aStep->GetPreStepPoint()->GetPosition());
00269 truHit->SetMomentum(aStep->GetPreStepPoint()->GetMomentum());
00270
00271 } else {
00272
00273 BesSensitiveManager *sensitiveManager = BesSensitiveManager::GetSensitiveManager();
00274 std::vector<BesTruthTrack*> *trackList = sensitiveManager->GetTrackList();
00275
00276 for(unsigned i=0;i<trackList->size();i++) {
00277 BesTruthTrack *truthTrack=(*trackList)[i];
00278
00279 if(trackIndex==truthTrack->GetIndex()) {
00280 truHit->SetHitEmc(0);
00281 truHit->SetPDGCode(truthTrack->GetPDGCode());
00282 truHit->SetPDGCharge(truthTrack->GetPDGCharge());
00283 truHit->SetParticleName(truthTrack->GetParticleName());
00284
00285 if(truthTrack->GetTerminalVertex()) {
00286 truHit->SetTime(truthTrack->GetTerminalVertex()->GetTime());
00287 truHit->SetPosition(truthTrack->GetTerminalVertex()->GetPosition());
00288 } else {
00289 truHit->SetTime(-99);
00290 truHit->SetPosition(G4ThreeVector(-99,-99,-99));
00291 }
00292
00293 break;
00294 }
00295
00296 }
00297 }
00298
00299 CalTruthList->insert(truHit);
00300 }
00301 }
00302 }
00303
00304 return true;
00305 }
00306
00307
00308
00309 void BesEmcSD::ComputeThetaPhi(G4int id, G4int sector, G4int nb)
00310 {
00311 if((sector>=0)&&(sector<3))
00312 sector+=16;
00313
00314 {
00315 if((nb>=0)&&(nb<4))
00316 {
00317 CryNumberTheta = 0;
00318 CryNumberPhi = (sector-3)*4+nb;
00319 }
00320 else if((nb>=4)&&(nb<8))
00321 {
00322 CryNumberTheta = 1;
00323 CryNumberPhi = (sector-3)*4+nb-4;
00324 }
00325 else if((nb>=8)&&(nb<13))
00326 {
00327 CryNumberTheta = 2;
00328 CryNumberPhi = (sector-3)*5+nb-8;
00329 }
00330 else if((nb>=13)&&(nb<18))
00331 {
00332 CryNumberTheta = 3;
00333 CryNumberPhi = (sector-3)*5+nb-13;
00334 }
00335 else if((nb>=18)&&(nb<24))
00336 {
00337 CryNumberTheta = 4;
00338 CryNumberPhi = (sector-3)*6+nb-18;
00339 }
00340 else if((nb>=24)&&(nb<30))
00341 {
00342 CryNumberTheta = 5;
00343 CryNumberPhi = (sector-3)*6+nb-24;
00344 }
00345 }
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380 if(id==2)
00381 {
00382 switch(CryNumberTheta)
00383 {
00384 case 0:
00385 if(CryNumberPhi<32)
00386 CryNumberPhi = 31-CryNumberPhi;
00387 else
00388 CryNumberPhi = 95-CryNumberPhi;
00389 break;
00390 case 1:
00391 if(CryNumberPhi<32)
00392 CryNumberPhi = 31-CryNumberPhi;
00393 else
00394 CryNumberPhi = 95-CryNumberPhi;
00395 break;
00396 case 2:
00397 if(CryNumberPhi<40)
00398 CryNumberPhi = 39-CryNumberPhi;
00399 else
00400 CryNumberPhi = 119-CryNumberPhi;
00401 break;
00402 case 3:
00403 if(CryNumberPhi<40)
00404 CryNumberPhi = 39-CryNumberPhi;
00405 else
00406 CryNumberPhi = 119-CryNumberPhi;
00407 break;
00408 case 4:
00409 if(CryNumberPhi<48)
00410 CryNumberPhi = 47-CryNumberPhi;
00411 else
00412 CryNumberPhi = 143-CryNumberPhi;
00413 break;
00414 case 5:
00415 if(CryNumberPhi<48)
00416 CryNumberPhi = 47-CryNumberPhi;
00417 else
00418 CryNumberPhi = 143-CryNumberPhi;
00419 default:
00420 ;
00421 }
00422 }
00423 }
00424
00425 G4int BesEmcSD::EndPhiNumberForGDML(G4int copyNb)
00426 {
00427 if(copyNb<0||copyNb>15) {
00428 G4Exception("Wrong copy number of EMC Endcap Phi Volume!");
00429 }
00430
00431 if(copyNb==0||copyNb==1) {
00432 return 15-copyNb*8;
00433 } else if(copyNb==2||copyNb==3) {
00434 return 30-copyNb*8;
00435 } else if(copyNb<=9) {
00436 return 17-copyNb;
00437 } else {
00438 return 15-copyNb;
00439 }
00440 }
00441
00442 G4int BesEmcSD::ComputeEndCopyNb(G4int num)
00443 {
00444 G4int copyNb;
00445 switch(num){
00446 case 30:
00447 copyNb = 5;
00448 break;
00449 case 31:
00450 copyNb = 6;
00451 break;
00452 case 32:
00453 copyNb = 14;
00454 break;
00455 case 33:
00456 copyNb = 15;
00457 break;
00458 case 34:
00459 copyNb = 16;
00460 break;
00461 default:
00462 copyNb = num;
00463 break;
00464 }
00465 return copyNb;
00466 }
00467
00468
00469
00470 void BesEmcSD::EndOfEvent(G4HCofThisEvent* HCE)
00471 {
00472 static G4int HCID = -1;
00473 if(HCID<0)
00474 { HCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]); }
00475 HCE->AddHitsCollection(HCID,CalCollection);
00476 }
00477
00478
00479
00480 void BesEmcSD::clear()
00481 {}
00482
00483
00484
00485 void BesEmcSD::DrawAll()
00486 {}
00487
00488
00489
00490 void BesEmcSD::PrintAll()
00491 {}
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511