Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 26 //GEANT4 - Depth-of-Interaction enabled Positr 27 28 //Authors and contributors 29 30 // Author list to be updated, with names of co 31 32 // Abdella M. Ahmed (1, 2), Andrew Chacon (1, 33 // Hideaki Tashima (3), Go Akamatsu (3), Akram 34 // Susanna Guatelli (2), and Mitra Safavi-Naei 35 36 // (1) Australian Nuclear Science and Technolo 37 // (2) University of Wollongong, Australia 38 // (3) National Institute of Radiological Scie 39 40 //Implemetation of the doiPETAnalysis.cc class 41 //This implementation mimics readout (or digit 42 //parameters are given in inputParameters.txt 43 //each detector block and axially multiplexed 44 //blurring parameters are in keV (for energy) 45 //program quits. In this class, an ideal PMT i 46 //(obtained by G4) is used to determine the di 47 //distance from the PMTs. Light sharing method 48 //error message will be displayed and the even 49 //processed into coinsidence list-mode data. A 50 //Explanation is given for the methods provide 51 52 53 #include "doiPETAnalysis.hh" 54 #include "G4SystemOfUnits.hh" 55 #include "G4PhysicalConstants.hh" 56 #include <iomanip> 57 #include "Randomize.hh" 58 #include "G4SPSRandomGenerator.hh" 59 #include "doiPETAnalysisMessenger.hh" 60 #include "G4UnitsTable.hh" 61 #include "globals.hh" 62 63 doiPETAnalysis* doiPETAnalysis::instance=0; 64 65 /////////// Constructor ////////////////////// 66 doiPETAnalysis::doiPETAnalysis() 67 { 68 factoryOn = false; //G4ROOT 69 // Initialization ntuple 70 for (G4int k=0; k<MaxNtCol; k++) fNtColId[ 71 72 fAnalysisMessenger = new doiPETAnalysisMesse 73 74 //Set energy window 75 lowerThreshold = 400*keV; 76 upperThreshold = 600*keV; 77 triggerEnergy = 50*eV; 78 79 //give default initial activity. Activity st 80 InitialActivity = 1000000*becquerel; 81 82 //In NEMA NU2, all test is done with F-18 83 halfLife = 109.771*60 * s;//The halfLife of 84 85 // 86 totalTime = 0 * s; 87 prev_totalTime = 0 * s; 88 prev_eventID = 0; 89 90 // 91 //Initialize crystal ID 92 crystalIDNew_tan = -1; 93 crystalIDNew_axial = -1; 94 crystalIDNew_DOI = -1; 95 96 // 97 scatterIndex = 0; 98 99 // 100 numberOfPixel_tan = 2*numberOfCrystal_tangen 101 numberOfPixel_axial = 2*numberOfCrystal_axia 102 103 //Default value for deadtime. 104 block_DeadTime = 256*ns; 105 module_DeadTime = 0*ns; 106 // 107 108 //Crystal blurring parameters. One detector 109 //So, a range of energy resolution is applie 110 //The energy resolution can be set in the in 111 crystalResolutionMin = 0.13;//13% 112 crystalResolutionMax = 0.17;//17% 113 114 crystalEnergyRef = 511 * keV;//Energy of ref 115 116 //The quantum efficiency models the probabil 117 //The quantum efficiency can be set in the i 118 crystalQuantumEfficiency = 1;//100% 119 // 120 121 //intialize deadtime for blocks and modules 122 numberOfBlocks_total = numberOfRings * numbe 123 blockTime = new double[numberOfBlocks_total] 124 moduleTime = new double[numberOfBlocks_total 125 126 //Initialize the deadtime for each detector 127 for(G4int i = 0; i<numberOfBlocks_total; i++ 128 blockTime [i] = 0.0; 129 moduleTime [i] = 0.0; 130 } 131 132 //Initialize type of output. The default out 133 getSinglesData = false; //default value 134 getCoincidenceData = false; 135 136 ApplyAngerLogic = true; 137 isDOIlookUpTablePrepared = false; 138 139 numberOfHit = 0; 140 141 //This value is based on the assumption that 142 shiftCoeff = 0.5; 143 144 // 145 PMTblurring_tan = 0.0; 146 PMTblurring_axial = 0.0; 147 } 148 ////////// Destructor //////////////////////// 149 doiPETAnalysis::~doiPETAnalysis() 150 { 151 delete fAnalysisMessenger; 152 delete [] blockTime; 153 delete [] moduleTime; 154 } 155 156 ////////// GetInstance /////////////////////// 157 doiPETAnalysis* doiPETAnalysis::GetInstance() 158 { 159 if(instance==0) instance = new doiPETAnalysi 160 return instance; 161 } 162 void doiPETAnalysis::Delete() 163 { 164 delete instance; 165 } 166 167 //If there is energy deposition in the phantom 168 //Use this for checking 169 void doiPETAnalysis::SetScatterIndexInPhantom( 170 scatterIndex = sci; 171 } 172 173 //Get the source position if the process is an 174 //Use this for checking 175 void doiPETAnalysis::SetSourcePosition(G4Three 176 spositionX = spos.x(); 177 spositionY = spos.y(); 178 spositionZ = spos.z(); 179 } 180 181 182 //Set the event ID 183 //Use this for checking. eventID should not be 184 void doiPETAnalysis::SetEventID(G4int evID){ 185 eventID = evID; 186 } 187 188 // 189 190 void doiPETAnalysis::GetSizeOfDetector(G4doubl 191 sizeOfDetector_DOI = detSizeDoi; 192 sizeOfDetector_axial = detSizeTan; 193 sizeOfDetector_tangential = detSizeAxial; 194 } 195 196 // 197 void doiPETAnalysis::SetActivity(G4double newA 198 InitialActivity = newActivity; 199 G4cout<<"Initial activity: "<<InitialActivit 200 } 201 void doiPETAnalysis::SetIsotopeHalfLife(G4doub 202 halfLife = newHalfLife; 203 G4cout<<"Half life of the isotope "<<halfLif 204 } 205 206 //The time is based on random time intervals b 207 //This time mimics acquisition time of a PET s 208 void doiPETAnalysis::CalulateAcquisitionTime() 209 //Calculate the strength of activity at t=to 210 activityNow = InitialActivity * std::exp(-(( 211 212 //Activity based time interval. 213 timeInterval = -std::log(G4UniformRand())*(1 214 totalTime = timeInterval+prev_totalTime; 215 prev_totalTime = totalTime; 216 } 217 218 //Apply energy blurring on the crystals. The v 219 G4double doiPETAnalysis::QuantumEffifciency(G4 220 { 221 if(fixedResolution){ 222 crystalResolution = energyResolution_fixed 223 } 224 else{ 225 crystalResolution = energyResolution_cryDe 226 } 227 crystalCoeff = crystalResolution * std::sqrt 228 229 G4double QE = G4UniformRand(); 230 231 //The quantum efficiency models the probabil 232 if(QE <= crystalQuantumEfficiency) 233 { 234 edep_AfterCrystalBlurring = G4RandGauss::s 235 } 236 else { 237 //not detected by the photodetector, event 238 edep_AfterCrystalBlurring = 0 *keV; 239 } 240 return edep_AfterCrystalBlurring; 241 } 242 243 ///////// ReadOut //////////////////////////// 244 245 void doiPETAnalysis::ReadOut(G4int blkID, G4in 246 { 247 blockID = blkID; 248 crystalID = cryID; 249 interactionTime = interTime; 250 time_annihil = timeAnnih; 251 interactionPos = interPos; 252 totalEdep = edep; 253 254 //Get the time of flight. This is the durati 255 time_tof = interactionTime - time_annihil; 256 257 //time of the event when detected (timerTag) 258 timeStamp = totalTime + time_tof; 259 260 //triggerEnergy is the energy deposited in t 261 if(totalEdep<triggerEnergy)return; 262 263 //************************************** App 264 //Apply paralizable dead-time in the block b 265 if(std::fabs(timeStamp - blockTime[blockID]) 266 blockTime[blockID] = timeStamp; 267 } 268 else { 269 //If the time difference is less than the 270 blockTime[blockID] = timeStamp; 271 272 //the event is then rejected 273 //continue; 274 return; 275 } 276 277 //Apply Non-paralyzable dead-time on axially 278 //If the time difference is less than the pr 279 if(std::fabs(timeStamp - moduleTime[blockID] 280 281 //The following finds the block id's of fo 282 for (G4int r_ring = 0; r_ring < numberOfRi 283 if (blockID >= r_ring*numberOfDetector_p 284 for (G4int m_module = 0; m_module < nu 285 286 //Set the time of the module (four b 287 moduleTime[blockID + (m_module - r_r 288 } 289 } 290 } 291 } 292 else return; 293 294 295 ///////////////////////// Write qualified 296 297 if(totalEdep>lowerThreshold && totalEdep<upp 298 299 //identifiy the layer 300 DOI_ID = G4int(crystalID/(numberOfCrystal 301 302 //identify the crystal id for each Layer. 303 crystalID_2D = crystalID - (DOI_ID*numberO 304 305 //identify the crystal ID in the tangentia 306 crystalID_axial = crystalID_2D/numberOfCry 307 crystalID_tangential = crystalID_2D%number 308 309 intPosX = interactionPos.x(); 310 intPosY = interactionPos.y(); 311 intPosZ = interactionPos.z(); 312 313 314 if(ApplyAngerLogic){ 315 //shiftCoeff = 0.5 is used. This value i 316 AngerLogic(intPosX, intPosY, intPosZ, to 317 } 318 319 //Single event output. Coincidence events 320 if(getSinglesData) WriteOutput(); 321 322 //Coincidence output 323 if(getCoincidenceData){ 324 eventID_coin.push_back(eventID); 325 blockID_coin.push_back(blockID); 326 cryID_axial_coin.push_back(crystalID_axi 327 cryID_tan_coin.push_back(crystalID_tange 328 edep_coin.push_back(totalEdep); 329 cryDOI_coin.push_back(DOI_ID); 330 time_coin.push_back(timeStamp); 331 332 numberOfHit++; 333 334 if(numberOfHit == 2){ //two events withi 335 WriteOutput(); 336 ResetNumberOfHits(); 337 } 338 } 339 } 340 } 341 342 343 344 ////////// Clear ///////////////////////////// 345 void doiPETAnalysis::ResetNumberOfHits() 346 { 347 numberOfHit = 0; 348 eventID_coin.clear(); 349 blockID_coin.clear(); 350 cryID_axial_coin.clear(); 351 cryID_tan_coin.clear(); 352 edep_coin.clear(); 353 cryDOI_coin.clear(); 354 time_coin.clear(); 355 356 } 357 358 // 359 void doiPETAnalysis::Open(G4String fileName) 360 { 361 if(getSinglesData){ 362 asciiFileName = fileName + "Singles.data"; 363 rootFileName = fileName+"Singles.root"; 364 } 365 if(getCoincidenceData){ 366 asciiFileName = fileName + "Coincidence.da 367 rootFileName = fileName+"Coincidence.root" 368 } 369 370 ofs.open(asciiFileName.c_str()); 371 if(!ofs.is_open()){ 372 G4cerr<<"=== \n File opening Error to writ 373 exit(0); 374 } 375 376 } 377 378 void doiPETAnalysis::WriteOutput(){ 379 if(getSinglesData){ 380 ofs<<eventID<<" "<<blockID<<" "<<std::setp 381 //ofs<<eventID<<" "<<blockID<<" "<<crystal 382 383 } 384 if(getCoincidenceData){ 385 //2 singles will qualify to be in coincide 386 for(G4int i=0; i<2; i++){ 387 388 //First Single 389 if(i==0){ 390 eventID0 = eventID_coin[0]; 391 blockID0 = blockID_coin[0]; 392 crystalID_axial0 = cryID_axial_coin 393 crystalID_tangential0 = cryID_tan_coin 394 DOI_ID0 = cryDOI_coin[0]; 395 timeStamp0 = time_coin[0]; 396 totalEdep0 = edep_coin[0]; 397 } 398 if(i==1){ 399 //Second Single 400 eventID1 = eventID_coin[1]; 401 blockID1 = blockID_coin[1]; 402 crystalID_axial1 = cryID_axial_coin 403 crystalID_tangential1 = cryID_tan_coin 404 DOI_ID1 = cryDOI_coin[1]; 405 timeStamp1 = time_coin[1]; 406 totalEdep1 = edep_coin[1]; 407 } 408 } 409 410 ofs<<eventID0<<" "<<blockID0<<" "<<crystal 411 <<eventID1<<" "<<blockID1<<" "<<crystalI 412 <<spositionX<<" "<<spositionY<<" "<<spos 413 414 } 415 #ifdef ANALYSIS_USE 416 FillListModeEvent(); 417 #endif 418 419 } 420 421 // 422 ///////// Close ////////////////////////////// 423 void doiPETAnalysis::Close() 424 { 425 //close ascii file 426 ofs.close(); 427 428 } 429 430 //Place the photomultiplier tube (PMT) at each 431 //All the PMTs are placed at the same doi (x) 432 433 //The PMT is placed at each corner of the crys 434 //The signal (energy deposition) of each PMT d 435 void doiPETAnalysis::PMTPosition(){ 436 437 sizeOfDetector_DOI = (numberOfCrystal_DOI * 438 sizeOfDetector_axial = (numberOfCrystal_axia 439 sizeOfDetector_tangential = (numberOfCrystal 440 441 //Position of PMT1. 442 posPMT1x = sizeOfDetector_DOI/2;//mm 443 posPMT1y = -sizeOfDetector_tangential/2; 444 posPMT1z = -sizeOfDetector_axial/2; 445 446 //Position of PMT2 447 posPMT2x = sizeOfDetector_DOI/2; 448 posPMT2y = sizeOfDetector_tangential/2; 449 posPMT2z = -sizeOfDetector_axial/2; 450 451 //Position of PMT3 452 posPMT3x = sizeOfDetector_DOI/2; 453 posPMT3y = -sizeOfDetector_tangential/2; 454 posPMT3z = sizeOfDetector_axial/2; 455 456 //Position of PMT4 457 posPMT4x = sizeOfDetector_DOI/2; 458 posPMT4y = sizeOfDetector_tangential/2; 459 posPMT4z = sizeOfDetector_axial/2; 460 461 G4cout<<"PMT positions: "<<G4endl; 462 G4cout<<"PMT1 (mm) ("<<posPMT1x<<", "<<posPM 463 G4cout<<"PMT2 (mm) ("<<posPMT2x<<", "<<posPM 464 G4cout<<"PMT3 (mm) ("<<posPMT3x<<", "<<posPM 465 G4cout<<"PMT4 (mm) ("<<posPMT4x<<", "<<posPM 466 467 } 468 469 //The blurring parameters are given and can be 470 void doiPETAnalysis::BlurringParameters(){ 471 char inputChar[256]; 472 std::string inputLine; 473 G4String value[7]; 474 std::string filename = "inputParameter.txt" 475 ifs.open(filename.c_str()); 476 if(!ifs.good()){ 477 G4cerr<<"File opening Error: Could not ope 478 exit(0); 479 } 480 while(!ifs.eof()){ 481 ifs.getline(inputChar,256); 482 inputLine = inputChar; 483 if(inputChar[0]!='#' && inputLine.length() 484 if( (std::string::size_type)inputLine.fi 485 std::istringstream tmpStream(inputLine 486 tmpStream >> value[0] >> value[1] >> v 487 block_DeadTime = atof(value[1].c_str() 488 if(value[2] != "ns"){ 489 G4cerr<<" Dead time unit is not in n 490 exit(0); 491 } 492 block_DeadTime = block_DeadTime*ns; 493 G4cout<<"Dead time of the detector: "< 494 } 495 if( (std::string::size_type)inputLine.fi 496 std::istringstream tmpStream(inputLine 497 tmpStream >> value[0] >> value[1] >> v 498 module_DeadTime = atof(value[1].c_str( 499 if(value[2] != "ns"){ 500 G4cerr<<" Dead time unit is not in n 501 exit(0); 502 } 503 module_DeadTime = module_DeadTime*ns; 504 G4cout<<"Dead time of the module (axia 505 } 506 // 507 if( (std::string::size_type)inputLine.fi 508 std::istringstream tmpStream(inputLine 509 tmpStream >> value[0] >> value[1]; 510 crystalResolutionMin = atof(value[1].c 511 G4cout<<"crystal Resolution (Min.): "< 512 } 513 if( (std::string::size_type)inputLine.fi 514 std::istringstream tmpStream(inputLine 515 tmpStream >> value[0] >> value[1]; 516 crystalResolutionMax = atof(value[1].c 517 G4cout<<"crystal Resolution (Max.): "< 518 } 519 520 // 521 if( (std::string::size_type)inputLine.fi 522 std::istringstream tmpStream(inputLine 523 tmpStream >> value[0] >> value[1]; 524 if(value[1]=="true"){ 525 fixedResolution = true; 526 energyResolution_fixed = (crystalRes 527 G4cout<<"Fixed crystal resolution is 528 } 529 else { 530 fixedResolution = false; 531 //Store into a file if needed. 532 std::string fname = "crystalDependen 533 std::ofstream outFname(fname.c_str() 534 535 G4cout<<" \n Crystal dependent resol 536 energyResolution_cryDependent.resize 537 for(G4int i_blk = 0; i_blk < numberO 538 for(G4int i_cry = 0; i_cry < numbe 539 energyResolution_cryDependent[i_ 540 //store into a file 541 outFname<<i_blk<<" "<<i_cry<<" " 542 } 543 } 544 G4cout<<"Done. \n"<<G4endl; 545 outFname.close(); 546 } 547 548 } 549 // 550 551 if( (std::string::size_type)inputLine.fi 552 std::istringstream tmpStream(inputLine 553 tmpStream >> value[0] >> value[1] >> v 554 crystalEnergyRef = atof(value[1].c_str 555 if(value[2] != "keV"){ 556 G4cerr<<" The unit of reference ener 557 exit(0); 558 } 559 crystalEnergyRef = crystalEnergyRef*ke 560 G4cout<<"Energy of refernce: "<<crysta 561 } 562 if( (std::string::size_type)inputLine.fi 563 std::istringstream tmpStream(inputLine 564 tmpStream >> value[0] >> value[1]; 565 crystalQuantumEfficiency = atof(value[ 566 G4cout<<"Quantum Efficiency "<<crystal 567 } 568 if( (std::string::size_type)inputLine.fi 569 std::istringstream tmpStream(inputLine 570 tmpStream >> value[0] >> value[1] >> v 571 lowerThreshold = atof(value[1].c_str() 572 if(value[2] != "keV"){ 573 G4cerr<<" The unit of Lower energy t 574 exit(0); 575 } 576 lowerThreshold = lowerThreshold*keV; 577 G4cout<<"Lower energy threshold: "<<lo 578 579 } 580 if( (std::string::size_type)inputLine.fi 581 std::istringstream tmpStream(inputLine 582 tmpStream >> value[0] >> value[1] >> v 583 upperThreshold = atof(value[1].c_str() 584 if(value[2] != "keV"){ 585 G4cerr<<" The unit of Upper energy t 586 exit(0); 587 } 588 upperThreshold = upperThreshold*keV; 589 G4cout<<"Upper energy threshold: "<<up 590 591 } 592 593 // 594 if( (std::string::size_type)inputLine.fi 595 std::istringstream tmpStream(inputLine 596 tmpStream >> value[0] >> value[1] >> v 597 triggerEnergy = atof(value[1].c_str()) 598 if(value[2] != "keV"){ 599 G4cerr<<" The unit of Trigger energy 600 exit(0); 601 } 602 triggerEnergy = triggerEnergy*keV; 603 G4cout<<"Trigger energy threshold: "<< 604 605 } 606 607 //Option to apply AngerLogic 608 if( (std::string::size_type)inputLine.fi 609 std::istringstream tmpStream(inputLine 610 tmpStream >> value[0] >> value[1]; 611 if(value[1]=="true"){ 612 ApplyAngerLogic = true; 613 G4cout<<"Angler Logic calculation is 614 } 615 else if(value[1]=="false") { 616 ApplyAngerLogic = false; 617 G4cout<<"Angler Logic calculation is 618 } 619 else { 620 ApplyAngerLogic = true; 621 G4cout<<"Angler Logic calculation is 622 } 623 } 624 625 //PMT position calculation blurring at F 626 if( (std::string::size_type)inputLine.fi 627 std::istringstream tmpStream(inputLine 628 tmpStream >> value[0] >> value[1]>>val 629 PMTblurring_tan = atof(value[1].c_str( 630 PMTblurring_axial = atof(value[2].c_st 631 G4cout<<"PMTblurring position response 632 } 633 634 // 635 if( (std::string::size_type)inputLine.fi 636 std::istringstream tmpStream(inputLine 637 tmpStream >> value[0] >> value[1]; 638 if(value[1]=="singlesOutput"){ 639 getSinglesData = true; 640 G4cout<<"Single mode output enabled. 641 } 642 else if(value[1]=="coincidenceOutput") 643 getCoincidenceData = true; 644 G4cout<<"Coicidence mode output enab 645 } 646 647 } 648 649 } 650 } 651 ifs.close(); 652 } 653 654 //The following function reads the reflector p 655 //For defualt reflector pattern, see https://l 656 //The patter of the reflectors can be changed 657 //The pattern is given as 0 and 1. If there is 658 659 void doiPETAnalysis::ReadReflectorPattern(){ 660 G4cout<<" Reflector pattern is being read "< 661 // 662 std::vector<std::string> stringReflectorValu 663 // 664 char inputChar[256]; 665 std::string inputLine; 666 667 //open inputParameter.txt to read reflector 668 std::string filename = "inputParameter.txt" 669 670 G4String refValue; 671 672 ifs.open(filename.c_str()); 673 if(!ifs.good()){ 674 G4cerr<<"File opening Error: Could not ope 675 exit(0); 676 } 677 while(!ifs.eof()){ 678 ifs.getline(inputChar,256); 679 inputLine = inputChar; 680 681 //The reflector patter in read from the in 682 if(inputChar[0]!='#' && inputLine.length() 683 684 //Reflector patter for Layer1 in the tan 685 if( (std::string::size_type)inputLine.fi 686 std::istringstream tmpStream(inputLine 687 while(tmpStream >> refValue){ 688 stringReflectorValue.push_back(refVa 689 if(stringReflectorValue.size()>1){ 690 G4int tmp_value = atoi(stringRefle 691 ireflectorLayer1_Tangential.push_b 692 } 693 } 694 stringReflectorValue.clear(); 695 } 696 697 698 //Reflector patter for Layer1 in the axi 699 if( (std::string::size_type)inputLine.fi 700 std::istringstream tmpStream(inputLine 701 while(tmpStream >> refValue){ 702 stringReflectorValue.push_back(refVa 703 if(stringReflectorValue.size()>1){ 704 G4int tmp_value = atoi(stringRefle 705 ireflectorLayer1_Axial.push_back(t 706 } 707 } 708 stringReflectorValue.clear(); 709 } 710 711 //Reflector patter for Layer2 in the tan 712 if( (std::string::size_type)inputLine.fi 713 std::istringstream tmpStream(inputLine 714 while(tmpStream >> refValue){ 715 stringReflectorValue.push_back(refVa 716 if(stringReflectorValue.size()>1){ 717 G4int tmp_value = atoi(stringRefle 718 ireflectorLayer2_Tangential.push_b 719 } 720 } 721 stringReflectorValue.clear(); 722 } 723 724 //Reflector patter for Layer2 in the axi 725 if( (std::string::size_type)inputLine.fi 726 std::istringstream tmpStream(inputLine 727 while(tmpStream >> refValue){ 728 stringReflectorValue.push_back(refVa 729 if(stringReflectorValue.size()>1){ 730 G4int tmp_value = atoi(stringRefle 731 ireflectorLayer2_Axial.push_back(t 732 } 733 } 734 stringReflectorValue.clear(); 735 } 736 737 //Reflector patter for Layer3 in the tan 738 if( (std::string::size_type)inputLine.fi 739 std::istringstream tmpStream(inputLine 740 while(tmpStream >> refValue){ 741 stringReflectorValue.push_back(refVa 742 if(stringReflectorValue.size()>1){ 743 G4int tmp_value = atoi(stringRefle 744 ireflectorLayer3_Tangential.push_b 745 } 746 } 747 stringReflectorValue.clear(); 748 } 749 750 //Reflector patter for Layer3 in the axi 751 if( (std::string::size_type)inputLine.fi 752 std::istringstream tmpStream(inputLine 753 while(tmpStream >> refValue){ 754 stringReflectorValue.push_back(refVa 755 if(stringReflectorValue.size()>1){ 756 G4int tmp_value = atoi(stringRefle 757 ireflectorLayer3_Axial.push_back(t 758 } 759 } 760 stringReflectorValue.clear(); 761 } 762 763 //Reflector patter for Layer4 in the tan 764 if( (std::string::size_type)inputLine.fi 765 std::istringstream tmpStream(inputLine 766 while(tmpStream >> refValue){ 767 stringReflectorValue.push_back(refVa 768 if(stringReflectorValue.size()>1){ 769 G4int tmp_value = atoi(stringRefle 770 ireflectorLayer4_Tangential.push_b 771 } 772 } 773 stringReflectorValue.clear(); 774 } 775 776 //Reflector patter for Layer4 in the axi 777 if( (std::string::size_type)inputLine.fi 778 std::istringstream tmpStream(inputLine 779 while(tmpStream >> refValue){ 780 stringReflectorValue.push_back(refVa 781 if(stringReflectorValue.size()>1){ 782 G4int tmp_value = atoi(stringRefle 783 ireflectorLayer4_Axial.push_back(t 784 } 785 } 786 stringReflectorValue.clear(); 787 } 788 }//# 789 }//while(eof) 790 791 // 792 //for debug 793 G4cout<<"\n========= Reflector Pattern ===== 794 G4cout<<"Layer 1"<<G4endl; 795 for(unsigned int i = 0; i<ireflectorLayer1_T 796 G4cout<<ireflectorLayer1_Tangential[i]<<" 797 }G4cout<<G4endl; 798 for(unsigned int i = 0; i<ireflectorLayer1_A 799 G4cout<<ireflectorLayer1_Axial[i]<<" "; 800 }G4cout<<G4endl; 801 G4cout<<"Layer 2"<<G4endl; 802 for(unsigned int i = 0; i<ireflectorLayer2_T 803 G4cout<<ireflectorLayer2_Tangential[i]<<" 804 }G4cout<<G4endl; 805 for(unsigned int i = 0; i<ireflectorLayer2_A 806 G4cout<<ireflectorLayer2_Axial[i]<<" "; 807 }G4cout<<G4endl; 808 G4cout<<"Layer 3"<<G4endl; 809 for(unsigned int i = 0; i<ireflectorLayer3_T 810 G4cout<<ireflectorLayer3_Tangential[i]<<" 811 }G4cout<<G4endl; 812 for(unsigned int i = 0; i<ireflectorLayer3_A 813 G4cout<<ireflectorLayer3_Axial[i]<<" "; 814 }G4cout<<G4endl; 815 G4cout<<"Layer 4"<<G4endl; 816 for(unsigned int i = 0; i<ireflectorLayer4_T 817 G4cout<<ireflectorLayer4_Tangential[i]<<" 818 }G4cout<<G4endl; 819 for(unsigned int i = 0; i<ireflectorLayer4_A 820 G4cout<<ireflectorLayer4_Axial[i]<<" "; 821 }G4cout<<G4endl; 822 G4cout<<"========= Reflector Pattern Ended = 823 824 //Read DOI look-up-table. This look-up-table 825 G4int index_doi = 0, doiID; 826 doi_table.resize(numberOfCrystal_tangential* 827 828 std::string LUT_FileName = "look_up_table_DO 829 std::ifstream ifs_doiLUT; 830 std::ofstream ofs_doiLUT; 831 ifs_doiLUT.open(LUT_FileName.c_str()); 832 if(ifs_doiLUT.is_open()){ 833 G4cout<<" DOI Look-up table found and used 834 //Read from file 835 while(ifs_doiLUT>>index_doi>>doiID && inde 836 doi_table[index_doi] = doiID; 837 } 838 if(index_doi==int(doi_table.size())){ 839 G4cout<<"!!!Warning: The DOI table index 840 PrepareDOILookUpTable(LUT_FileName); 841 } 842 isDOIlookUpTablePrepared = true; // 843 } 844 else 845 { 846 PrepareDOILookUpTable(LUT_FileName); 847 } 848 //Write into a file. 849 ofs_doiLUT.open(LUT_FileName.c_str()); 850 if(!ofs_doiLUT.is_open()){ 851 G4cerr<<"Unable to open file to write doi_ 852 exit(0); 853 } 854 for(G4int i=0;i<int(doi_table.size()); i++){ 855 ofs_doiLUT<<i<<"\t"<<doi_table[i]<<G4endl; 856 } 857 ifs_doiLUT.close(); 858 859 860 ifs.close(); 861 } 862 void doiPETAnalysis::PrepareDOILookUpTable(G4S 863 isDOIlookUpTablePrepared = false; 864 G4cout<<"Preparing DOI look-up table... "< 865 std::string outputFileName = "_check_2Dpos 866 std::ofstream outFile(outputFileName.c_str 867 868 G4double crystalPositionX; 869 G4double crystalPositionY; 870 G4double crystalPositionZ; 871 //doi_table.resize(numberOfCrystal_tangent 872 873 for(G4int i_DOI = 0; i_DOI<numberOfCrystal 874 crystalPositionX=(i_DOI-((float)numberOf 875 for(G4int i_axial=0; i_axial< numberOfCr 876 crystalPositionZ = (i_axial-((float)nu 877 for(G4int i_tan=0; i_tan<numberOfCryst 878 crystalPositionY=(i_tan-((float)numb 879 AngerLogic(crystalPositionX, crystal 880 outFile<<PositionAngerZ<<" "<<Positi 881 doi_table[crystalID_in2D_posHist]=i_ 882 883 } 884 } 885 } 886 G4cout<<"done."<<G4endl; 887 isDOIlookUpTablePrepared = true; 888 } 889 890 891 //Based on ideal photomultiplier tube (PMT) pl 892 //The reflectors shifts the response by some d 893 //From this 2D position histogram, the new cry 894 //If the crystal ID after Anger method apllied 895 void doiPETAnalysis::AngerLogic(G4double posCr 896 { 897 G4double crystalPitch_DOI = sizeOfCrystal_DO 898 G4double crystalPitch_tan = sizeOfCrystal_ta 899 G4double crystalPitch_axial = sizeOfCrystal_ 900 901 //The crystal ID are calculated based on the 902 G4int i_doi = posCrystalX/crystalPitch_DOI + 903 G4int i_tan = posCrystalY/crystalPitch_tan + 904 G4int i_axial = posCrystalZ/crystalPitch_axi 905 906 //position of interaction is shifted the ce 907 posCrystalX = (i_doi-((float)numberOfCrystal 908 posCrystalY = (i_tan-((float)numberOfCrystal 909 posCrystalZ = (i_axial-((float)numberOfCryst 910 911 //1z and 2z are at the same z distance; 3z a 912 //The signal (the energy deposition) is devi 913 914 //The following is based on symetrical placm 915 916 //Calculate the axial (z) distance from the 917 dist1z = std::fabs(posPMT1z - posCrystalZ); 918 dist2z = std::fabs(posPMT2z - posCrystalZ); 919 dist3z = std::fabs(posPMT3z - posCrystalZ); 920 dist4z = std::fabs(posPMT4z - posCrystalZ); 921 922 //Calculate the resultant distance 923 //dist1z = dist2z, and dist3z = dist4z, so o 924 //distz = ((dist1z + dist2z) + (dist3z + dis 925 distz = dist1z + dist3z; 926 927 //1y and 3y are at the same y distance; and 928 //Calculate the tangential (y) distance fro 929 dist1y = std::fabs(posPMT1y - posCrystalY); 930 dist2y = std::fabs(posPMT2y - posCrystalY); 931 dist3y = std::fabs(posPMT3y - posCrystalY); 932 dist4y = std::fabs(posPMT4y - posCrystalY); 933 934 //Calculate the resultant distance 935 //dist1y = dist3y, and dist2y = dist4y, so o 936 //disty = ((dist1y + dist3y) + (dist2y+dist4 937 disty = dist1y + dist2y; 938 939 //signalPMT1z = signalPMT2z, and signalPMT3z 940 signalPMT1z = Edep * dist3z/(dist1z + dist3z 941 signalPMT3z = Edep * dist1z/(dist1z + dist3z 942 943 signalPMT2z = Edep * dist4z/(dist2z + dist4z 944 signalPMT4z = Edep * dist2z/(dist2z + dist4z 945 946 947 //signalPMT1y = signalPMT3y, and signalPMT2y 948 signalPMT1y = Edep * dist2y/(dist1y + dist2y 949 signalPMT2y = Edep * dist1y/(dist1y + dist2y 950 951 signalPMT3y = Edep * dist4y/(dist3y + dist4y 952 signalPMT4y = Edep * dist3y/(dist3y + dist4y 953 954 //Calculate the signal on each PMT from the 955 signalPMT1 = (signalPMT1z + signalPMT1y)*0. 956 signalPMT2 = (signalPMT2z + signalPMT2y)*0. 957 signalPMT3 = (signalPMT3z + signalPMT3y)*0. 958 signalPMT4 = (signalPMT4z + signalPMT4y)*0. 959 960 961 signalZplus = (signalPMT3 + signalPMT4); 962 signalZminus = (signalPMT1 + signalPMT2); 963 signalYplus = (signalPMT2 + signalPMT4); 964 signalYminus = (signalPMT1 + signalPMT3); 965 966 967 //Position of interaction is calculated base 968 //To get the position by Anger calculation, 969 PositionAngerZ = (signalZplus - signalZminus 970 PositionAngerY = (signalYplus - signalYminus 971 972 973 //For detectors with reflector insertion (li 974 //Here, it is assumed that the shift of the 975 976 //If reflector is only in the left side of t 977 //If reflector is only in the right side of 978 979 //Response shift for 1st Layer 980 if(i_doi == 0){ 981 //If reflector is only in one (left) side 982 if(ireflectorLayer1_Tangential[i_tan] == 1 983 984 //If reflector is only in one (right) side 985 if(ireflectorLayer1_Tangential[i_tan] == 0 986 987 if(ireflectorLayer1_Axial[i_axial] == 1 && 988 if(ireflectorLayer1_Axial[i_axial] == 0 && 989 } 990 if(i_doi == 1){ //Response shift for 2nd Lay 991 if(ireflectorLayer2_Tangential[i_tan] == 1 992 if(ireflectorLayer2_Tangential[i_tan] == 0 993 994 if(ireflectorLayer2_Axial[i_axial] == 1 && 995 if(ireflectorLayer2_Axial[i_axial] == 0 && 996 } 997 if(i_doi == 2){ //Response shift for 3rd Lay 998 if(ireflectorLayer3_Tangential[i_tan] == 1 999 if(ireflectorLayer3_Tangential[i_tan] == 0 1000 1001 if(ireflectorLayer3_Axial[i_axial] == 1 & 1002 if(ireflectorLayer3_Axial[i_axial] == 0 & 1003 } 1004 if(i_doi == 3){ //Response shift for 4th La 1005 if(ireflectorLayer4_Tangential[i_tan] == 1006 if(ireflectorLayer4_Tangential[i_tan] == 1007 1008 if(ireflectorLayer4_Axial[i_axial] == 1 & 1009 if(ireflectorLayer4_Axial[i_axial] == 0 & 1010 } 1011 1012 //Blur the 2D position (obtained by ANger 1013 if(isDOI_LUT){ 1014 PositionAngerZ = G4RandGauss::shoot(Posit 1015 PositionAngerY = G4RandGauss::shoot(Posit 1016 } 1017 //The main purpose of shifting the response 1018 //by comparing with a look-up-table which i 1019 1020 //The crystal ID in 2D position histogram a 1021 crystalID_in2D_posHist_axial = (G4int)(Posi 1022 1023 //The crystal ID in 2D position histogram a 1024 crystalID_in2D_posHist_tan = (G4int)(Posi 1025 1026 //continuous crystal ID in the 2D position 1027 crystalID_in2D_posHist = crystalID_in2D_pos 1028 1029 1030 //Now, lets find the crystal ID in 3D after 1031 1032 //Crystal ID along the tangential diraction 1033 crystalIDNew_tan = (G4int)(crystalID_in2D_p 1034 1035 //Crystal ID along the axial diraction afte 1036 crystalIDNew_axial = (G4int)(crystalID_in2D 1037 1038 ////Crystal ID along the DOI diraction afte 1039 if(crystalID_in2D_posHist>numberOfCrystal_t 1040 crystalIDNew_DOI = doi_table[crystalID_in2D 1041 1042 //If the crsytal ID is beyond the given the 1043 if(crystalIDNew_tan < 0 || crystalIDNew_axi 1044 crystalIDNew_tan >= numberOfCrystal_tange 1045 return; 1046 } 1047 1048 CrystalIDAfterAngerLogic(crystalIDNew_tan,c 1049 } 1050 1051 ///// 1052 void doiPETAnalysis::CrystalIDAfterAngerLogic 1053 crystalID_tangential = i_tan; 1054 crystalID_axial = i_axial; 1055 DOI_ID = i_doi; 1056 } 1057 1058 void doiPETAnalysis::book() 1059 { 1060 auto manager = G4AnalysisManager::Instance( 1061 1062 //manager->SetVerboseLevel(2); 1063 1064 G4bool fileOpen = manager->OpenFile(rootFil 1065 if (!fileOpen) { 1066 G4cout << "\n---> HistoManager::book(): c 1067 << rootFileName << G4endl; 1068 return; 1069 } 1070 // Create directories 1071 //manager->SetNtupleDirectoryName("ListMode 1072 1073 manager->SetFirstNtupleId(1); 1074 1075 if(getSinglesData){ 1076 manager -> CreateNtuple("Singles", "Singles 1077 fNtColId[0] = manager -> CreateNtupleIColum 1078 fNtColId[1] = manager -> CreateNtupleIColum 1079 //fNtColId[2] = manager -> CreateNtupleDCol 1080 //fNtColId[3] = manager -> CreateNtupleDCol 1081 //fNtColId[4] = manager -> CreateNtupleDCol 1082 fNtColId[2] = manager -> CreateNtupleDColum 1083 fNtColId[3] = manager -> CreateNtupleDColum 1084 1085 //Interaction position of the photon with t 1086 fNtColId[4] = manager -> CreateNtupleDColum 1087 fNtColId[5] = manager -> CreateNtupleDColum 1088 fNtColId[6] = manager -> CreateNtupleDColum 1089 1090 ////source position (annihilation position) 1091 fNtColId[7] = manager -> CreateNtupleDColum 1092 fNtColId[8] = manager -> CreateNtupleDColum 1093 fNtColId[9] = manager -> CreateNtupleDColum 1094 1095 manager -> FinishNtuple(); 1096 } 1097 if(getCoincidenceData){ 1098 manager -> CreateNtuple("Coincidence", "C 1099 fNtColId[0] = manager -> CreateNtupleIColum 1100 fNtColId[1] = manager -> CreateNtupleIColum 1101 fNtColId[2] = manager -> CreateNtupleIColum 1102 fNtColId[3] = manager -> CreateNtupleIColum 1103 fNtColId[4] = manager -> CreateNtupleIColum 1104 fNtColId[5] = manager -> CreateNtupleDColum 1105 fNtColId[6] = manager -> CreateNtupleDColum 1106 1107 fNtColId[7] = manager -> CreateNtupleIColum 1108 fNtColId[8] = manager -> CreateNtupleIColum 1109 fNtColId[9] = manager -> CreateNtupleIColum 1110 fNtColId[10] = manager -> CreateNtupleIColu 1111 fNtColId[11] = manager -> CreateNtupleIColu 1112 fNtColId[12] = manager -> CreateNtupleDColu 1113 fNtColId[13] = manager -> CreateNtupleDColu 1114 1115 //source position 1116 fNtColId[14] = manager -> CreateNtupleDColu 1117 fNtColId[15] = manager -> CreateNtupleDColu 1118 fNtColId[16] = manager -> CreateNtupleDColu 1119 1120 manager -> FinishNtuple(); 1121 1122 } 1123 1124 1125 factoryOn = true; 1126 } 1127 void doiPETAnalysis::FillListModeEvent() 1128 { 1129 1130 auto manager = G4AnalysisManager::Instance( 1131 if(getSinglesData){ 1132 manager -> FillNtupleIColumn(1, fNtColId[ 1133 manager -> FillNtupleIColumn(1, fNtColId[ 1134 //manager -> FillNtupleDColumn(1, fNtColI 1135 //manager -> FillNtupleDColumn(1, fNtColI 1136 //manager -> FillNtupleDColumn(1, fNtColI 1137 manager -> FillNtupleDColumn(1, fNtColId[ 1138 manager -> FillNtupleDColumn(1, fNtColId[ 1139 1140 //Interaction position of the photon in t 1141 manager -> FillNtupleDColumn(1, fNtColId[ 1142 manager -> FillNtupleDColumn(1, fNtColId[ 1143 manager -> FillNtupleDColumn(1, fNtColId[ 1144 1145 // 1146 //Add source position 1147 manager -> FillNtupleDColumn(1, fNtColId[ 1148 manager -> FillNtupleDColumn(1, fNtColId[ 1149 manager -> FillNtupleDColumn(1, fNtColId[ 1150 1151 manager -> AddNtupleRow(1); 1152 } 1153 1154 if(getCoincidenceData){ 1155 //First Single 1156 manager -> FillNtupleIColumn(1, fNtColId[ 1157 manager -> FillNtupleIColumn(1, fNtColId[ 1158 manager -> FillNtupleIColumn(1, fNtColId[ 1159 manager -> FillNtupleIColumn(1, fNtColId[ 1160 manager -> FillNtupleIColumn(1, fNtColId[ 1161 manager -> FillNtupleDColumn(1, fNtColId[ 1162 manager -> FillNtupleDColumn(1, fNtColId[ 1163 1164 //Second Single 1165 manager -> FillNtupleIColumn(1, fNtColId[ 1166 manager -> FillNtupleIColumn(1, fNtColId[ 1167 manager -> FillNtupleIColumn(1, fNtColId[ 1168 manager -> FillNtupleIColumn(1, fNtColId[ 1169 manager -> FillNtupleIColumn(1, fNtColId[ 1170 manager -> FillNtupleDColumn(1, fNtColId[ 1171 manager -> FillNtupleDColumn(1, fNtColId[ 1172 1173 //Add source position 1174 manager -> FillNtupleDColumn(1, fNtColId[ 1175 manager -> FillNtupleDColumn(1, fNtColId[ 1176 manager -> FillNtupleDColumn(1, fNtColId[ 1177 1178 manager -> AddNtupleRow(1); 1179 } 1180 1181 } 1182 void doiPETAnalysis::finish() 1183 { 1184 if (factoryOn) 1185 { 1186 auto manager = G4AnalysisManager::Instanc 1187 manager -> Write(); 1188 manager -> CloseFile(); 1189 1190 // delete G4AnalysisManager::Instance(); 1191 factoryOn = false; 1192 } 1193 }