Message: Problem with GetH1()->mean() Not Logged In (login)
 Next-in-Thread Next-in-Thread
 Next-in-Forum Next-in-Forum

Question Problem with GetH1()->mean() 

Forum: Analysis
Date: 29 Sep, 2014
From: Leonardo Ghizoni <Leonardo Ghizoni>

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.79472e-06 Rad 161.353 MeV

Target 1.00005e-05 Rad 6.09191 MeV

============================================================

(from RE02)

But when I output the statistics, I get this: ---> Print Histograms Statistics

 Deposited Dose in Sphere: mean = 8.79472e-09 Rad | rms = 2.81776e-08 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.00005e-08 Rad | rms = 1.13753e-07 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:
 1 1
 All All
Outline Depth:
 1 1
 2 2
 All All
Add message: (add)

1 None: Re: Problem with GetH1()->mean()   (Tsukasa Aso - 29 Sep, 2014)
(_ Question: Re: Problem with GetH1()->mean()   (Leonardo Ghizoni - 30 Sep, 2014)
 Add Message Add Message
to: "Problem with GetH1()->mean()"

 Subscribe Subscribe

This site runs SLAC HyperNews version 1.11-slac-98, derived from the original HyperNews


[ Geant 4 Home | Geant 4 HyperNews | Search | Request New Forum | Feedback ]