#include "AnalysisManager.hh" #include "AnalysisMessenger.hh" #include "G4UnitsTable.hh" #include "G4IonTable.hh" #include "TH1D.h" #include "TH2D.h" #include "TTree.h" #include "TFile.h" AnalysisManager::AnalysisManager(): fileName("Default.root") { for(int i=0;i<2;i++) { totE_Ge[i] = 0; totEGauss_Ge[i] = 0; } MomSrcGamma[0] = G4ThreeVector(-999,-999,-999); MomSrcGamma[1] = G4ThreeVector(-999,-999,-999); for(int i=0;i<100;i++) { MomCapGamma[i] = G4ThreeVector(-999,-999,-999); } MomIon[0] = G4ThreeVector(-999,-999,-999); MomIon[1] = G4ThreeVector(-999,-999,-999); MomCapIon[0] = G4ThreeVector(-999,-999,-999); analysisMessenger = new AnalysisMessenger(this); } AnalysisManager::~AnalysisManager() { delete analysisMessenger; } void AnalysisManager::ClearVariables() { for(int i=0;i<2;i++) { totE_Ge[i] = 0; totEGauss_Ge[i] = 0; } MomSrcGamma[0] = G4ThreeVector(-999,-999,-999); MomSrcGamma[1] = G4ThreeVector(-999,-999,-999); for(int i=0;i<100;i++) { MomCapGamma[i] = G4ThreeVector(-999,-999,-999); } MomIon[0] = G4ThreeVector(-999,-999,-999); MomIon[1] = G4ThreeVector(-999,-999,-999); MomCapIon[0] = G4ThreeVector(-999,-999,-999); } void AnalysisManager::FillhGeTotEdep(G4double feDep,G4double feDepGauss,G4int fcopyNo,G4int fNofHits,G4int fparticleCode) { totE_Ge[fcopyNo] += feDep; totEGauss_Ge[fcopyNo] += feDepGauss; hCodes->Fill(fparticleCode-1E9); } void AnalysisManager::FillHistogramsAndTree() { for(int i=0;i<2;i++) { if(totE_Ge[i] != 0) hGeTotEdep[i]->Fill(totE_Ge[i]); if(totEGauss_Ge[i] != 0) hGeTotEdepGauss[i]->Fill(totEGauss_Ge[i]); } if(totE_Ge[0]!=0 || totE_Ge[1]!=0) tree1->Fill(); if(MomSrcGamma[0]!=G4ThreeVector(-999,-999,-999) && MomSrcGamma[1]!=G4ThreeVector(-999,-999,-999)) hSrcGGCosTheta->Fill(cos(MomSrcGamma[0].angle(MomSrcGamma[1]))); if(MomIon[0]!=G4ThreeVector(-999,-999,-999) && MomIon[1]!=G4ThreeVector(-999,-999,-999)) { hKinEFragIon1VSKinEFragIon2->Fill( (sqrt(MomIon[0].mag2()+DefIon[0]->GetPDGMass()*DefIon[0]->GetPDGMass())-DefIon[0]->GetPDGMass())/keV,(sqrt(MomIon[1].mag2()+DefIon[1]->GetPDGMass()*DefIon[1]->GetPDGMass())-DefIon[1]->GetPDGMass())/keV); hthetaFragIon1VSKinEFragIon1->Fill(MomIon[0].theta()/deg,(sqrt(MomIon[0].mag2()+DefIon[0]->GetPDGMass()*DefIon[0]->GetPDGMass())-DefIon[0]->GetPDGMass())/keV); hthetaFragIon2VSKinEFragIon2->Fill(MomIon[1].theta()/deg,(sqrt(MomIon[1].mag2()+DefIon[1]->GetPDGMass()*DefIon[1]->GetPDGMass())-DefIon[1]->GetPDGMass())/keV); } //G4cout << "MomCapGamma: " << MomCapGamma[0] << " MomCapIon: " << MomCapIon[0] << G4endl; if(MomCapGamma[0]!=G4ThreeVector(-999,-999,-999) && MomCapIon[0]!=G4ThreeVector(-999,-999,-999)) { hKinECapGammaVSKinECapFragIon->Fill( (sqrt(MomCapGamma[0].mag2()+G4Gamma::Gamma()->GetPDGMass()*G4Gamma::Gamma()->GetPDGMass())-G4Gamma::Gamma()->GetPDGMass())/keV,(sqrt(MomCapIon[0].mag2()+DefCapIon[0]->GetPDGMass()*DefCapIon[0]->GetPDGMass())-DefCapIon[0]->GetPDGMass())/keV); hthetaCapGammaVSKinECapGamma->Fill(MomCapGamma[0].theta()/deg,(sqrt(MomCapGamma[0].mag2()+G4Gamma::Gamma()->GetPDGMass()*G4Gamma::Gamma()->GetPDGMass())-G4Gamma::Gamma()->GetPDGMass())/keV); hthetaCapFragIonVSKinECapFragIon->Fill(MomCapIon[0].theta()/deg,(sqrt(MomCapIon[0].mag2()+DefCapIon[0]->GetPDGMass()*DefCapIon[0]->GetPDGMass())-DefCapIon[0]->GetPDGMass())/keV); } } void AnalysisManager::FillhPrimaries(G4ThreeVector fPos,G4double fkinE,G4ThreeVector fmom) { hPrimaryXvsY->Fill(fPos.x(),fPos.y()); hPrimaryZ->Fill(fPos.z()); hPrimarykinE->Fill(fkinE); hPrimaryDiv->Fill(fmom.theta()/deg); } void AnalysisManager::FillhSecondaries(G4int fNofSec, G4double *fSecKinE, G4ThreeVector *fSecMom, G4ParticleDefinition **fSecPartDef, G4double *fParKinE, G4ThreeVector *fParMom, G4ThreeVector *fParPos, G4ParticleDefinition **fParPartDef) { //G4cout << "AnalysisManager:" << G4endl; //G4cout << "NofSec: " << fNofSec << G4endl; G4int NofGamma = 0; G4int NofIon = 0; for(int i=0;iFill(fSecKinE[i]/keV); if(fParPartDef[i] == G4Proton::Proton()) hSecGammaKinEVSParKinE->Fill(fSecKinE[i]/keV,fParKinE[i]/keV); if(fParPartDef[i]->GetAtomicMass() == 22 && fParPartDef[i]->GetAtomicNumber() == 11) { // ------- here it is hParPosZVSSecKinE->Fill(fParPos[i].z(),fSecKinE[i]/keV); hParPosZ->Fill(fParPos[i].z()); // ---------------- } if(fSecKinE[i] > Gamma1-0.01 && fSecKinE[i] < Gamma1+0.01) MomSrcGamma[0] = fSecMom[i]; else if(fSecKinE[i] > Gamma2-0.01 && fSecKinE[i] < Gamma2+0.01) MomSrcGamma[1] = fSecMom[i]; MomCapGamma[NofGamma] = fSecMom[i]; NofGamma++; } else if(fSecPartDef[i]->GetAtomicMass()>0 && fSecPartDef[i]->GetAtomicNumber()>0) { if(fSecPartDef[i]->GetAtomicMass() == FragIonA1 && fSecPartDef[i]->GetAtomicNumber() == FragIonZ1) { MomIon[0] = fSecMom[i]; DefIon[0] = fSecPartDef[i]; //G4cout << "MomIon1: " << MomIon[0] << G4endl; } else if(fSecPartDef[i]->GetAtomicMass() == FragIonA2 && fSecPartDef[i]->GetAtomicNumber() == FragIonZ2) { MomIon[1] = fSecMom[i]; DefIon[1] = fSecPartDef[i]; //G4cout << "MomIon2: " << MomIon[1] << G4endl; } else if(fSecPartDef[i]->GetAtomicMass() == CapFragIonA && fSecPartDef[i]->GetAtomicNumber() == CapFragIonZ) { MomCapIon[0] = fSecMom[i]; DefCapIon[0] = fSecPartDef[i]; //G4cout << "MomIon2: " << MomIon[1] << G4endl; } NofIon++; } } } void AnalysisManager::BeginOfRun() { rootFile = new TFile(fileName,"RECREATE"); // EventAction char histname_Ge[100]; const char *prehist_Ge = "hGeTotE"; char histnameGauss_Ge[100]; const char *prehistGauss_Ge = "hGeTotEGauss"; // EventAction hCodes = new TH1D("ParticleCodes","",1000000,0,1000000); for(int i=0;i<2;i++) { std::stringstream sstri; sstri << i; std::string str1i = sstri.str(); const char *iStr = str1i.c_str(); strcpy(histname_Ge,prehist_Ge); strcat(histname_Ge,iStr); strcpy(histnameGauss_Ge,prehistGauss_Ge); strcat(histnameGauss_Ge,iStr); hGeTotEdep[i] = new TH1D(histname_Ge,"",16384,1,16384); hGeTotEdepGauss[i] = new TH1D(histnameGauss_Ge,"",16384,1,16384); } hSecGammaKinE = new TH1D("SecGammaKinE","",16384,1,16384); hSecGammaKinEVSParKinE = new TH2D("SecGammaKinEVSParKinE","",500,0,16384,500,0,6000); hParPosZVSSecKinE = new TH2D("ParPosZVSSecKinE","",500,-200,200,500,0,6000); hParPosZ = new TH1D("ParPosZ","",600,-300,300); hSrcGGCosTheta = new TH1D("SrcGGCosTheta","",1000,-1,1); hKinEFragIon1VSKinEFragIon2 = new TH2D("KinEFragIon1VSKinEFragIon2","",500,0,5000,500,0,5000); hthetaFragIon1VSKinEFragIon1 = new TH2D("thetaFragIon1VSKinEFragIon1","",500,0,180,500,0,5000); hthetaFragIon2VSKinEFragIon2 = new TH2D("thetaFragIon2VSKinEFragIon2","",500,0,180,500,0,5000); hKinECapGammaVSKinECapFragIon = new TH2D("KinECapGammaVSKinECapFragIon","",500,0,15000,500,0,5000); hthetaCapGammaVSKinECapGamma = new TH2D("thetaCapGammaVSKinECapGamma","",500,0,180,500,0,15000); hthetaCapFragIonVSKinECapFragIon = new TH2D("thetaCapFragIonVSKinECapFragIon","",500,0,180,500,0,5000); // PrimaryGeneratorAction hPrimaryXvsY = new TH2D("PrimaryXvsY","",500,-50,50,500,-50,50); hPrimaryZ = new TH1D("PrimaryZ","",2000,-1000,1000); hPrimarykinE = new TH1D("PrimarykinE","",6000,0,6000); hPrimaryDiv = new TH1D("PrimaryDiv","",100,0,10); // Tree tree1 = new TTree("Tree1", "RawData"); tree1->Branch("EdepRaw", &totE_Ge); tree1->Branch("zint", &fParPos.z()); } void AnalysisManager::EndOfRun() { if(rootFile) { rootFile->Write(); rootFile->Close(); } }