Message: Problem with GetH1()>mean()  Not Logged In (login) 
Hello,
I'm building a simulation where I use the primitive scorers to get the total energy deposit and total dose deposit. I studied the RE02 example, and together with th B4d example I'm able to get the output of the sum of the events and the statistics of the events. I already studied a lot my code and cannot find the bug, but when I generate 1000 proton beams I get the total values: # DEVICE # Ddep Edep Sphere 8.79472e06 Rad 161.353 MeV Target 1.00005e05 Rad 6.09191 MeV ============================================================ (from RE02) But when I output the statistics, I get this: > Print Histograms Statistics
Deposited Dose in Sphere: mean = 8.79472e09 Rad  rms = 2.81776e08 Rad Deposited Energy in Sphere: mean = 13.1466 keV  rms = 95.42 keV Track Length in Sphere: mean = 427.324 um  rms = 1.34248 mm Deposited Dose in Target: mean = 1.00005e08 Rad  rms = 1.13753e07 Rad Deposited Energy in Target: mean = 4.64628 keV  rms = 52.1214 keV Track Length in Target: mean = 17.5213 um  rms = 193.357 um The mean of the deposited energy on the sphere should be 161.353 keV, since it was 1000 beams, but I always get this strange value. This is the piece if my code:
void SCDRunAction::BeginOfRunAction(const G4Run* run) { G4cout<<"### Run "<<run>GetRunID()+1<<" Start."<<G4endl;
// Create Analysis Manager // The choice of analysis manager is done via selection of a namespace // in SCDAnalysis.hh G4AnalysisManager* analysisManager = G4AnalysisManager::Instance(); G4cout<<"Using "<<analysisManager>GetType()<<" analysis manager"<<G4endl;
// Open an output file // G4String fileName = "SCD"; analysisManager>OpenFile(fileName); analysisManager>SetFirstHistoId(1);
// Creating Histograms // // Sphere analysisManager>CreateH1("1","Deposited Dose in Sphere",10000,0.,1*gray); analysisManager>CreateH1("2","Deposited Energy in Sphere",100,0.,1*MeV); analysisManager>CreateH1("3","Track Length in Sphere",100,0.,1*cm); // Target analysisManager>CreateH1("4","Deposited Dose in Target",10000,0.,1*gray); analysisManager>CreateH1("5","Deposited Energy in Target",100,0.,1*MeV); analysisManager>CreateH1("6","Track Length in Target",100,0.,1*cm);
// Creating nTuples // // Sphere analysisManager>CreateNtuple("SCD","Ddep, Edep and TrackLength"); analysisManager>CreateNtupleDColumn("sphereDdep"); analysisManager>CreateNtupleDColumn("sphereEdep"); analysisManager>CreateNtupleDColumn("sphereTrackLength"); // Target analysisManager>CreateNtupleDColumn("targetDdep"); analysisManager>CreateNtupleDColumn("targetEdep"); analysisManager>CreateNtupleDColumn("targetTrackLength"); analysisManager>FinishNtuple(); }
void SCDRunAction::EndOfRunAction(const G4Run* aRun) { // SCDRun object. SCDRun* scdrun = (SCDRun*)aRun; // Dump all scored quantities involved in SCDRun. scdrun>DumpAllScorer();
// // Dump accumulated quantities for this RUN. // G4THitsMap<G4double>* spheretotDdep = scdrun>GetHitsMap("sphereD/Ddep"); G4THitsMap<G4double>* spheretotEdep = scdrun>GetHitsMap("sphereD/Edep");
G4THitsMap<G4double>* targettotDdep = scdrun>GetHitsMap("targetD/Ddep"); G4THitsMap<G4double>* targettotEdep = scdrun>GetHitsMap("targetD/Edep");
G4cout<<"============================================================"<<G4endl; G4cout<<" RESULTS "<<G4endl; G4cout<<"============================================================"<<G4endl; G4cout<<std::setw(10)<<"# DEVICE #"; G4cout<<std::setw(16)<<spheretotDdep>GetName(); G4cout<<std::setw(16)<<spheretotEdep>GetName()<<G4endl;
G4double* sphereTotD = (*spheretotDdep)[0]; G4double* sphereTotE = (*spheretotEdep)[0]; G4double* targetTotD = (*targettotDdep)[0]; G4double* targetTotE = (*targettotEdep)[0]; if(!sphereTotD) sphereTotD = new G4double(0.0); if(!sphereTotE) sphereTotE = new G4double(0.0); if(!targetTotD) targetTotD = new G4double(0.0); if(!targetTotE) targetTotE = new G4double(0.0); G4cout<<std::setw( 6)<<"Sphere"<<" " <<std::setw(12)<<(*sphereTotD)*1e+14<<" Rad" <<std::setw(12)<<G4BestUnit(*sphereTotE,"Energy") <<G4endl; G4cout<<std::setw( 6)<<"Target"<<" " <<std::setw(12)<<(*targetTotD)*1e+14<<" Rad" <<std::setw(12)<<G4BestUnit(*targetTotE,"Energy") <<G4endl;
G4cout<<"============================================================"<<G4endl;
G4int nofEvents = aRun>GetNumberOfEvent(); if(nofEvents==0) return;
// Print Histogram Statistics // G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
if(analysisManager>GetH1(1)){ G4cout<<"\n > Print Histograms Statistics \n"<<G4endl; // Sphere G4cout<<" Deposited Dose in Sphere: mean = "<<(analysisManager>GetH1(1)>mean())*1e+14<<" Rad" <<"  rms = "<<(analysisManager>GetH1(1)>rms())*1e+14<<" Rad" <<G4endl; G4cout<<" Deposited Energy in Sphere: mean = "<<G4BestUnit(analysisManager>GetH1(2)>mean(),"Energy") <<"  rms = "<<G4BestUnit(analysisManager>GetH1(2)>rms(),"Energy") <<G4endl; G4cout<<" Track Length in Sphere: mean = "<<G4BestUnit(analysisManager>GetH1(3)>mean(),"Length") <<"  rms = "<<G4BestUnit(analysisManager>GetH1(3)>rms(),"Length") <<G4endl; // Target G4cout<<" Deposited Dose in Target: mean = "<<(analysisManager>GetH1(4)>mean())*1e+14<<" Rad" <<"  rms = "<<(analysisManager>GetH1(4)>rms())*1e+14<<" Rad" <<G4endl; G4cout<<" Deposited Energy in Target: mean = "<<G4BestUnit(analysisManager>GetH1(5)>mean(),"Energy") <<"  rms = "<<G4BestUnit(analysisManager>GetH1(5)>rms(),"Energy") <<G4endl; G4cout<<" Track Length in Target: mean = "<<G4BestUnit(analysisManager>GetH1(6)>mean(),"Length") <<"  rms = "<<G4BestUnit(analysisManager>GetH1(6)>rms(),"Length") <<G4endl; G4cout<<"########################################################################"<<G4endl; } // Save Histograms // analysisManager>Write(); analysisManager>CloseFile();
// Complete CleanUp // delete G4AnalysisManager::Instance(); } I also generated 100 beams and did the sum one by one to check, and the total values are the right ones. Could you please help? I don't know what I'm doing wrong here. Thank you very much. Best Regards Leonardo. ps: don't know if it changes something, but I'm using the 4.9.6 version and the ROOT analysis manager.

Inline Depth:  Outline Depth:  Add message: 
to: 