diff --git a/Examples/Example1/ShowResults.C b/Examples/Example1/ShowResults.C index 4a6f99d23035e0e6c11d8635981fcf1c4d2dcaa4..9168863302ea8510128edcb2b9a27ed53e3affe8 100644 --- a/Examples/Example1/ShowResults.C +++ b/Examples/Example1/ShowResults.C @@ -7,69 +7,69 @@ TCanvas* c1 = NULL; //////////////////////////////////////////////////////////////////////////////// void LoadCuts(){ -TFile* File_ETOF = new TFile("cuts/ETOF.root","READ"); -ETOF = (TCutG*) File_ETOF->FindObjectAny("ETOF"); - -TFile* File_EDE = new TFile("cuts/EDE.root","READ"); -EDE= (TCutG*) File_EDE->FindObjectAny("EDE"); + TFile* File_ETOF = new TFile("cuts/ETOF.root","READ"); + ETOF = (TCutG*) File_ETOF->FindObjectAny("ETOF"); + + TFile* File_EDE = new TFile("cuts/EDE.root","READ"); + EDE= (TCutG*) File_EDE->FindObjectAny("EDE"); } //////////////////////////////////////////////////////////////////////////////// void LoadChain(){ -chain = new TChain("PhysicsTree"); -chain->Add("../../Outputs/Analysis/Example1.root"); + chain = new TChain("PhysicsTree"); + chain->Add("../../Outputs/Analysis/Example1.root"); } //////////////////////////////////////////////////////////////////////////////// void ShowResults(){ -LoadChain(); -LoadCuts(); - -c1 = new TCanvas("Example1","Example1",0,0,600,600); -c1->Divide(2,2); - -// Light Particle ID // -// E-DE -c1->cd(1); -chain->Draw("SSSD.Energy:MUST2.Si_E>>hIDE(1000,0,35,1000,0,5)","MUST2.CsI_E<0 && MUST2.TelescopeNumber<5","colz"); -EDE->Draw("same"); - -// E-TOF -c1->cd(2); -chain->Draw("-MUST2.Si_T:SSSD.Energy+MUST2.Si_E>>hIDT(1000,0,35,1000,-15,0)","MUST2.CsI_E<0 && MUST2.TelescopeNumber<5","colz"); -ETOF->Draw("same"); - -// Kinematical Line // -c1->cd(3); -//chain->Draw("ELab:ThetaLab>>hKine(500,0,45,400,0,40)","MUST2.CsI_E<0 && MUST2.TelescopeNumber<5 && EDE && ETOF","colz"); -chain->Draw("ELab:ThetaLab>>h(1000,0,90,1000,0,30)","MUST2.CsI_E<0 && MUST2.TelescopeNumber<5 && EDE && ETOF","colz"); - -NPL::Reaction r("11Li(d,3He)10He@553"); -r.SetExcitationHeavy(1.4); -TGraph* Kine = r.GetKinematicLine3(); -Kine->SetLineWidth(2); -Kine->SetLineColor(kOrange-3); -Kine->Draw("c"); - -// Excitation Energy // -c1->cd(4); -int bin=50; -double Emin = -5; -double Emax = 5; - -chain->Draw(Form("Ex>>hEx(%d,%f,%f)",bin,Emin,Emax),"MUST2.CsI_E<0 && MUST2.TelescopeNumber<5 && EDE && ETOF"); -TH1F* hEx = (TH1F*) gDirectory->FindObjectAny("hEx"); -hEx->GetYaxis()->SetTitle(Form("counts / %d keV",(int) (1000*(Emax-Emin)/bin))); -hEx->GetXaxis()->SetTitle("E_{10He}"); -hEx->SetFillStyle(1001); -hEx->SetLineColor(kAzure+7); -hEx->SetFillColor(kAzure+7); - -hEx->Fit("gaus"); - -TF1* f = hEx->GetFunction("gaus"); -f->SetLineWidth(2); -f->SetLineColor(kOrange-3); -f->SetNpx(1000); - + LoadChain(); + LoadCuts(); + + c1 = new TCanvas("Example1","Example1",0,0,600,600); + c1->Divide(2,2); + + // Light Particle ID // + // E-DE + c1->cd(1); + chain->Draw("SSSD.Energy:MUST2.Si_E>>hIDE(1000,0,35,1000,0,5)","MUST2.CsI_E<0 && MUST2.TelescopeNumber<5","colz"); + EDE->Draw("same"); + + // E-TOF + c1->cd(2); + chain->Draw("-MUST2.Si_T:SSSD.Energy+MUST2.Si_E>>hIDT(1000,0,35,1000,-15,0)","MUST2.CsI_E<0 && MUST2.TelescopeNumber<5","colz"); + ETOF->Draw("same"); + + // Kinematical Line // + c1->cd(3); + //chain->Draw("ELab:ThetaLab>>hKine(500,0,45,400,0,40)","MUST2.CsI_E<0 && MUST2.TelescopeNumber<5 && EDE && ETOF","colz"); + chain->Draw("ELab:ThetaLab>>h(1000,0,90,1000,0,30)","MUST2.CsI_E<0 && MUST2.TelescopeNumber<5 && EDE && ETOF","colz"); + + NPL::Reaction r("11Li(d,3He)10He@553"); + r.SetExcitationHeavy(1.4); + TGraph* Kine = r.GetKinematicLine3(); + Kine->SetLineWidth(2); + Kine->SetLineColor(kOrange-3); + Kine->Draw("c"); + + // Excitation Energy // + c1->cd(4); + int bin=50; + double Emin = -5; + double Emax = 5; + + chain->Draw(Form("Ex>>hEx(%d,%f,%f)",bin,Emin,Emax),"MUST2.CsI_E<0 && MUST2.TelescopeNumber<5 && EDE && ETOF"); + TH1F* hEx = (TH1F*) gDirectory->FindObjectAny("hEx"); + hEx->GetYaxis()->SetTitle(Form("counts / %d keV",(int) (1000*(Emax-Emin)/bin))); + hEx->GetXaxis()->SetTitle("E_{10He}"); + hEx->SetFillStyle(1001); + hEx->SetLineColor(kAzure+7); + hEx->SetFillColor(kAzure+7); + + hEx->Fit("gaus"); + + TF1* f = hEx->GetFunction("gaus"); + f->SetLineWidth(2); + f->SetLineColor(kOrange-3); + f->SetNpx(1000); + } diff --git a/Inputs/CrossSection/flat-0_20.txt b/Inputs/CrossSection/flat-0_20.txt new file mode 100644 index 0000000000000000000000000000000000000000..b1f1d9cfe6591bb2b710170cd0675cc5845b7300 --- /dev/null +++ b/Inputs/CrossSection/flat-0_20.txt @@ -0,0 +1,1801 @@ +0 1 +0.1 1 +0.2 1 +0.3 1 +0.4 1 +0.5 1 +0.6 1 +0.7 1 +0.8 1 +0.9 1 +1 1 +1.1 1 +1.2 1 +1.3 1 +1.4 1 +1.5 1 +1.6 1 +1.7 1 +1.8 1 +1.9 1 +2 1 +2.1 1 +2.2 1 +2.3 1 +2.4 1 +2.5 1 +2.6 1 +2.7 1 +2.8 1 +2.9 1 +3 1 +3.1 1 +3.2 1 +3.3 1 +3.4 1 +3.5 1 +3.6 1 +3.7 1 +3.8 1 +3.9 1 +4 1 +4.1 1 +4.2 1 +4.3 1 +4.4 1 +4.5 1 +4.6 1 +4.7 1 +4.8 1 +4.9 1 +5 1 +5.1 1 +5.2 1 +5.3 1 +5.4 1 +5.5 1 +5.6 1 +5.7 1 +5.8 1 +5.9 1 +6 1 +6.1 1 +6.2 1 +6.3 1 +6.4 1 +6.5 1 +6.6 1 +6.7 1 +6.8 1 +6.9 1 +7 1 +7.1 1 +7.2 1 +7.3 1 +7.4 1 +7.5 1 +7.6 1 +7.7 1 +7.8 1 +7.9 1 +8 1 +8.1 1 +8.2 1 +8.3 1 +8.4 1 +8.5 1 +8.6 1 +8.7 1 +8.8 1 +8.9 1 +9 1 +9.1 1 +9.2 1 +9.3 1 +9.4 1 +9.5 1 +9.6 1 +9.7 1 +9.8 1 +9.9 1 +10 1 +10.1 1 +10.2 1 +10.3 1 +10.4 1 +10.5 1 +10.6 1 +10.7 1 +10.8 1 +10.9 1 +11 1 +11.1 1 +11.2 1 +11.3 1 +11.4 1 +11.5 1 +11.6 1 +11.7 1 +11.8 1 +11.9 1 +12 1 +12.1 1 +12.2 1 +12.3 1 +12.4 1 +12.5 1 +12.6 1 +12.7 1 +12.8 1 +12.9 1 +13 1 +13.1 1 +13.2 1 +13.3 1 +13.4 1 +13.5 1 +13.6 1 +13.7 1 +13.8 1 +13.9 1 +14 1 +14.1 1 +14.2 1 +14.3 1 +14.4 1 +14.5 1 +14.6 1 +14.7 1 +14.8 1 +14.9 1 +15 1 +15.1 1 +15.2 1 +15.3 1 +15.4 1 +15.5 1 +15.6 1 +15.7 1 +15.8 1 +15.9 1 +16 1 +16.1 1 +16.2 1 +16.3 1 +16.4 1 +16.5 1 +16.6 1 +16.7 1 +16.8 1 +16.9 1 +17 1 +17.1 1 +17.2 1 +17.3 1 +17.4 1 +17.5 1 +17.6 1 +17.7 1 +17.8 1 +17.9 1 +18 1 +18.1 1 +18.2 1 +18.3 1 +18.4 1 +18.5 1 +18.6 1 +18.7 1 +18.8 1 +18.9 1 +19 1 +19.1 1 +19.2 1 +19.3 1 +19.4 1 +19.5 1 +19.6 1 +19.7 1 +19.8 1 +19.9 1 +20 1 +20.1 0 +20.2 0 +20.3 0 +20.4 0 +20.5 0 +20.6 0 +20.7 0 +20.8 0 +20.9 0 +21 0 +21.1 0 +21.2 0 +21.3 0 +21.4 0 +21.5 0 +21.6 0 +21.7 0 +21.8 0 +21.9 0 +22 0 +22.1 0 +22.2 0 +22.3 0 +22.4 0 +22.5 0 +22.6 0 +22.7 0 +22.8 0 +22.9 0 +23 0 +23.1 0 +23.2 0 +23.3 0 +23.4 0 +23.5 0 +23.6 0 +23.7 0 +23.8 0 +23.9 0 +24 0 +24.1 0 +24.2 0 +24.3 0 +24.4 0 +24.5 0 +24.6 0 +24.7 0 +24.8 0 +24.9 0 +25 0 +25.1 0 +25.2 0 +25.3 0 +25.4 0 +25.5 0 +25.6 0 +25.7 0 +25.8 0 +25.9 0 +26 0 +26.1 0 +26.2 0 +26.3 0 +26.4 0 +26.5 0 +26.6 0 +26.7 0 +26.8 0 +26.9 0 +27 0 +27.1 0 +27.2 0 +27.3 0 +27.4 0 +27.5 0 +27.6 0 +27.7 0 +27.8 0 +27.9 0 +28 0 +28.1 0 +28.2 0 +28.3 0 +28.4 0 +28.5 0 +28.6 0 +28.7 0 +28.8 0 +28.9 0 +29 0 +29.1 0 +29.2 0 +29.3 0 +29.4 0 +29.5 0 +29.6 0 +29.7 0 +29.8 0 +29.9 0 +30 0 +30.1 0 +30.2 0 +30.3 0 +30.4 0 +30.5 0 +30.6 0 +30.7 0 +30.8 0 +30.9 0 +31 0 +31.1 0 +31.2 0 +31.3 0 +31.4 0 +31.5 0 +31.6 0 +31.7 0 +31.8 0 +31.9 0 +32 0 +32.1 0 +32.2 0 +32.3 0 +32.4 0 +32.5 0 +32.6 0 +32.7 0 +32.8 0 +32.9 0 +33 0 +33.1 0 +33.2 0 +33.3 0 +33.4 0 +33.5 0 +33.6 0 +33.7 0 +33.8 0 +33.9 0 +34 0 +34.1 0 +34.2 0 +34.3 0 +34.4 0 +34.5 0 +34.6 0 +34.7 0 +34.8 0 +34.9 0 +35 0 +35.1 0 +35.2 0 +35.3 0 +35.4 0 +35.5 0 +35.6 0 +35.7 0 +35.8 0 +35.9 0 +36 0 +36.1 0 +36.2 0 +36.3 0 +36.4 0 +36.5 0 +36.6 0 +36.7 0 +36.8 0 +36.9 0 +37 0 +37.1 0 +37.2 0 +37.3 0 +37.4 0 +37.5 0 +37.6 0 +37.7 0 +37.8 0 +37.9 0 +38 0 +38.1 0 +38.2 0 +38.3 0 +38.4 0 +38.5 0 +38.6 0 +38.7 0 +38.8 0 +38.9 0 +39 0 +39.1 0 +39.2 0 +39.3 0 +39.4 0 +39.5 0 +39.6 0 +39.7 0 +39.8 0 +39.9 0 +40 0 +40.1 0 +40.2 0 +40.3 0 +40.4 0 +40.5 0 +40.6 0 +40.7 0 +40.8 0 +40.9 0 +41 0 +41.1 0 +41.2 0 +41.3 0 +41.4 0 +41.5 0 +41.6 0 +41.7 0 +41.8 0 +41.9 0 +42 0 +42.1 0 +42.2 0 +42.3 0 +42.4 0 +42.5 0 +42.6 0 +42.7 0 +42.8 0 +42.9 0 +43 0 +43.1 0 +43.2 0 +43.3 0 +43.4 0 +43.5 0 +43.6 0 +43.7 0 +43.8 0 +43.9 0 +44 0 +44.1 0 +44.2 0 +44.3 0 +44.4 0 +44.5 0 +44.6 0 +44.7 0 +44.8 0 +44.9 0 +45 0 +45.1 0 +45.2 0 +45.3 0 +45.4 0 +45.5 0 +45.6 0 +45.7 0 +45.8 0 +45.9 0 +46 0 +46.1 0 +46.2 0 +46.3 0 +46.4 0 +46.5 0 +46.6 0 +46.7 0 +46.8 0 +46.9 0 +47 0 +47.1 0 +47.2 0 +47.3 0 +47.4 0 +47.5 0 +47.6 0 +47.7 0 +47.8 0 +47.9 0 +48 0 +48.1 0 +48.2 0 +48.3 0 +48.4 0 +48.5 0 +48.6 0 +48.7 0 +48.8 0 +48.9 0 +49 0 +49.1 0 +49.2 0 +49.3 0 +49.4 0 +49.5 0 +49.6 0 +49.7 0 +49.8 0 +49.9 0 +50 0 +50.1 0 +50.2 0 +50.3 0 +50.4 0 +50.5 0 +50.6 0 +50.7 0 +50.8 0 +50.9 0 +51 0 +51.1 0 +51.2 0 +51.3 0 +51.4 0 +51.5 0 +51.6 0 +51.7 0 +51.8 0 +51.9 0 +52 0 +52.1 0 +52.2 0 +52.3 0 +52.4 0 +52.5 0 +52.6 0 +52.7 0 +52.8 0 +52.9 0 +53 0 +53.1 0 +53.2 0 +53.3 0 +53.4 0 +53.5 0 +53.6 0 +53.7 0 +53.8 0 +53.9 0 +54 0 +54.1 0 +54.2 0 +54.3 0 +54.4 0 +54.5 0 +54.6 0 +54.7 0 +54.8 0 +54.9 0 +55 0 +55.1 0 +55.2 0 +55.3 0 +55.4 0 +55.5 0 +55.6 0 +55.7 0 +55.8 0 +55.9 0 +56 0 +56.1 0 +56.2 0 +56.3 0 +56.4 0 +56.5 0 +56.6 0 +56.7 0 +56.8 0 +56.9 0 +57 0 +57.1 0 +57.2 0 +57.3 0 +57.4 0 +57.5 0 +57.6 0 +57.7 0 +57.8 0 +57.9 0 +58 0 +58.1 0 +58.2 0 +58.3 0 +58.4 0 +58.5 0 +58.6 0 +58.7 0 +58.8 0 +58.9 0 +59 0 +59.1 0 +59.2 0 +59.3 0 +59.4 0 +59.5 0 +59.6 0 +59.7 0 +59.8 0 +59.9 0 +60 0 +60.1 0 +60.2 0 +60.3 0 +60.4 0 +60.5 0 +60.6 0 +60.7 0 +60.8 0 +60.9 0 +61 0 +61.1 0 +61.2 0 +61.3 0 +61.4 0 +61.5 0 +61.6 0 +61.7 0 +61.8 0 +61.9 0 +62 0 +62.1 0 +62.2 0 +62.3 0 +62.4 0 +62.5 0 +62.6 0 +62.7 0 +62.8 0 +62.9 0 +63 0 +63.1 0 +63.2 0 +63.3 0 +63.4 0 +63.5 0 +63.6 0 +63.7 0 +63.8 0 +63.9 0 +64 0 +64.1 0 +64.2 0 +64.3 0 +64.4 0 +64.5 0 +64.6 0 +64.7 0 +64.8 0 +64.9 0 +65 0 +65.1 0 +65.2 0 +65.3 0 +65.4 0 +65.5 0 +65.6 0 +65.7 0 +65.8 0 +65.9 0 +66 0 +66.1 0 +66.2 0 +66.3 0 +66.4 0 +66.5 0 +66.6 0 +66.7 0 +66.8 0 +66.9 0 +67 0 +67.1 0 +67.2 0 +67.3 0 +67.4 0 +67.5 0 +67.6 0 +67.7 0 +67.8 0 +67.9 0 +68 0 +68.1 0 +68.2 0 +68.3 0 +68.4 0 +68.5 0 +68.6 0 +68.7 0 +68.8 0 +68.9 0 +69 0 +69.1 0 +69.2 0 +69.3 0 +69.4 0 +69.5 0 +69.6 0 +69.7 0 +69.8 0 +69.9 0 +70 0 +70.1 0 +70.2 0 +70.3 0 +70.4 0 +70.5 0 +70.6 0 +70.7 0 +70.8 0 +70.9 0 +71 0 +71.1 0 +71.2 0 +71.3 0 +71.4 0 +71.5 0 +71.6 0 +71.7 0 +71.8 0 +71.9 0 +72 0 +72.1 0 +72.2 0 +72.3 0 +72.4 0 +72.5 0 +72.6 0 +72.7 0 +72.8 0 +72.9 0 +73 0 +73.1 0 +73.2 0 +73.3 0 +73.4 0 +73.5 0 +73.6 0 +73.7 0 +73.8 0 +73.9 0 +74 0 +74.1 0 +74.2 0 +74.3 0 +74.4 0 +74.5 0 +74.6 0 +74.7 0 +74.8 0 +74.9 0 +75 0 +75.1 0 +75.2 0 +75.3 0 +75.4 0 +75.5 0 +75.6 0 +75.7 0 +75.8 0 +75.9 0 +76 0 +76.1 0 +76.2 0 +76.3 0 +76.4 0 +76.5 0 +76.6 0 +76.7 0 +76.8 0 +76.9 0 +77 0 +77.1 0 +77.2 0 +77.3 0 +77.4 0 +77.5 0 +77.6 0 +77.7 0 +77.8 0 +77.9 0 +78 0 +78.1 0 +78.2 0 +78.3 0 +78.4 0 +78.5 0 +78.6 0 +78.7 0 +78.8 0 +78.9 0 +79 0 +79.1 0 +79.2 0 +79.3 0 +79.4 0 +79.5 0 +79.6 0 +79.7 0 +79.8 0 +79.9 0 +80 0 +80.1 0 +80.2 0 +80.3 0 +80.4 0 +80.5 0 +80.6 0 +80.7 0 +80.8 0 +80.9 0 +81 0 +81.1 0 +81.2 0 +81.3 0 +81.4 0 +81.5 0 +81.6 0 +81.7 0 +81.8 0 +81.9 0 +82 0 +82.1 0 +82.2 0 +82.3 0 +82.4 0 +82.5 0 +82.6 0 +82.7 0 +82.8 0 +82.9 0 +83 0 +83.1 0 +83.2 0 +83.3 0 +83.4 0 +83.5 0 +83.6 0 +83.7 0 +83.8 0 +83.9 0 +84 0 +84.1 0 +84.2 0 +84.3 0 +84.4 0 +84.5 0 +84.6 0 +84.7 0 +84.8 0 +84.9 0 +85 0 +85.1 0 +85.2 0 +85.3 0 +85.4 0 +85.5 0 +85.6 0 +85.7 0 +85.8 0 +85.9 0 +86 0 +86.1 0 +86.2 0 +86.3 0 +86.4 0 +86.5 0 +86.6 0 +86.7 0 +86.8 0 +86.9 0 +87 0 +87.1 0 +87.2 0 +87.3 0 +87.4 0 +87.5 0 +87.6 0 +87.7 0 +87.8 0 +87.9 0 +88 0 +88.1 0 +88.2 0 +88.3 0 +88.4 0 +88.5 0 +88.6 0 +88.7 0 +88.8 0 +88.9 0 +89 0 +89.1 0 +89.2 0 +89.3 0 +89.4 0 +89.5 0 +89.6 0 +89.7 0 +89.8 0 +89.9 0 +90 0 +90.1 0 +90.2 0 +90.3 0 +90.4 0 +90.5 0 +90.6 0 +90.7 0 +90.8 0 +90.9 0 +91 0 +91.1 0 +91.2 0 +91.3 0 +91.4 0 +91.5 0 +91.6 0 +91.7 0 +91.8 0 +91.9 0 +92 0 +92.1 0 +92.2 0 +92.3 0 +92.4 0 +92.5 0 +92.6 0 +92.7 0 +92.8 0 +92.9 0 +93 0 +93.1 0 +93.2 0 +93.3 0 +93.4 0 +93.5 0 +93.6 0 +93.7 0 +93.8 0 +93.9 0 +94 0 +94.1 0 +94.2 0 +94.3 0 +94.4 0 +94.5 0 +94.6 0 +94.7 0 +94.8 0 +94.9 0 +95 0 +95.1 0 +95.2 0 +95.3 0 +95.4 0 +95.5 0 +95.6 0 +95.7 0 +95.8 0 +95.9 0 +96 0 +96.1 0 +96.2 0 +96.3 0 +96.4 0 +96.5 0 +96.6 0 +96.7 0 +96.8 0 +96.9 0 +97 0 +97.1 0 +97.2 0 +97.3 0 +97.4 0 +97.5 0 +97.6 0 +97.7 0 +97.8 0 +97.9 0 +98 0 +98.1 0 +98.2 0 +98.3 0 +98.4 0 +98.5 0 +98.6 0 +98.7 0 +98.8 0 +98.9 0 +99 0 +99.1 0 +99.2 0 +99.3 0 +99.4 0 +99.5 0 +99.6 0 +99.7 0 +99.8 0 +99.9 0 +100 0 +100.1 0 +100.2 0 +100.3 0 +100.4 0 +100.5 0 +100.6 0 +100.7 0 +100.8 0 +100.9 0 +101 0 +101.1 0 +101.2 0 +101.3 0 +101.4 0 +101.5 0 +101.6 0 +101.7 0 +101.8 0 +101.9 0 +102 0 +102.1 0 +102.2 0 +102.3 0 +102.4 0 +102.5 0 +102.6 0 +102.7 0 +102.8 0 +102.9 0 +103 0 +103.1 0 +103.2 0 +103.3 0 +103.4 0 +103.5 0 +103.6 0 +103.7 0 +103.8 0 +103.9 0 +104 0 +104.1 0 +104.2 0 +104.3 0 +104.4 0 +104.5 0 +104.6 0 +104.7 0 +104.8 0 +104.9 0 +105 0 +105.1 0 +105.2 0 +105.3 0 +105.4 0 +105.5 0 +105.6 0 +105.7 0 +105.8 0 +105.9 0 +106 0 +106.1 0 +106.2 0 +106.3 0 +106.4 0 +106.5 0 +106.6 0 +106.7 0 +106.8 0 +106.9 0 +107 0 +107.1 0 +107.2 0 +107.3 0 +107.4 0 +107.5 0 +107.6 0 +107.7 0 +107.8 0 +107.9 0 +108 0 +108.1 0 +108.2 0 +108.3 0 +108.4 0 +108.5 0 +108.6 0 +108.7 0 +108.8 0 +108.9 0 +109 0 +109.1 0 +109.2 0 +109.3 0 +109.4 0 +109.5 0 +109.6 0 +109.7 0 +109.8 0 +109.9 0 +110 0 +110.1 0 +110.2 0 +110.3 0 +110.4 0 +110.5 0 +110.6 0 +110.7 0 +110.8 0 +110.9 0 +111 0 +111.1 0 +111.2 0 +111.3 0 +111.4 0 +111.5 0 +111.6 0 +111.7 0 +111.8 0 +111.9 0 +112 0 +112.1 0 +112.2 0 +112.3 0 +112.4 0 +112.5 0 +112.6 0 +112.7 0 +112.8 0 +112.9 0 +113 0 +113.1 0 +113.2 0 +113.3 0 +113.4 0 +113.5 0 +113.6 0 +113.7 0 +113.8 0 +113.9 0 +114 0 +114.1 0 +114.2 0 +114.3 0 +114.4 0 +114.5 0 +114.6 0 +114.7 0 +114.8 0 +114.9 0 +115 0 +115.1 0 +115.2 0 +115.3 0 +115.4 0 +115.5 0 +115.6 0 +115.7 0 +115.8 0 +115.9 0 +116 0 +116.1 0 +116.2 0 +116.3 0 +116.4 0 +116.5 0 +116.6 0 +116.7 0 +116.8 0 +116.9 0 +117 0 +117.1 0 +117.2 0 +117.3 0 +117.4 0 +117.5 0 +117.6 0 +117.7 0 +117.8 0 +117.9 0 +118 0 +118.1 0 +118.2 0 +118.3 0 +118.4 0 +118.5 0 +118.6 0 +118.7 0 +118.8 0 +118.9 0 +119 0 +119.1 0 +119.2 0 +119.3 0 +119.4 0 +119.5 0 +119.6 0 +119.7 0 +119.8 0 +119.9 0 +120 0 +120.1 0 +120.2 0 +120.3 0 +120.4 0 +120.5 0 +120.6 0 +120.7 0 +120.8 0 +120.9 0 +121 0 +121.1 0 +121.2 0 +121.3 0 +121.4 0 +121.5 0 +121.6 0 +121.7 0 +121.8 0 +121.9 0 +122 0 +122.1 0 +122.2 0 +122.3 0 +122.4 0 +122.5 0 +122.6 0 +122.7 0 +122.8 0 +122.9 0 +123 0 +123.1 0 +123.2 0 +123.3 0 +123.4 0 +123.5 0 +123.6 0 +123.7 0 +123.8 0 +123.9 0 +124 0 +124.1 0 +124.2 0 +124.3 0 +124.4 0 +124.5 0 +124.6 0 +124.7 0 +124.8 0 +124.9 0 +125 0 +125.1 0 +125.2 0 +125.3 0 +125.4 0 +125.5 0 +125.6 0 +125.7 0 +125.8 0 +125.9 0 +126 0 +126.1 0 +126.2 0 +126.3 0 +126.4 0 +126.5 0 +126.6 0 +126.7 0 +126.8 0 +126.9 0 +127 0 +127.1 0 +127.2 0 +127.3 0 +127.4 0 +127.5 0 +127.6 0 +127.7 0 +127.8 0 +127.9 0 +128 0 +128.1 0 +128.2 0 +128.3 0 +128.4 0 +128.5 0 +128.6 0 +128.7 0 +128.8 0 +128.9 0 +129 0 +129.1 0 +129.2 0 +129.3 0 +129.4 0 +129.5 0 +129.6 0 +129.7 0 +129.8 0 +129.9 0 +130 0 +130.1 0 +130.2 0 +130.3 0 +130.4 0 +130.5 0 +130.6 0 +130.7 0 +130.8 0 +130.9 0 +131 0 +131.1 0 +131.2 0 +131.3 0 +131.4 0 +131.5 0 +131.6 0 +131.7 0 +131.8 0 +131.9 0 +132 0 +132.1 0 +132.2 0 +132.3 0 +132.4 0 +132.5 0 +132.6 0 +132.7 0 +132.8 0 +132.9 0 +133 0 +133.1 0 +133.2 0 +133.3 0 +133.4 0 +133.5 0 +133.6 0 +133.7 0 +133.8 0 +133.9 0 +134 0 +134.1 0 +134.2 0 +134.3 0 +134.4 0 +134.5 0 +134.6 0 +134.7 0 +134.8 0 +134.9 0 +135 0 +135.1 0 +135.2 0 +135.3 0 +135.4 0 +135.5 0 +135.6 0 +135.7 0 +135.8 0 +135.9 0 +136 0 +136.1 0 +136.2 0 +136.3 0 +136.4 0 +136.5 0 +136.6 0 +136.7 0 +136.8 0 +136.9 0 +137 0 +137.1 0 +137.2 0 +137.3 0 +137.4 0 +137.5 0 +137.6 0 +137.7 0 +137.8 0 +137.9 0 +138 0 +138.1 0 +138.2 0 +138.3 0 +138.4 0 +138.5 0 +138.6 0 +138.7 0 +138.8 0 +138.9 0 +139 0 +139.1 0 +139.2 0 +139.3 0 +139.4 0 +139.5 0 +139.6 0 +139.7 0 +139.8 0 +139.9 0 +140 0 +140.1 0 +140.2 0 +140.3 0 +140.4 0 +140.5 0 +140.6 0 +140.7 0 +140.8 0 +140.9 0 +141 0 +141.1 0 +141.2 0 +141.3 0 +141.4 0 +141.5 0 +141.6 0 +141.7 0 +141.8 0 +141.9 0 +142 0 +142.1 0 +142.2 0 +142.3 0 +142.4 0 +142.5 0 +142.6 0 +142.7 0 +142.8 0 +142.9 0 +143 0 +143.1 0 +143.2 0 +143.3 0 +143.4 0 +143.5 0 +143.6 0 +143.7 0 +143.8 0 +143.9 0 +144 0 +144.1 0 +144.2 0 +144.3 0 +144.4 0 +144.5 0 +144.6 0 +144.7 0 +144.8 0 +144.9 0 +145 0 +145.1 0 +145.2 0 +145.3 0 +145.4 0 +145.5 0 +145.6 0 +145.7 0 +145.8 0 +145.9 0 +146 0 +146.1 0 +146.2 0 +146.3 0 +146.4 0 +146.5 0 +146.6 0 +146.7 0 +146.8 0 +146.9 0 +147 0 +147.1 0 +147.2 0 +147.3 0 +147.4 0 +147.5 0 +147.6 0 +147.7 0 +147.8 0 +147.9 0 +148 0 +148.1 0 +148.2 0 +148.3 0 +148.4 0 +148.5 0 +148.6 0 +148.7 0 +148.8 0 +148.9 0 +149 0 +149.1 0 +149.2 0 +149.3 0 +149.4 0 +149.5 0 +149.6 0 +149.7 0 +149.8 0 +149.9 0 +150 0 +150.1 0 +150.2 0 +150.3 0 +150.4 0 +150.5 0 +150.6 0 +150.7 0 +150.8 0 +150.9 0 +151 0 +151.1 0 +151.2 0 +151.3 0 +151.4 0 +151.5 0 +151.6 0 +151.7 0 +151.8 0 +151.9 0 +152 0 +152.1 0 +152.2 0 +152.3 0 +152.4 0 +152.5 0 +152.6 0 +152.7 0 +152.8 0 +152.9 0 +153 0 +153.1 0 +153.2 0 +153.3 0 +153.4 0 +153.5 0 +153.6 0 +153.7 0 +153.8 0 +153.9 0 +154 0 +154.1 0 +154.2 0 +154.3 0 +154.4 0 +154.5 0 +154.6 0 +154.7 0 +154.8 0 +154.9 0 +155 0 +155.1 0 +155.2 0 +155.3 0 +155.4 0 +155.5 0 +155.6 0 +155.7 0 +155.8 0 +155.9 0 +156 0 +156.1 0 +156.2 0 +156.3 0 +156.4 0 +156.5 0 +156.6 0 +156.7 0 +156.8 0 +156.9 0 +157 0 +157.1 0 +157.2 0 +157.3 0 +157.4 0 +157.5 0 +157.6 0 +157.7 0 +157.8 0 +157.9 0 +158 0 +158.1 0 +158.2 0 +158.3 0 +158.4 0 +158.5 0 +158.6 0 +158.7 0 +158.8 0 +158.9 0 +159 0 +159.1 0 +159.2 0 +159.3 0 +159.4 0 +159.5 0 +159.6 0 +159.7 0 +159.8 0 +159.9 0 +160 0 +160.1 0 +160.2 0 +160.3 0 +160.4 0 +160.5 0 +160.6 0 +160.7 0 +160.8 0 +160.9 0 +161 0 +161.1 0 +161.2 0 +161.3 0 +161.4 0 +161.5 0 +161.6 0 +161.7 0 +161.8 0 +161.9 0 +162 0 +162.1 0 +162.2 0 +162.3 0 +162.4 0 +162.5 0 +162.6 0 +162.7 0 +162.8 0 +162.9 0 +163 0 +163.1 0 +163.2 0 +163.3 0 +163.4 0 +163.5 0 +163.6 0 +163.7 0 +163.8 0 +163.9 0 +164 0 +164.1 0 +164.2 0 +164.3 0 +164.4 0 +164.5 0 +164.6 0 +164.7 0 +164.8 0 +164.9 0 +165 0 +165.1 0 +165.2 0 +165.3 0 +165.4 0 +165.5 0 +165.6 0 +165.7 0 +165.8 0 +165.9 0 +166 0 +166.1 0 +166.2 0 +166.3 0 +166.4 0 +166.5 0 +166.6 0 +166.7 0 +166.8 0 +166.9 0 +167 0 +167.1 0 +167.2 0 +167.3 0 +167.4 0 +167.5 0 +167.6 0 +167.7 0 +167.8 0 +167.9 0 +168 0 +168.1 0 +168.2 0 +168.3 0 +168.4 0 +168.5 0 +168.6 0 +168.7 0 +168.8 0 +168.9 0 +169 0 +169.1 0 +169.2 0 +169.3 0 +169.4 0 +169.5 0 +169.6 0 +169.7 0 +169.8 0 +169.9 0 +170 0 +170.1 0 +170.2 0 +170.3 0 +170.4 0 +170.5 0 +170.6 0 +170.7 0 +170.8 0 +170.9 0 +171 0 +171.1 0 +171.2 0 +171.3 0 +171.4 0 +171.5 0 +171.6 0 +171.7 0 +171.8 0 +171.9 0 +172 0 +172.1 0 +172.2 0 +172.3 0 +172.4 0 +172.5 0 +172.6 0 +172.7 0 +172.8 0 +172.9 0 +173 0 +173.1 0 +173.2 0 +173.3 0 +173.4 0 +173.5 0 +173.6 0 +173.7 0 +173.8 0 +173.9 0 +174 0 +174.1 0 +174.2 0 +174.3 0 +174.4 0 +174.5 0 +174.6 0 +174.7 0 +174.8 0 +174.9 0 +175 0 +175.1 0 +175.2 0 +175.3 0 +175.4 0 +175.5 0 +175.6 0 +175.7 0 +175.8 0 +175.9 0 +176 0 +176.1 0 +176.2 0 +176.3 0 +176.4 0 +176.5 0 +176.6 0 +176.7 0 +176.8 0 +176.9 0 +177 0 +177.1 0 +177.2 0 +177.3 0 +177.4 0 +177.5 0 +177.6 0 +177.7 0 +177.8 0 +177.9 0 +178 0 +178.1 0 +178.2 0 +178.3 0 +178.4 0 +178.5 0 +178.6 0 +178.7 0 +178.8 0 +178.9 0 +179 0 +179.1 0 +179.2 0 +179.3 0 +179.4 0 +179.5 0 +179.6 0 +179.7 0 +179.8 0 +179.9 0 +180 0 diff --git a/Inputs/DetectorConfiguration/actar_test.detector b/Inputs/DetectorConfiguration/actar_test.detector new file mode 100644 index 0000000000000000000000000000000000000000..c46dd5f1d2c540b66ccdc643207035dff16ce266 --- /dev/null +++ b/Inputs/DetectorConfiguration/actar_test.detector @@ -0,0 +1,19 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +Target + THICKNESS= 1 micrometer + ANGLE= 0 deg + RADIUS= 7.5 mm + MATERIAL= CH2 + X= 0 mm + Y= 0 mm + Z= -110 mm +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +Actar + POS= 0 0 0 mm + Shape= Square + Si= 1 + CsI= 1 + GasMaterial= iC4H10 + GasFraction= 100 + Temperature= 295 kelvin + Pressure= 0.4 bar diff --git a/Inputs/EventGenerator/18Opp.reaction b/Inputs/EventGenerator/18Opp.reaction new file mode 100644 index 0000000000000000000000000000000000000000..488eed07eb95c1528e653d8bda9377052c44ecc1 --- /dev/null +++ b/Inputs/EventGenerator/18Opp.reaction @@ -0,0 +1,28 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%% Reaction file for 11Li(d,3He)10He reaction %%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%Beam energy given in MeV ; Excitation in MeV +Beam + Particle= 18O + Energy= 137.88 + SigmaEnergy= 0 + SigmaThetaX= 0 + SigmaPhiY= 0 + SigmaX= 0 + SigmaY= 0 + MeanThetaX= 0 + MeanPhiY= 0 + MeanX= 0 + MeanY= 0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +TwoBodyReaction + Beam= 18O + Target= 1H + Light= 1H + Heavy= 18O + ExcitationEnergyLight= 0.0 + ExcitationEnergyHeavy= 0.0 + CrossSectionPath= flat.txt CS + ShootLight= 1 + ShootHeavy= 1 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/Inputs/EventGenerator/alpha.source b/Inputs/EventGenerator/alpha.source index 6cb81d2825f67c2a17088cd685af319ea22cd7a2..ea19e4bc82f69f9e898e26f1beadd9a70658a35b 100644 --- a/Inputs/EventGenerator/alpha.source +++ b/Inputs/EventGenerator/alpha.source @@ -4,13 +4,13 @@ % Energy are given in MeV , Position in mm % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Isotropic - EnergyLow= 30 MeV - EnergyHigh= 30 MeV + EnergyLow= 5.5 MeV + EnergyHigh= 5.5 MeV HalfOpenAngleMin= 0 deg - HalfOpenAngleMax= 180 deg + HalfOpenAngleMax= 90 deg x0= 0 mm y0= 0 mm - z0= 0 mm + z0= -110 mm Particle= alpha ExcitationEnergy= 0 MeV diff --git a/Inputs/EventGenerator/proton.source b/Inputs/EventGenerator/proton.source index 59809fe9e40efbcfc60679c7961af628121436ba..7b62f807a43297d6ba4566711765a7c559e00578 100644 --- a/Inputs/EventGenerator/proton.source +++ b/Inputs/EventGenerator/proton.source @@ -4,16 +4,16 @@ % Energy are given in MeV , Position in mm % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Isotropic - EnergyLow= 0 - EnergyHigh= 200 + EnergyLow= 30 + EnergyHigh= 80 HalfOpenAngleMin= 0 - HalfOpenAngleMax= 180 + HalfOpenAngleMax= 1 x0= 0 y0= 0 z0= 0 SigmaX= 0 SigmaY= 0 - Multiplicity= 4 + Multiplicity= 1 Particle= proton %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Supported particle type: proton, neutron, deuton, triton, He3 , alpha diff --git a/NPLib/Core/NPDetectorManager.cxx b/NPLib/Core/NPDetectorManager.cxx index 55609deff2f53937ebaaed1d38cd2e6ac7d499de..e07f8b5427605ecaddd8720f581b8c36f06ba415 100644 --- a/NPLib/Core/NPDetectorManager.cxx +++ b/NPLib/Core/NPDetectorManager.cxx @@ -191,22 +191,28 @@ void NPL::DetectorManager::ReadConfigurationFile(string Path) { void NPL::DetectorManager::BuildPhysicalEvent(){ #if __cplusplus > 199711L // add new job +//cout << "TEST0a" << endl; map<string,VDetector*>::iterator it; unsigned int i = 0; for (it = m_Detector.begin(); it != m_Detector.end(); ++it) { +//cout << "TEST0" << endl; m_Ready[i++]=true; } - +//cout << "TEST1" << endl; { // aquire the sub thread lock std::unique_lock<std::mutex> lk(m_ThreadMtx); } m_CV.notify_all(); +//cout << "TEST2" << endl; while(!IsDone()){ - // this_thread::yield(); +//cout << "TEST2a" << endl; + //this_thread::yield(); } +//cout << "TEST2b" << endl; #else +//cout << "TEST3" << endl; map<string,VDetector*>::iterator it; for (it = m_Detector.begin(); it != m_Detector.end(); ++it) { (it->second->*m_ClearEventPhysicsPtr)(); @@ -375,6 +381,7 @@ void NPL::DetectorManager::StartThread(NPL::VDetector* det,unsigned int id){ vector<bool>::iterator it = m_Ready.begin()+id; while(true){ { // Aquire the lock +////cout << "WWWW" << endl; std::unique_lock<std::mutex> lk(m_ThreadMtx); // wait for work to be given while(!m_Ready[id]){ @@ -409,9 +416,16 @@ void NPL::DetectorManager::StopThread(){ } //////////////////////////////////////////////////////////////////////////////// bool NPL::DetectorManager::IsDone(){ +int ijk=0; +//cout << m_Ready.size() << " !" << endl; for(vector<bool>::iterator it = m_Ready.begin() ; it!=m_Ready.end() ; it++){ if((*it)) +{ +ijk++; +//cout << *it << endl; +//cout << ijk << endl; return false; +} } return true; } diff --git a/NPLib/Core/NPGlobalSystemOfUnits.h b/NPLib/Core/NPGlobalSystemOfUnits.h index 3e53094b3b924d2d22d61250101af91455b42b08..3eed6c7b82185adccea7175599234e7c9a078f54 100644 --- a/NPLib/Core/NPGlobalSystemOfUnits.h +++ b/NPLib/Core/NPGlobalSystemOfUnits.h @@ -106,6 +106,7 @@ using NPUNITS::watt; using NPUNITS::newton; using NPUNITS::hep_pascal; using NPUNITS::bar; +using NPUNITS::torr; using NPUNITS::atmosphere; using NPUNITS::ampere; using NPUNITS::milliampere; diff --git a/NPLib/Core/NPSystemOfUnits.h b/NPLib/Core/NPSystemOfUnits.h index d6513b20982f7c315e2de09104b38929edf6a0f0..11feef095077ee3418546397fee2a5ad41aa54c0 100644 --- a/NPLib/Core/NPSystemOfUnits.h +++ b/NPLib/Core/NPSystemOfUnits.h @@ -25,246 +25,247 @@ #define NP_SYSTEM_OF_UNITS_H namespace NPUNITS { - - // - // Length [L] - // - static const double millimeter = 1.; - static const double millimeter2 = millimeter*millimeter; - static const double millimeter3 = millimeter*millimeter*millimeter; - - static const double centimeter = 10.*millimeter; - static const double centimeter2 = centimeter*centimeter; - static const double centimeter3 = centimeter*centimeter*centimeter; - - static const double meter = 1000.*millimeter; - static const double meter2 = meter*meter; - static const double meter3 = meter*meter*meter; - - static const double kilometer = 1000.*meter; - static const double kilometer2 = kilometer*kilometer; - static const double kilometer3 = kilometer*kilometer*kilometer; - - static const double parsec = 3.0856775807e+16*meter; - - static const double micrometer = 1.e-6 *meter; - static const double nanometer = 1.e-9 *meter; - static const double angstrom = 1.e-10*meter; - static const double fermi = 1.e-15*meter; - - static const double barn = 1.e-28*meter2; - static const double millibarn = 1.e-3 *barn; - static const double microbarn = 1.e-6 *barn; - static const double nanobarn = 1.e-9 *barn; - static const double picobarn = 1.e-12*barn; - - // symbols - static const double nm = nanometer; - static const double um = micrometer; - - static const double mm = millimeter; - static const double mm2 = millimeter2; - static const double mm3 = millimeter3; - - static const double cm = centimeter; - static const double cm2 = centimeter2; - static const double cm3 = centimeter3; - - static const double m = meter; - static const double m2 = meter2; - static const double m3 = meter3; - - static const double km = kilometer; - static const double km2 = kilometer2; - static const double km3 = kilometer3; - - static const double pc = parsec; - - // - // Angle - // - static const double radian = 1.; - static const double milliradian = 1.e-3*radian; - static const double degree = (3.14159265358979323846/180.0)*radian; - - static const double steradian = 1.; - - // symbols - static const double rad = radian; - static const double mrad = milliradian; - static const double sr = steradian; - static const double deg = degree; - - // - // Time [T] - // - static const double nanosecond = 1.; - static const double second = 1.e+9 *nanosecond; - static const double millisecond = 1.e-3 *second; - static const double microsecond = 1.e-6 *second; - static const double picosecond = 1.e-12*second; - - static const double hertz = 1./second; - static const double kilohertz = 1.e+3*hertz; - static const double megahertz = 1.e+6*hertz; - - // symbols - static const double ns = nanosecond; - static const double us = microsecond; - static const double s = second; - static const double ms = millisecond; - - // - // Electric charge [Q] - // - static const double eplus = 1. ;// positron charge - static const double e_SI = 1.602176487e-19;// positron charge in coulomb - static const double coulomb = eplus/e_SI;// coulomb = 6.24150 e+18 * eplus - - // - // Energy [E] - // - static const double megaelectronvolt = 1. ; - static const double electronvolt = 1.e-6*megaelectronvolt; - static const double kiloelectronvolt = 1.e-3*megaelectronvolt; - static const double gigaelectronvolt = 1.e+3*megaelectronvolt; - static const double teraelectronvolt = 1.e+6*megaelectronvolt; - static const double petaelectronvolt = 1.e+9*megaelectronvolt; - - static const double joule = electronvolt/e_SI;// joule = 6.24150 e+12 * MeV - - // symbols - static const double MeV = megaelectronvolt; - static const double eV = electronvolt; - static const double keV = kiloelectronvolt; - static const double GeV = gigaelectronvolt; - static const double TeV = teraelectronvolt; - static const double PeV = petaelectronvolt; - - // - // Mass [E][T^2][L^-2] - // - static const double kilogram = joule*second*second/(meter*meter); - static const double gram = 1.e-3*kilogram; - static const double milligram = 1.e-3*gram; - - // symbols - static const double kg = kilogram; - static const double g = gram; - static const double mg = milligram; - - // - // Power [E][T^-1] - // - static const double watt = joule/second;// watt = 6.24150 e+3 * MeV/ns - - // - // Force [E][L^-1] - // - static const double newton = joule/meter;// newton = 6.24150 e+9 * MeV/mm - - // - // Pressure [E][L^-3] - // -//#define pascal hep_pascal // a trick to avoid warnings - static const double hep_pascal = newton/m2; // pascal = 6.24150 e+3 * MeV/mm3 - static const double pascal = newton/m2; // pascal = 6.24150 e+3 * MeV/mm3 - static const double bar = 100000*pascal; // bar = 6.24150 e+8 * MeV/mm3 - static const double atmosphere = 101325*pascal; // atm = 6.32420 e+8 * MeV/mm3 - - // - // Electric current [Q][T^-1] - // - static const double ampere = coulomb/second; // ampere = 6.24150 e+9 * eplus/ns - static const double milliampere = 1.e-3*ampere; - static const double microampere = 1.e-6*ampere; - static const double nanoampere = 1.e-9*ampere; - - // - // Electric potential [E][Q^-1] - // - static const double megavolt = megaelectronvolt/eplus; - static const double kilovolt = 1.e-3*megavolt; - static const double volt = 1.e-6*megavolt; - - // - // Electric resistance [E][T][Q^-2] - // - static const double ohm = volt/ampere;// ohm = 1.60217e-16*(MeV/eplus)/(eplus/ns) - - // - // Electric capacitance [Q^2][E^-1] - // - static const double farad = coulomb/volt;// farad = 6.24150e+24 * eplus/Megavolt - static const double millifarad = 1.e-3*farad; - static const double microfarad = 1.e-6*farad; - static const double nanofarad = 1.e-9*farad; - static const double picofarad = 1.e-12*farad; - - // - // Magnetic Flux [T][E][Q^-1] - // - static const double weber = volt*second;// weber = 1000*megavolt*ns - - // - // Magnetic Field [T][E][Q^-1][L^-2] - // - static const double tesla = volt*second/meter2;// tesla =0.001*megavolt*ns/mm2 - - static const double gauss = 1.e-4*tesla; - static const double kilogauss = 1.e-1*tesla; - - // - // Inductance [T^2][E][Q^-2] - // - static const double henry = weber/ampere;// henry = 1.60217e-7*MeV*(ns/eplus)**2 - - // - // Temperature - // - static const double kelvin = 1.; - - // - // Amount of substance - // - static const double mole = 1.; - - // - // Activity [T^-1] - // - static const double becquerel = 1./second ; - static const double curie = 3.7e+10 * becquerel; - - // - // Absorbed dose [L^2][T^-2] - // - static const double gray = joule/kilogram ; - static const double kilogray = 1.e+3*gray; - static const double milligray = 1.e-3*gray; - static const double microgray = 1.e-6*gray; - - // - // Luminous intensity [I] - // - static const double candela = 1.; - - // - // Luminous flux [I] - // - static const double lumen = candela*steradian; - - // - // Illuminance [I][L^-2] - // - static const double lux = lumen/meter2; - - // - // Miscellaneous - // - static const double perCent = 0.01 ; - static const double perThousand = 0.001; - static const double perMillion = 0.000001; - + + // + // Length [L] + // + static const double millimeter = 1.; + static const double millimeter2 = millimeter*millimeter; + static const double millimeter3 = millimeter*millimeter*millimeter; + + static const double centimeter = 10.*millimeter; + static const double centimeter2 = centimeter*centimeter; + static const double centimeter3 = centimeter*centimeter*centimeter; + + static const double meter = 1000.*millimeter; + static const double meter2 = meter*meter; + static const double meter3 = meter*meter*meter; + + static const double kilometer = 1000.*meter; + static const double kilometer2 = kilometer*kilometer; + static const double kilometer3 = kilometer*kilometer*kilometer; + + static const double parsec = 3.0856775807e+16*meter; + + static const double micrometer = 1.e-6 *meter; + static const double nanometer = 1.e-9 *meter; + static const double angstrom = 1.e-10*meter; + static const double fermi = 1.e-15*meter; + + static const double barn = 1.e-28*meter2; + static const double millibarn = 1.e-3 *barn; + static const double microbarn = 1.e-6 *barn; + static const double nanobarn = 1.e-9 *barn; + static const double picobarn = 1.e-12*barn; + + // symbols + static const double nm = nanometer; + static const double um = micrometer; + + static const double mm = millimeter; + static const double mm2 = millimeter2; + static const double mm3 = millimeter3; + + static const double cm = centimeter; + static const double cm2 = centimeter2; + static const double cm3 = centimeter3; + + static const double m = meter; + static const double m2 = meter2; + static const double m3 = meter3; + + static const double km = kilometer; + static const double km2 = kilometer2; + static const double km3 = kilometer3; + + static const double pc = parsec; + + // + // Angle + // + static const double radian = 1.; + static const double milliradian = 1.e-3*radian; + static const double degree = (3.14159265358979323846/180.0)*radian; + + static const double steradian = 1.; + + // symbols + static const double rad = radian; + static const double mrad = milliradian; + static const double sr = steradian; + static const double deg = degree; + + // + // Time [T] + // + static const double nanosecond = 1.; + static const double second = 1.e+9 *nanosecond; + static const double millisecond = 1.e-3 *second; + static const double microsecond = 1.e-6 *second; + static const double picosecond = 1.e-12*second; + + static const double hertz = 1./second; + static const double kilohertz = 1.e+3*hertz; + static const double megahertz = 1.e+6*hertz; + + // symbols + static const double ns = nanosecond; + static const double us = microsecond; + static const double s = second; + static const double ms = millisecond; + + // + // Electric charge [Q] + // + static const double eplus = 1. ;// positron charge + static const double e_SI = 1.602176487e-19;// positron charge in coulomb + static const double coulomb = eplus/e_SI;// coulomb = 6.24150 e+18 * eplus + + // + // Energy [E] + // + static const double megaelectronvolt = 1. ; + static const double electronvolt = 1.e-6*megaelectronvolt; + static const double kiloelectronvolt = 1.e-3*megaelectronvolt; + static const double gigaelectronvolt = 1.e+3*megaelectronvolt; + static const double teraelectronvolt = 1.e+6*megaelectronvolt; + static const double petaelectronvolt = 1.e+9*megaelectronvolt; + + static const double joule = electronvolt/e_SI;// joule = 6.24150 e+12 * MeV + + // symbols + static const double MeV = megaelectronvolt; + static const double eV = electronvolt; + static const double keV = kiloelectronvolt; + static const double GeV = gigaelectronvolt; + static const double TeV = teraelectronvolt; + static const double PeV = petaelectronvolt; + + // + // Mass [E][T^2][L^-2] + // + static const double kilogram = joule*second*second/(meter*meter); + static const double gram = 1.e-3*kilogram; + static const double milligram = 1.e-3*gram; + + // symbols + static const double kg = kilogram; + static const double g = gram; + static const double mg = milligram; + + // + // Power [E][T^-1] + // + static const double watt = joule/second;// watt = 6.24150 e+3 * MeV/ns + + // + // Force [E][L^-1] + // + static const double newton = joule/meter;// newton = 6.24150 e+9 * MeV/mm + + // + // Pressure [E][L^-3] + // + //#define pascal hep_pascal // a trick to avoid warnings + static const double hep_pascal = newton/m2; // pascal = 6.24150 e+3 * MeV/mm3 + static const double pascal = newton/m2; // pascal = 6.24150 e+3 * MeV/mm3 + static const double bar = 100000*pascal; // bar = 6.24150 e+8 * MeV/mm3 + static const double atmosphere = 101325*pascal; // atm = 6.32420 e+8 * MeV/mm3 + static const double torr = 0.00133322*bar; + + // + // Electric current [Q][T^-1] + // + static const double ampere = coulomb/second; // ampere = 6.24150 e+9 * eplus/ns + static const double milliampere = 1.e-3*ampere; + static const double microampere = 1.e-6*ampere; + static const double nanoampere = 1.e-9*ampere; + + // + // Electric potential [E][Q^-1] + // + static const double megavolt = megaelectronvolt/eplus; + static const double kilovolt = 1.e-3*megavolt; + static const double volt = 1.e-6*megavolt; + + // + // Electric resistance [E][T][Q^-2] + // + static const double ohm = volt/ampere;// ohm = 1.60217e-16*(MeV/eplus)/(eplus/ns) + + // + // Electric capacitance [Q^2][E^-1] + // + static const double farad = coulomb/volt;// farad = 6.24150e+24 * eplus/Megavolt + static const double millifarad = 1.e-3*farad; + static const double microfarad = 1.e-6*farad; + static const double nanofarad = 1.e-9*farad; + static const double picofarad = 1.e-12*farad; + + // + // Magnetic Flux [T][E][Q^-1] + // + static const double weber = volt*second;// weber = 1000*megavolt*ns + + // + // Magnetic Field [T][E][Q^-1][L^-2] + // + static const double tesla = volt*second/meter2;// tesla =0.001*megavolt*ns/mm2 + + static const double gauss = 1.e-4*tesla; + static const double kilogauss = 1.e-1*tesla; + + // + // Inductance [T^2][E][Q^-2] + // + static const double henry = weber/ampere;// henry = 1.60217e-7*MeV*(ns/eplus)**2 + + // + // Temperature + // + static const double kelvin = 1.; + + // + // Amount of substance + // + static const double mole = 1.; + + // + // Activity [T^-1] + // + static const double becquerel = 1./second ; + static const double curie = 3.7e+10 * becquerel; + + // + // Absorbed dose [L^2][T^-2] + // + static const double gray = joule/kilogram ; + static const double kilogray = 1.e+3*gray; + static const double milligray = 1.e-3*gray; + static const double microgray = 1.e-6*gray; + + // + // Luminous intensity [I] + // + static const double candela = 1.; + + // + // Luminous flux [I] + // + static const double lumen = candela*steradian; + + // + // Illuminance [I][L^-2] + // + static const double lux = lumen/meter2; + + // + // Miscellaneous + // + static const double perCent = 0.01 ; + static const double perThousand = 0.001; + static const double perMillion = 0.000001; + } // namespace NPUNITS #endif diff --git a/NPLib/Detectors/Actar/CMakeLists.txt b/NPLib/Detectors/Actar/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..a1b027f6992044326fbeefd56d43a1e2c8150d56 --- /dev/null +++ b/NPLib/Detectors/Actar/CMakeLists.txt @@ -0,0 +1,6 @@ +add_custom_command(OUTPUT TActarPhysicsDict.cxx COMMAND ../../scripts/build_dict.sh TActarPhysics.h TActarPhysicsDict.cxx TActarPhysics.rootmap libNPActar.dylib DEPENDS TActarPhysics.h) +add_custom_command(OUTPUT TActarDataDict.cxx COMMAND ../../scripts/build_dict.sh TActarData.h TActarDataDict.cxx TActarData.rootmap libNPActar.dylib DEPENDS TActarData.h) +add_library(NPActar SHARED TActarSpectra.cxx TActarData.cxx TActarPhysics.cxx TActarDataDict.cxx TActarPhysicsDict.cxx ) +target_link_libraries(NPActar ${ROOT_LIBRARIES} NPCore) +install(FILES TActarData.h TActarPhysics.h TActarSpectra.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY}) + diff --git a/NPLib/Detectors/Actar/TActarData.cxx b/NPLib/Detectors/Actar/TActarData.cxx new file mode 100644 index 0000000000000000000000000000000000000000..a9462d0ee2a68c4ca7dc012ee52835fb712e9d1d --- /dev/null +++ b/NPLib/Detectors/Actar/TActarData.cxx @@ -0,0 +1,88 @@ +/***************************************************************************** + * Copyright (C) 2009-2017 this file is part of the NPTool Project * + * * + * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * + * For the list of contributors see $NPTOOL/Licence/Contributors * + *****************************************************************************/ + +/***************************************************************************** + * Original Author: Pierre Morfouace contact address: morfouac@nscl.msu.edu * + * * + * Creation Date : September 2017 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class hold Actar Raw data * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ +#include "TActarData.h" + +#include <iostream> +#include <fstream> +#include <sstream> +#include <string> +using namespace std; + +ClassImp(TActarData) + + +////////////////////////////////////////////////////////////////////// +TActarData::TActarData() { +} + + + +////////////////////////////////////////////////////////////////////// +TActarData::~TActarData() { +} + + + +////////////////////////////////////////////////////////////////////// +void TActarData::Clear() { + // Charge + fActar_PadNumber.clear(); + fActar_PadRow.clear(); + fActar_PadColumn.clear(); + fActar_PadCharge.clear(); + fActar_PadTime.clear(); + fActar_dig_Charge.clear(); + fActar_dig_Time.clear(); + + fSilicon_Energy.clear(); + fSilicon_Time.clear(); + fSilicon_DetectorNumber.clear(); + + fCsI_Energy.clear(); + fCsI_CrystalNumber.clear(); +} + +//////////////////////////////////////////////////////////////////////////////// +TGraph* TActarData::GetChargeAsGraph(){ + TGraph* res = new TGraph (fActar_dig_Charge.size(),&fActar_dig_Time[0],&fActar_dig_Charge[0]); + res->Sort(); + return res; + +} + + +////////////////////////////////////////////////////////////////////// +void TActarData::Dump() const { + // This method is very useful for debuging and worth the dev. + cout << "XXXXXXXXXXXXXXXXXXXXXXXX New Event [TActarData::Dump()] XXXXXXXXXXXXXXXXX" << endl; + + // Charge + size_t mysize = fActar_PadNumber.size(); + cout << "Actar_Mult: " << mysize << endl; + + for (size_t i = 0 ; i < mysize ; i++){ + cout << "Pad Number: " << fActar_PadNumber[i] + << " Charge: " << fActar_PadCharge[i] + << " Time: " << fActar_PadTime[i]; + } + +} diff --git a/NPLib/Detectors/Actar/TActarData.h b/NPLib/Detectors/Actar/TActarData.h new file mode 100644 index 0000000000000000000000000000000000000000..385512cd1e87e73157fb9ea0fdcc3238cd4b4408 --- /dev/null +++ b/NPLib/Detectors/Actar/TActarData.h @@ -0,0 +1,166 @@ +#ifndef __ActarDATA__ +#define __ActarDATA__ +/***************************************************************************** + * Copyright (C) 2009-2017 this file is part of the NPTool Project * + * * + * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * + * For the list of contributors see $NPTOOL/Licence/Contributors * + *****************************************************************************/ + +/***************************************************************************** + * Original Author: Pierre Morfouace contact address: morfouac@nscl.msu.edu * + * * + * Creation Date : September 2017 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class hold Actar Raw data * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +// STL +#include <vector> +using namespace std; + +// ROOT +#include "TObject.h" +#include "TGraph.h" + +class TActarData : public TObject { + ////////////////////////////////////////////////////////////// + // data members are hold into vectors in order + // to allow multiplicity treatment + private: + // Charge + vector<int> fActar_PadNumber; + vector<int> fActar_PadTime; + vector<int> fActar_PadRow; + vector<int> fActar_PadColumn; + vector<double> fActar_PadCharge; + + vector<double> fActar_dig_Time; + vector<double> fActar_dig_Charge; + + vector<double> fSilicon_Energy; + vector<int> fSilicon_DetectorNumber; + vector<double> fSilicon_Time; + + vector<double> fCsI_Energy; + vector<int> fCsI_CrystalNumber; + + + ////////////////////////////////////////////////////////////// + // Constructor and destructor + public: + TActarData(); + ~TActarData(); + + + ////////////////////////////////////////////////////////////// + // Inherited from TObject and overriden to avoid warnings + public: + void Clear(); + void Clear(const Option_t*) {}; + void Dump() const; + + + ////////////////////////////////////////////////////////////// + // Getters and Setters + // Prefer inline declaration to avoid unnecessary called of + // frequently used methods + // add //! to avoid ROOT creating dictionnary for the methods + public: + ////////////////////// SETTERS //////////////////////// + // Charge + inline void SetCharge(const Double_t& Charge){ + fActar_PadCharge.push_back(Charge); + };//! + // Pad + inline void SetPadNumber(const UShort_t& PadNbr){ + fActar_PadNumber.push_back(PadNbr); + };//! + // Row + inline void SetRowNumber(const UShort_t& RowNbr){ + fActar_PadRow.push_back(RowNbr); + } + // Column + inline void SetColumnNumber(const UShort_t& ColumnNbr){ + fActar_PadColumn.push_back(ColumnNbr); + } + // Time + inline void SetTime(const Double_t& Time) { + fActar_PadTime.push_back(Time); + };//! + + // Digitizer + void AddEnergyPoint(double& E,double& T){ + fActar_dig_Charge.push_back(E); + fActar_dig_Time.push_back(T); + } + + //Silicon + inline void SetSiliconEnergy(const Double_t& Energy){ + fSilicon_Energy.push_back(Energy); + } + inline void SetSiliconDetectorNumber(const Int_t& Det){ + fSilicon_DetectorNumber.push_back(Det); + } + inline void SetSiliconTime(const Double_t& Time){ + fSilicon_Time.push_back(Time); + } + + //CsI + inline void SetCsIEnergy(const Double_t& Energy){ + fCsI_Energy.push_back(Energy); + } + inline void SetCsICrystalNumber(const Int_t& Det){ + fCsI_CrystalNumber.push_back(Det); + } + + + ////////////////////// GETTERS //////////////////////// + // Charge + inline UShort_t GetMultCharge() const + {return fActar_PadNumber.size();} + inline UShort_t Get_PadNumber(const unsigned int &i) const + {return fActar_PadNumber[i];}//! + inline UShort_t Get_PadRow(const unsigned int &i) const + {return fActar_PadRow[i];}//! + inline UShort_t Get_PadColumn(const unsigned int &i) const + {return fActar_PadColumn[i];}//! + inline Double_t Get_Charge(const unsigned int &i) const + {return fActar_PadCharge[i];}//! + + // Time + inline Double_t Get_Time(const unsigned int &i) const + {return fActar_PadTime[i];}//! + + inline vector<double> Get_dig_Charge() {return fActar_dig_Charge;} + TGraph* GetChargeAsGraph(); + + // Silicon + inline Double_t Get_SiliconEnergy(const unsigned int &i) const + {return fSilicon_Energy[i];}//! + inline Double_t Get_SiliconTime(const unsigned int &i) const + {return fSilicon_Time[i];}//! + inline Int_t Get_SiliconDetectorNumber(const unsigned int &i) const + {return fSilicon_DetectorNumber[i];}//! + + // CsI + inline Double_t Get_CsIEnergy(const unsigned int &i) const + {return fCsI_Energy[i];}//! + inline Int_t Get_CsICrystalNumber(const unsigned int &i) const + {return fCsI_CrystalNumber[i];}//! + + + + ////////////////////////////////////////////////////////////// + // Required for ROOT dictionnary + ClassDef(TActarData,1) // ActarData structure +}; + +#endif diff --git a/NPLib/Detectors/Actar/TActarPhysics.cxx b/NPLib/Detectors/Actar/TActarPhysics.cxx new file mode 100644 index 0000000000000000000000000000000000000000..c1b79c0d872d73695783ca9e8e88e67e239148ac --- /dev/null +++ b/NPLib/Detectors/Actar/TActarPhysics.cxx @@ -0,0 +1,374 @@ +/***************************************************************************** + * Copyright (C) 2009-2017 this file is part of the NPTool Project * + * * + * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * + * For the list of contributors see $NPTOOL/Licence/Contributors * + *****************************************************************************/ + +/***************************************************************************** + * Original Author: Pierre Morfouace contact address: morfouac@nscl.msu.edu * + * * + * Creation Date : September 2017 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class hold Actar Treated data * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +#include "TActarPhysics.h" + +// STL +#include <sstream> +#include <iostream> +#include <cmath> +#include <stdlib.h> +#include <limits> +using namespace std; + +// NPLNPSystemOfUnits.h +#include "RootInput.h" +#include "RootOutput.h" +#include "NPDetectorFactory.h" +#include "NPOptionManager.h" +#include "NPSystemOfUnits.h" + +// ROOT +#include "TChain.h" +#include "TH2F.h" +#include "TCanvas.h" + +ClassImp(TActarPhysics) + + +/////////////////////////////////////////////////////////////////////////// +TActarPhysics::TActarPhysics() + : m_EventData(new TActarData), + m_PreTreatedData(new TActarData), + m_EventPhysics(this), + m_Spectra(0), + m_E_RAW_Threshold(0), // adc channels + m_E_Threshold(0), // MeV + m_NumberOfDetectors(0) { +} + +/////////////////////////////////////////////////////////////////////////// +/// A usefull method to bundle all operation to add a detector +void TActarPhysics::AddDetector(TVector3 , string ){ + // In That simple case nothing is done + // Typically for more complex detector one would calculate the relevant + // positions (stripped silicon) or angles (gamma array) + m_NumberOfDetectors++; +} + +/////////////////////////////////////////////////////////////////////////// +void TActarPhysics::AddDetector(double R, double Theta, double Phi, string shape){ + // Compute the TVector3 corresponding + TVector3 Pos(R*sin(Theta)*cos(Phi),R*sin(Theta)*sin(Phi),R*cos(Theta)); + // Call the cartesian method + AddDetector(Pos,shape); +} + +/////////////////////////////////////////////////////////////////////////// +void TActarPhysics::BuildSimplePhysicalEvent() { + BuildPhysicalEvent(); +} + + + +/////////////////////////////////////////////////////////////////////////// +void TActarPhysics::BuildPhysicalEvent() { + // apply thresholds and calibration + + PreTreat(); + + // match Charge and time together + unsigned int mysizeE = m_PreTreatedData->GetMultCharge(); + + for (UShort_t e = 0; e < mysizeE ; e++) { + PadNumber.push_back(m_PreTreatedData->Get_PadNumber(e)); + PadRow.push_back(m_PreTreatedData->Get_PadRow(e)); + PadColumn.push_back(m_PreTreatedData->Get_PadColumn(e)); + PadCharge.push_back(m_PreTreatedData->Get_Charge(e)); + PadTime.push_back(m_PreTreatedData->Get_Time(e)); + } + + HoughTransform(PadRow, PadColumn, PadTime); +} + +/////////////////////////////////////////////////////////////////////////// +void TActarPhysics::PreTreat() { + // This method typically applies thresholds and calibrations + // Might test for disabled channels for more complex detector + + // clear pre-treated object + ClearPreTreatedData(); + + // instantiate CalibrationManager + static CalibrationManager* Cal = CalibrationManager::getInstance(); + + // Charge + unsigned int mysize = m_EventData->GetMultCharge(); + + for (UShort_t i = 0; i < mysize ; ++i) { + if (m_EventData->Get_Charge(i) > m_E_RAW_Threshold) { + //Double_t Charge = Cal->ApplyCalibration("Actar/CHARGE"+NPL::itoa(m_EventData->Get_PadNumber(i)),m_EventData->Get_Charge(i)); + //Double_t Time= Cal->ApplyCalibration("Actar/TIME"+NPL::itoa(m_EventData->Get_PadNumber(i)),m_EventData->Get_Time(i)); + if (m_EventData->Get_Charge(i) > m_E_Threshold) { + m_PreTreatedData->SetCharge(m_EventData->Get_Charge(i)); + m_PreTreatedData->SetPadNumber(m_EventData->Get_PadNumber(i)); + m_PreTreatedData->SetRowNumber(m_EventData->Get_PadRow(i)); + m_PreTreatedData->SetColumnNumber(m_EventData->Get_PadColumn(i)); + m_PreTreatedData->SetTime(m_EventData->Get_Time(i)); + } + } + } +} + +/////////////////////////////////////////////////////////////////////////// +void TActarPhysics::HoughTransform(vector<int> v1, vector<int> v2, vector<int> v3){ + for(unsigned int i=0; i<v1.size(); i++){ + for(int itheta=0; itheta<180; itheta++){ + HoughAngle.push_back(itheta); + double tmp= v1[i]*cos(itheta*NPUNITS::deg)+v2[i]*sin(itheta*NPUNITS::deg); + HoughRadius.push_back(tmp); + } + } +} + +/////////////////////////////////////////////////////////////////////////// +void TActarPhysics::ReadAnalysisConfig() { + bool ReadingStatus = false; + + // path to file + string FileName = "./configs/ConfigActar.dat"; + + // open analysis config file + ifstream AnalysisConfigFile; + AnalysisConfigFile.open(FileName.c_str()); + + if (!AnalysisConfigFile.is_open()) { + cout << " No ConfigActar.dat found: Default parameter loaded for Analayis " << FileName << endl; + return; + } + cout << " Loading user parameter for Analysis from ConfigActar.dat " << endl; + + // Save it in a TAsciiFile + TAsciiFile* asciiConfig = RootOutput::getInstance()->GetAsciiFileAnalysisConfig(); + asciiConfig->AppendLine("%%% ConfigActar.dat %%%"); + asciiConfig->Append(FileName.c_str()); + asciiConfig->AppendLine(""); + // read analysis config file + string LineBuffer,DataBuffer,whatToDo; + while (!AnalysisConfigFile.eof()) { + // Pick-up next line + getline(AnalysisConfigFile, LineBuffer); + + // search for "header" + string name = "ConfigActar"; + if (LineBuffer.compare(0, name.length(), name) == 0) + ReadingStatus = true; + + // loop on tokens and data + while (ReadingStatus ) { + whatToDo=""; + AnalysisConfigFile >> whatToDo; + + // Search for comment symbol (%) + if (whatToDo.compare(0, 1, "%") == 0) { + AnalysisConfigFile.ignore(numeric_limits<streamsize>::max(), '\n' ); + } + + else if (whatToDo=="E_RAW_THRESHOLD") { + AnalysisConfigFile >> DataBuffer; + m_E_RAW_Threshold = atof(DataBuffer.c_str()); + cout << whatToDo << " " << m_E_RAW_Threshold << endl; + } + + else if (whatToDo=="E_THRESHOLD") { + AnalysisConfigFile >> DataBuffer; + m_E_Threshold = atof(DataBuffer.c_str()); + cout << whatToDo << " " << m_E_Threshold << endl; + } + + else { + ReadingStatus = false; + } + } + } +} + + + +/////////////////////////////////////////////////////////////////////////// +void TActarPhysics::Clear() { + PadNumber.clear(); + PadRow.clear(); + PadColumn.clear(); + PadCharge.clear(); + PadTime.clear(); + + HoughRadius.clear(); + HoughAngle.clear(); +} + + + +/////////////////////////////////////////////////////////////////////////// +void TActarPhysics::ReadConfiguration(NPL::InputParser parser) { + vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Actar"); + if(NPOptionManager::getInstance()->GetVerboseLevel()) + cout << "//// " << blocks.size() << " detectors found " << endl; + + vector<string> cart = {"POS","Shape"}; + vector<string> sphe = {"R","Theta","Phi","Shape"}; + + for(unsigned int i = 0 ; i < blocks.size() ; i++){ + if(blocks[i]->HasTokenList(cart)){ + if(NPOptionManager::getInstance()->GetVerboseLevel()) + cout << endl << "//// Actar " << i+1 << endl; + + TVector3 Pos = blocks[i]->GetTVector3("POS","mm"); + string Shape = blocks[i]->GetString("Shape"); + AddDetector(Pos,Shape); + } + else if(blocks[i]->HasTokenList(sphe)){ + if(NPOptionManager::getInstance()->GetVerboseLevel()) + cout << endl << "//// Actar " << i+1 << endl; + double R = blocks[i]->GetDouble("R","mm"); + double Theta = blocks[i]->GetDouble("Theta","deg"); + double Phi = blocks[i]->GetDouble("Phi","deg"); + string Shape = blocks[i]->GetString("Shape"); + AddDetector(R,Theta,Phi,Shape); + } + else{ + cout << "ERROR: check your input file formatting " << endl; + exit(1); + } + } +} + +/////////////////////////////////////////////////////////////////////////// +void TActarPhysics::InitSpectra() { + m_Spectra = new TActarSpectra(m_NumberOfDetectors); +} + + + +/////////////////////////////////////////////////////////////////////////// +void TActarPhysics::FillSpectra() { + m_Spectra -> FillRawSpectra(m_EventData); + m_Spectra -> FillPreTreatedSpectra(m_PreTreatedData); + m_Spectra -> FillPhysicsSpectra(m_EventPhysics); +} + + + +/////////////////////////////////////////////////////////////////////////// +void TActarPhysics::CheckSpectra() { + m_Spectra->CheckSpectra(); +} + + + +/////////////////////////////////////////////////////////////////////////// +void TActarPhysics::ClearSpectra() { + // To be done +} + + + +/////////////////////////////////////////////////////////////////////////// +map< string , TH1*> TActarPhysics::GetSpectra() { + if(m_Spectra) + return m_Spectra->GetMapHisto(); + else{ + map< string , TH1*> empty; + return empty; + } +} + + + +//////////////////////////////////////////////////////////////////////////////// +vector<TCanvas*> TActarPhysics::GetCanvas() { + //if(m_Spectra) + //return m_Spectra->GetCanvas(); + //else{ + //vector<TCanvas*> empty; + //return empty; + //} +} + + + +/////////////////////////////////////////////////////////////////////////// +void TActarPhysics::WriteSpectra() { + m_Spectra->WriteSpectra(); +} + + + +/////////////////////////////////////////////////////////////////////////// +void TActarPhysics::AddParameterToCalibrationManager() { + CalibrationManager* Cal = CalibrationManager::getInstance(); + for (int i = 0; i < m_NumberOfDetectors; ++i) { + Cal->AddParameter("Actar", "D"+ NPL::itoa(i+1)+"_CHARGE","Actar_D"+ NPL::itoa(i+1)+"_CHARGE"); + Cal->AddParameter("Actar", "D"+ NPL::itoa(i+1)+"_TIME","Actar_D"+ NPL::itoa(i+1)+"_TIME"); + } +} + + + +/////////////////////////////////////////////////////////////////////////// +void TActarPhysics::InitializeRootInputRaw() { + TChain* inputChain = RootInput::getInstance()->GetChain(); + inputChain->SetBranchStatus("Actar", true ); + inputChain->SetBranchAddress("Actar", &m_EventData ); +} + + + +/////////////////////////////////////////////////////////////////////////// +void TActarPhysics::InitializeRootInputPhysics() { + TChain* inputChain = RootInput::getInstance()->GetChain(); + inputChain->SetBranchAddress("Actar", &m_EventPhysics); +} + + + +/////////////////////////////////////////////////////////////////////////// +void TActarPhysics::InitializeRootOutput() { + TTree* outputTree = RootOutput::getInstance()->GetTree(); + outputTree->Branch("Actar", "TActarPhysics", &m_EventPhysics); +} + + + +//////////////////////////////////////////////////////////////////////////////// +// Construct Method to be pass to the DetectorFactory // +//////////////////////////////////////////////////////////////////////////////// +NPL::VDetector* TActarPhysics::Construct() { + return (NPL::VDetector*) new TActarPhysics(); +} + + + +//////////////////////////////////////////////////////////////////////////////// +// Registering the construct method to the factory // +//////////////////////////////////////////////////////////////////////////////// +extern "C"{ +class proxy_Actar{ + public: + proxy_Actar(){ + NPL::DetectorFactory::getInstance()->AddToken("Actar","Actar"); + NPL::DetectorFactory::getInstance()->AddDetector("Actar",TActarPhysics::Construct); + } +}; + +proxy_Actar p_Actar; +} diff --git a/NPLib/Detectors/Actar/TActarPhysics.h b/NPLib/Detectors/Actar/TActarPhysics.h new file mode 100644 index 0000000000000000000000000000000000000000..6240921b06ad69c29bafb88f883dc1e19a93ed25 --- /dev/null +++ b/NPLib/Detectors/Actar/TActarPhysics.h @@ -0,0 +1,189 @@ +#ifndef TActarPHYSICS_H +#define TActarPHYSICS_H +/***************************************************************************** + * Copyright (C) 2009-2017 this file is part of the NPTool Project * + * * + * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * + * For the list of contributors see $NPTOOL/Licence/Contributors * + *****************************************************************************/ + +/***************************************************************************** + * Original Author: Pierre Morfouace contact address: morfouac@nscl.msu.edu * + * * + * Creation Date : September 2017 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class hold Actar Treated data * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +// C++ headers +#include <vector> +#include <map> +#include <string> +using namespace std; + +// ROOT headers +#include "TObject.h" +#include "TH1.h" +#include "TCanvas.h" +#include "TVector3.h" +// NPTool headers +#include "TActarData.h" +#include "TActarSpectra.h" +#include "NPCalibrationManager.h" +#include "NPVDetector.h" +#include "NPInputParser.h" +// forward declaration +class TActarSpectra; + + + +class TActarPhysics : public TObject, public NPL::VDetector { + ////////////////////////////////////////////////////////////// + // constructor and destructor + public: + TActarPhysics(); + ~TActarPhysics() {}; + + + ////////////////////////////////////////////////////////////// + // Inherited from TObject and overriden to avoid warnings + public: + void Clear(); + void Clear(const Option_t*) {}; + + + ////////////////////////////////////////////////////////////// + // data obtained after BuildPhysicalEvent() and stored in + // output ROOT file + public: + vector<int> PadNumber; + vector<int> PadRow; + vector<int> PadColumn; + vector<int> PadTime; + vector<double> PadCharge; + + vector<double> HoughRadius; + vector<double> HoughAngle; + + /// A usefull method to bundle all operation to add a detector + void AddDetector(TVector3 POS, string shape); + void AddDetector(double R, double Theta, double Phi, string shape); + + ////////////////////////////////////////////////////////////// + // methods inherited from the VDetector ABC class + public: + // read stream from ConfigFile to pick-up detector parameters + void ReadConfiguration(NPL::InputParser); + + // add parameters to the CalibrationManger + void AddParameterToCalibrationManager(); + + // method called event by event, aiming at extracting the + // physical information from detector + void BuildPhysicalEvent(); + + // same as BuildPhysicalEvent() method but with a simpler + // treatment + void BuildSimplePhysicalEvent(); + + // same as above but for online analysis + void BuildOnlinePhysicalEvent() {BuildPhysicalEvent();}; + + // activate raw data object and branches from input TChain + // in this method mother branches (Detector) AND daughter leaves + // (fDetector_parameter) have to be activated + void InitializeRootInputRaw(); + + // activate physics data object and branches from input TChain + // in this method mother branches (Detector) AND daughter leaves + // (fDetector_parameter) have to be activated + void InitializeRootInputPhysics(); + + // create branches of output ROOT file + void InitializeRootOutput(); + + // clear the raw and physical data objects event by event + void ClearEventPhysics() {Clear();} + void ClearEventData() {m_EventData->Clear();} + + // methods related to the TActarSpectra class + // instantiate the TActarSpectra class and + // declare list of histograms + void InitSpectra(); + + // fill the spectra + void FillSpectra(); + + // used for Online mainly, sanity check for histograms and + // change their color if issues are found, for example + void CheckSpectra(); + + // used for Online only, clear all the spectra + void ClearSpectra(); + + // write spectra to ROOT output file + void WriteSpectra(); + + + ////////////////////////////////////////////////////////////// + // specific methods to Actar array + public: + // remove bad channels, calibrate the data and apply thresholds + void PreTreat(); + + // clear the pre-treated object + void ClearPreTreatedData() {m_PreTreatedData->Clear();} + + // read the user configuration file. If no file is found, load standard one + void ReadAnalysisConfig(); + + // give and external TActarData object to TActarPhysics. + // needed for online analysis for example + void SetRawDataPointer(TActarData* rawDataPointer) {m_EventData = rawDataPointer;} + + void HoughTransform(vector<int> v1, vector<int> v2, vector<int> v3); + + // objects are not written in the TTree + private: + TActarData* m_EventData; //! + TActarData* m_PreTreatedData; //! + TActarPhysics* m_EventPhysics; //! + + // getters for raw and pre-treated data object + public: + TActarData* GetRawData() const {return m_EventData;} + TActarData* GetPreTreatedData() const {return m_PreTreatedData;} + + // parameters used in the analysis + private: + // thresholds + double m_E_RAW_Threshold; //! + double m_E_Threshold; //! + + // number of detectors + private: + int m_NumberOfDetectors; //! + + // spectra class + private: + TActarSpectra* m_Spectra; // ! + + // spectra getter + public: + map<string, TH1*> GetSpectra(); + vector<TCanvas*> GetCanvas(); + + // Static constructor to be passed to the Detector Factory + public: + static NPL::VDetector* Construct(); + + ClassDef(TActarPhysics,1) // ActarPhysics structure +}; +#endif diff --git a/NPLib/Detectors/Actar/TActarSpectra.cxx b/NPLib/Detectors/Actar/TActarSpectra.cxx new file mode 100644 index 0000000000000000000000000000000000000000..f37527d1a503c3e4e94b2baa6c3ca88ca01184b2 --- /dev/null +++ b/NPLib/Detectors/Actar/TActarSpectra.cxx @@ -0,0 +1,173 @@ +/***************************************************************************** + * Copyright (C) 2009-2017 this file is part of the NPTool Project * + * * + * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * + * For the list of contributors see $NPTOOL/Licence/Contributors * + *****************************************************************************/ + +/***************************************************************************** + * Original Author: Pierre Morfouace contact address: morfouac@nscl.msu.edu * + * * + * Creation Date : September 2017 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class hold Actar Spectra * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +// class header +#include "TActarSpectra.h" + +// STL +#include <iostream> +#include <string> +using namespace std; + +// NPTool header +#include "NPOptionManager.h" + + + +//////////////////////////////////////////////////////////////////////////////// +TActarSpectra::TActarSpectra() + : fNumberOfDetectors(0) { + SetName("Actar"); +} + + + +//////////////////////////////////////////////////////////////////////////////// +TActarSpectra::TActarSpectra(unsigned int NumberOfDetectors) { + if(NPOptionManager::getInstance()->GetVerboseLevel()>0) + cout << "************************************************" << endl + << "TActarSpectra : Initalizing control spectra for " + << NumberOfDetectors << " Detectors" << endl + << "************************************************" << endl ; + SetName("Actar"); + fNumberOfDetectors = NumberOfDetectors; + + InitRawSpectra(); + InitPreTreatedSpectra(); + InitPhysicsSpectra(); +} + + + +//////////////////////////////////////////////////////////////////////////////// +TActarSpectra::~TActarSpectra() { +} + + + +//////////////////////////////////////////////////////////////////////////////// +void TActarSpectra::InitRawSpectra() { + static string name; + for (unsigned int i = 0; i < fNumberOfDetectors; i++) { // loop on number of detectors + // Charge + name = "Actar"+NPL::itoa(i+1)+"_CHARGE_RAW"; + AddHisto1D(name, name, 4096, 0, 16384, "Actar/RAW"); + // Time + name = "Actar"+NPL::itoa(i+1)+"_TIME_RAW"; + AddHisto1D(name, name, 4096, 0, 16384, "Actar/RAW"); + } // end loop on number of detectors +} + + + +//////////////////////////////////////////////////////////////////////////////// +void TActarSpectra::InitPreTreatedSpectra() { + static string name; + for (unsigned int i = 0; i < fNumberOfDetectors; i++) { // loop on number of detectors + // Energy + name = "Actar"+NPL::itoa(i+1)+"_CHARGE_CAL"; + AddHisto1D(name, name, 500, 0, 25, "Actar/CAL"); + // Time + name = "Actar"+NPL::itoa(i+1)+"_TIME_CAL"; + AddHisto1D(name, name, 500, 0, 25, "Actar/CAL"); + + + } // end loop on number of detectors +} + + + +//////////////////////////////////////////////////////////////////////////////// +void TActarSpectra::InitPhysicsSpectra() { + static string name; + // Kinematic Plot + name = "Actar_CHARGE_TIME"; + AddHisto2D(name, name, 500, 0, 500, 500, 0, 50, "Actar/PHY"); +} + + + +//////////////////////////////////////////////////////////////////////////////// +void TActarSpectra::FillRawSpectra(TActarData* RawData) { + static string name; + static string family; + + // Charge + unsigned int sizeE = RawData->GetMultCharge(); + for (unsigned int i = 0; i < sizeE; i++) { + name = "Actar"+NPL::itoa(RawData->Get_PadNumber(i))+"_CHARGE_RAW"; + family = "Actar/RAW"; + + //GetHisto(family,name) -> Fill(RawData->Get_Charge(i)); + } + + // Time + unsigned int sizeT = RawData->GetMultCharge(); + for (unsigned int i = 0; i < sizeT; i++) { + name = "Actar"+NPL::itoa(RawData->Get_PadNumber(i))+"_TIME_RAW"; + family = "Actar/RAW"; + + //GetHisto(family,name) -> Fill(RawData->Get_Time(i)); + } +} + + + +//////////////////////////////////////////////////////////////////////////////// +void TActarSpectra::FillPreTreatedSpectra(TActarData* PreTreatedData) { + static string name; + static string family; + + // Energy + unsigned int sizeE = PreTreatedData->GetMultCharge(); + for (unsigned int i = 0; i < sizeE; i++) { + name = "Actar"+NPL::itoa(PreTreatedData->Get_PadNumber(i))+"_ENERGY_CAL"; + family = "Actar/CAL"; + + //GetHisto(family,name) -> Fill(PreTreatedData->Get_Charge(i)); + } + + // Time + unsigned int sizeT = PreTreatedData->GetMultCharge(); + for (unsigned int i = 0; i < sizeT; i++) { + name = "Actar"+NPL::itoa(PreTreatedData->Get_PadNumber(i))+"_TIME_CAL"; + family = "Actar/CAL"; + + //GetHisto(family,name) -> Fill(PreTreatedData->Get_Time(i)); + } +} + + + +//////////////////////////////////////////////////////////////////////////////// +void TActarSpectra::FillPhysicsSpectra(TActarPhysics* Physics) { + static string name; + static string family; + family= "Actar/PHY"; + + // Energy vs time + unsigned int sizeE = Physics->PadCharge.size(); + for(unsigned int i = 0 ; i < sizeE ; i++){ + name = "Actar_ENERGY_TIME"; + //GetHisto(family,name) -> Fill(Physics->Charge[i],Physics->Time[i]); + } +} diff --git a/NPLib/Detectors/Actar/TActarSpectra.h b/NPLib/Detectors/Actar/TActarSpectra.h new file mode 100644 index 0000000000000000000000000000000000000000..2534a52b9cddddff2f8a65390a48bd26f4c2ed95 --- /dev/null +++ b/NPLib/Detectors/Actar/TActarSpectra.h @@ -0,0 +1,62 @@ +#ifndef TActarSPECTRA_H +#define TActarSPECTRA_H +/***************************************************************************** + * Copyright (C) 2009-2017 this file is part of the NPTool Project * + * * + * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * + * For the list of contributors see $NPTOOL/Licence/Contributors * + *****************************************************************************/ + +/***************************************************************************** + * Original Author: Pierre Morfouace contact address: morfouac@nscl.msu.edu * + * * + * Creation Date : September 2017 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class hold Actar Spectra * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +// NPLib headers +#include "NPVSpectra.h" +#include "TActarData.h" +#include "TActarPhysics.h" + +// Forward Declaration +class TActarPhysics; + + +class TActarSpectra : public VSpectra { + ////////////////////////////////////////////////////////////// + // constructor and destructor + public: + TActarSpectra(); + TActarSpectra(unsigned int NumberOfDetectors); + ~TActarSpectra(); + + ////////////////////////////////////////////////////////////// + // Initialization methods + private: + void InitRawSpectra(); + void InitPreTreatedSpectra(); + void InitPhysicsSpectra(); + + ////////////////////////////////////////////////////////////// + // Filling methods + public: + void FillRawSpectra(TActarData*); + void FillPreTreatedSpectra(TActarData*); + void FillPhysicsSpectra(TActarPhysics*); + + ////////////////////////////////////////////////////////////// + // Detector parameters + private: + unsigned int fNumberOfDetectors; +}; + +#endif diff --git a/NPLib/Detectors/Hira/THiraPhysics.cxx b/NPLib/Detectors/Hira/THiraPhysics.cxx index ccefeec09b206e6e32b2241520db22bff81fa57f..1ebaa897401b1f59c2fc7bc2e679eec420f4008e 100644 --- a/NPLib/Detectors/Hira/THiraPhysics.cxx +++ b/NPLib/Detectors/Hira/THiraPhysics.cxx @@ -758,9 +758,8 @@ int THiraPhysics::CheckEvent(){ return 1 ; // Regular Event // INterstrip management is not coded, so waste of time to make this test - /* else if( m_PreTreatedData->GetMMStripXEMult() == m_PreTreatedData->GetMMStripYEMult()+1 - || m_PreTreatedData->GetMMStripXEMult() == m_PreTreatedData->GetMMStripYEMult()-1 ) - return 2 ; // Pseudo Event, potentially interstrip*/ + else if( fabs(m_PreTreatedData->GetHiraStripXEMult() - m_PreTreatedData->GetHiraStripYEMult() )< 3) + return 2 ; // Pseudo Event, potentially interstrip events else return -1 ; // Rejected Event diff --git a/NPLib/Detectors/Tiara/TTiaraHyballPhysics.cxx b/NPLib/Detectors/Tiara/TTiaraHyballPhysics.cxx index 624a7e4abb11a9bf2d9c7522c46f8181c82336f1..37e09190bf80c6804376f05330e79e2bf13e1cec 100644 --- a/NPLib/Detectors/Tiara/TTiaraHyballPhysics.cxx +++ b/NPLib/Detectors/Tiara/TTiaraHyballPhysics.cxx @@ -70,6 +70,10 @@ ClassImp(TTiaraHyballPhysics) m_StripSector_E_RAW_Threshold = 0 ; m_StripSector_E_Threshold = 0.300 ; // MeV +//by Shuya 171019 + m_Take_E_Ring_Sector_Average=false; + m_Sigma_Ring_E = 0.018; // MeV + m_Sigma_Sector_E = 0.022; // MeV //by Shuya 170523. //m_Take_E_Ring=false; m_Take_E_Ring=true; @@ -124,10 +128,38 @@ void TTiaraHyballPhysics::BuildPhysicalEvent(){ StripSector_E.push_back(Sector_E) ; StripSector_T.push_back(Sector_T) ; +//by Shuya 170703 +/* if(m_Take_E_Ring) Strip_E.push_back(Ring_E) ; else Strip_E.push_back(Sector_E) ; +*/ + //by Shuya 171019 + if(m_Take_E_Ring_Sector_Average) + { + if(Ring_E && Sector_E) + { + double Weighted_E, Sigma_RingE2, Sigma_SectorE2; + double Sigma_RingE = m_Sigma_Ring_E; + double Sigma_SectorE = m_Sigma_Sector_E; + Sigma_RingE2 = m_Sigma_Ring_E * m_Sigma_Ring_E; + Sigma_SectorE2 = m_Sigma_Sector_E * m_Sigma_Sector_E; + + Weighted_E = (Ring_E / Sigma_RingE2 + Sector_E / Sigma_SectorE2) / (1.0 / Sigma_RingE2 + 1.0 / Sigma_SectorE2); + Strip_E.push_back(Weighted_E) ; + } + else if(Ring_E) Strip_E.push_back(Ring_E) ; + else Strip_E.push_back(Sector_E) ; + } + //if(m_Take_E_Ring) + else if(m_Take_E_Ring) + { + if(Ring_E) Strip_E.push_back(Ring_E) ; + else Strip_E.push_back(Sector_E) ; + } + else + Strip_E.push_back(Sector_E) ; if(m_Take_T_Sector) Strip_T.push_back(Sector_T) ; @@ -277,7 +309,9 @@ void TTiaraHyballPhysics::ReadAnalysisConfig(){ getline(AnalysisConfigFile, LineBuffer); // search for "header" - if (LineBuffer.compare(0, 11, "ConfigTiaraHyball") == 0) ReadingStatus = true; +//by Shuya 171019 + //if (LineBuffer.compare(0, 11, "ConfigTiaraHyball") == 0) ReadingStatus = true; + if (LineBuffer.compare(0, 17, "ConfigTiaraHyball") == 0) ReadingStatus = true; // loop on tokens and data while (ReadingStatus ) { @@ -338,6 +372,26 @@ void TTiaraHyballPhysics::ReadAnalysisConfig(){ } + //by Shuya 171019 + else if (whatToDo=="TAKE_E_RING_SECTOR_Average") { + m_Take_E_Ring_Sector_Average = true; + cout << whatToDo << endl; + } + + //by Shuya 171019 + else if (whatToDo=="SIGMA_RING_E") { + AnalysisConfigFile >> DataBuffer; + m_Sigma_Ring_E = atof(DataBuffer.c_str() ); + cout << "STANDARD DEVIATION OF RING ENERGY " << m_Sigma_Ring_E <<endl; + } + + //by Shuya 171019 + else if (whatToDo=="SIGMA_SECTOR_E") { + AnalysisConfigFile >> DataBuffer; + m_Sigma_Sector_E = atof(DataBuffer.c_str() ); + cout << "STANDARD DEVIATION OF SECTOR ENERGY " << m_Sigma_Sector_E <<endl; + } + else if (whatToDo=="TAKE_E_RING") { m_Take_E_Ring = true; cout << whatToDo << endl; diff --git a/NPLib/Detectors/Tiara/TTiaraHyballPhysics.h b/NPLib/Detectors/Tiara/TTiaraHyballPhysics.h index b1da351c394c31eafbc402a80a281d8f6d266c6c..abcf2c1a6b88c306111910db48eb9a05dd466915 100644 --- a/NPLib/Detectors/Tiara/TTiaraHyballPhysics.h +++ b/NPLib/Detectors/Tiara/TTiaraHyballPhysics.h @@ -167,6 +167,10 @@ class TTiaraHyballPhysics : public TObject, public NPL::VDetector{ // By default take EX and TY. bool m_Take_E_Ring;//! bool m_Take_T_Sector;//! +//by Shuya 171019 + bool m_Take_E_Ring_Sector_Average;//! + double m_Sigma_Ring_E; + double m_Sigma_Sector_E; // Event over this value after pre-treatment are not treated / avoid long treatment time on spurious event unsigned int m_MaximumStripMultiplicityAllowed ;//! diff --git a/NPLib/Utility/npanalysis.cxx b/NPLib/Utility/npanalysis.cxx index 4a66050a999a2b2b6b3cdff2c0bbad43a8f2b408..061bfbb413cd0a1f3d65b64befe6c8deee6f8ae0 100644 --- a/NPLib/Utility/npanalysis.cxx +++ b/NPLib/Utility/npanalysis.cxx @@ -168,15 +168,19 @@ int main(int argc , char** argv){ if(!IsPhysics){ for (unsigned int i = 0 ; i < nentries; i++) { // Get the raw Data + //cout << "!" << endl; Chain -> GetEntry(i); - Chain -> GetEntry(i); + //cout << "!!" << endl; // Build the current event myDetector->BuildPhysicalEvent(); + //cout << "!!!" << endl; // User Analysis UserAnalysis->TreatEvent(); + //cout << "!!!!" << endl; // Fill the tree tree->Fill(); + //cout << "!!!!!" << endl; current_tree = Chain->GetTreeNumber()+1; ProgressDisplay(begin,end,treated,inter,nentries,mean_rate,displayed,current_tree,total_tree); diff --git a/NPSimulation/Core/MaterialManager.cc b/NPSimulation/Core/MaterialManager.cc index 29a0653db4a713880c185f7e8e6a59c056c6e437..29bdcad13612925cf79e1a5af4412854b96b4ed7 100644 --- a/NPSimulation/Core/MaterialManager.cc +++ b/NPSimulation/Core/MaterialManager.cc @@ -742,7 +742,10 @@ G4Material* MaterialManager::GetGasFromLibrary(string Name, double Pressure, dou string newName= oss.str(); map<string,G4Material*>::iterator it; it = m_Material.find(Name); - double density = 0 ; + double density = 0 ; + + G4double Vm=0.08206*Temperature*atmosphere/(Pressure*kelvin); + // The element is not found if(it==m_Material.end()){ if(Name == "CF4"){ // 52 torr @@ -756,7 +759,61 @@ G4Material* MaterialManager::GetGasFromLibrary(string Name, double Pressure, dou m_Material[Name]=material; return material; } + + if(Name == "He"){ + density = (4.0026/Vm)*mg/cm3; + G4Material* material = new G4Material("NPS_"+newName,density,1,kStateGas,Temperature,Pressure); + material->AddElement(GetElementFromLibrary("He"), 1); + m_Material[Name]=material; + return material; + } + + if(Name == "iC4H10" || Name == "Isobutane" || Name == "isobutane"){ + density = ((4*12.0107+10*1.00794)/Vm)*mg/cm3; + G4Material* material = new G4Material("NPS_"+newName,density,2,kStateGas,Temperature,Pressure); + material->AddElement(GetElementFromLibrary("C"), 4); + material->AddElement(GetElementFromLibrary("H"), 10); + m_Material[Name]=material; + return material; + } + + if(Name == "CH4"){ + density = ((12.0107+4*1.00794)/Vm)*mg/cm3; + G4Material* material = new G4Material("NPS_"+newName,density,2,kStateGas,Temperature,Pressure); + material->AddElement(GetElementFromLibrary("C"), 1); + material->AddElement(GetElementFromLibrary("H"), 4); + m_Material[Name]=material; + return material; + } + + if(Name == "CO2"){ + density = ((12.0107+2*16)/Vm)*mg/cm3; + G4Material* material = new G4Material("NPS_"+newName,density,2,kStateGas,Temperature,Pressure); + material->AddElement(GetElementFromLibrary("C"), 1); + material->AddElement(GetElementFromLibrary("O"), 2); + m_Material[Name]=material; + return material; + } + + if(Name == "H2"){ + density = (2*1.00794/Vm)*mg/cm3; + G4Material* material = new G4Material("NPS_"+newName,density,2,kStateGas,Temperature,Pressure); + material->AddElement(GetElementFromLibrary("H"), 1); + material->AddElement(GetElementFromLibrary("H"), 1); + m_Material[Name]=material; + return material; + } + + if(Name == "D2"){ + density = (2*2.0140/Vm)*mg/cm3; + G4Material* material = new G4Material("NPS_"+newName,density,2,kStateGas,Temperature,Pressure); + material->AddElement(GetElementFromLibrary("D"), 1); + material->AddElement(GetElementFromLibrary("D"), 1); + m_Material[Name]=material; + return material; + } + else{ exit(1); } diff --git a/NPSimulation/Core/ParticleStack.cc b/NPSimulation/Core/ParticleStack.cc index 1d7ac6c1711117c6d893bb7acba571708d44c4b7..236eab375ff5d25d7c3c03017806295546578e9f 100644 --- a/NPSimulation/Core/ParticleStack.cc +++ b/NPSimulation/Core/ParticleStack.cc @@ -106,8 +106,7 @@ void ParticleStack::AddBeamParticleToStack(Particle& particle){ m_ParticleStack.push_back(particle); // Incident beam parameter m_InitialConditions-> SetIncidentParticleName (particle.GetParticleDefinition()->GetParticleName()); - //m_InitialConditions-> SetIncidentInitialKineticEnergy (particle. GetParticleThetaCM()); - m_InitialConditions-> SetIncidentInitialKineticEnergy (particle. GetParticleKineticEnergy()); + m_InitialConditions-> SetIncidentFinalKineticEnergy (particle. GetParticleKineticEnergy()); G4ThreeVector U(1,0,0); G4ThreeVector V(0,1,0); @@ -118,7 +117,7 @@ void ParticleStack::AddBeamParticleToStack(Particle& particle){ m_InitialConditions-> SetIncidentEmittancePhi (particle.GetParticleMomentumDirection().phi()/deg); // Beam status at the initial interaction point - m_InitialConditions-> SetIncidentFinalKineticEnergy (particle. GetParticleKineticEnergy()); + m_InitialConditions-> SetIncidentInitialKineticEnergy (particle. GetParticleThetaCM()); m_InitialConditions-> SetIncidentPositionX (particle. GetParticlePosition().x()); m_InitialConditions-> SetIncidentPositionY (particle. GetParticlePosition().y()); m_InitialConditions-> SetIncidentPositionZ (particle. GetParticlePosition().z()); diff --git a/NPSimulation/Detectors/Actar/Actar.cc b/NPSimulation/Detectors/Actar/Actar.cc new file mode 100644 index 0000000000000000000000000000000000000000..741f247d34fc5ee8a031fd414568268824864a6c --- /dev/null +++ b/NPSimulation/Detectors/Actar/Actar.cc @@ -0,0 +1,691 @@ +/***************************************************************************** + * Copyright (C) 2009-2017 this file is part of the NPTool Project * + * * + * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * + * For the list of contributors see $NPTOOL/Licence/Contributors * + *****************************************************************************/ + +/***************************************************************************** + * Original Author: Pierre Morfouace contact address: morfouac@nscl.msu.edu * + * * + * Creation Date : September 2017 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class describe Actar simulation * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + *****************************************************************************/ + +// C++ headers +#include <sstream> +#include <cmath> +#include <limits> +//G4 Geometry object +#include "G4Tubs.hh" +#include "G4Box.hh" + +//G4 sensitive +#include "G4SDManager.hh" +#include "G4MultiFunctionalDetector.hh" + +//G4 various object +#include "G4Material.hh" +#include "G4Transform3D.hh" +#include "G4PVPlacement.hh" +#include "G4VisAttributes.hh" +#include "G4Colour.hh" + +// G4 Field +#include "G4FieldManager.hh" +#include "G4ElectricField.hh" +#include "G4UniformElectricField.hh" +#include "G4TransportationManager.hh" +#include "G4EqMagElectricField.hh" +#include "G4MagIntegratorStepper.hh" +#include "G4ClassicalRK4.hh" +#include "G4MagIntegratorDriver.hh" +#include "G4ChordFinder.hh" +#include "G4MaterialPropertiesTable.hh" + +// NPTool header +#include "Actar.hh" +#include "SiliconScorers.hh" +#include "DriftElectronScorers.hh" +#include "RootOutput.h" +#include "MaterialManager.hh" +#include "NPSDetectorFactory.hh" +#include "NPOptionManager.h" +#include "FastDriftElectron.hh" +#include "NPSHitsMap.hh" +#include "CalorimeterScorers.hh" + + +// CLHEP header +#include "CLHEP/Random/RandGauss.h" + +// ROOT +#include "TH1D.h" +#include "TF1.h" + +using namespace std; +using namespace CLHEP; + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +namespace Actar_NS{ + // Energy and time Resolution + const double ChargeThreshold = 0; + const double ResoTime = 0.1*ns ; + //const double ResoEnergy = 1.0*MeV ; + const double ChamberThickness = 55*cm ; + const double ChamberWidth = 50*cm ; + const double ChamberHeight = 40*cm ; + //const int NumberOfPads = 16384; + const int PadX = 128; + const int PadZ = 128; + + const double Nose_Rmin = 2.5*cm; + const double Nose_Rmax = 3.5*cm; + const double Nose_Length = 12*cm; + + const double Mylar_Rmax = 3.5*cm; + const double Mylar_Thickness = 7*micrometer; + + const double XGazVolume = 256.*mm; + const double YGazVolume = 256.*mm; + const double ZGazVolume = 256.*mm; + + const double SiliconHeight = 53.*mm; + const double SiliconWidth = 53.*mm; + const double SiliconThickness = 0.7*mm; + const double DistInterSi = 1.*mm; + const double Si_PosZ=15.*cm; + const double ResoSilicon = 0.80/2.35; + const double EnergyThreshold = 0.1; + + const double CsIThickness = 1.*cm; + const double CsIHeight = 2.5*cm; + const double CsIWidth = 2.5*cm; + const double DistInterCsI = 1.*mm; + const double CsI_PosZ = 16.*cm; + const double ResoCsI = 0.200/2.35; +} +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// Actar Specific Method +Actar::Actar(){ + m_Event = new TActarData() ; + m_ActarScorer = 0; + m_SquareDetector = 0; + + // RGB Color + Transparency + m_VisChamber = new G4VisAttributes(G4Colour(0.7, 0.7, 0.7, 0.3)); + m_VisWindows = new G4VisAttributes(G4Colour(1, 0, 0, 0.25)); + m_VisGas = new G4VisAttributes(G4Colour(0, 0.5, 0.5, 0.3)); + m_VisPads = new G4VisAttributes(G4Colour(255, 223, 50, 0.8)); + m_SiliconVisAtt = new G4VisAttributes(G4Colour(0.529412, 0.807843, 0.980392, 0.95)) ; + m_CsIVisAtt = new G4VisAttributes(G4Colour(0.429412, 0.607843, 0.780392, 0.95)); + m_VisPads->SetForceWireframe(true); + + m_build_Silicon=1; + m_build_CsI=1; +} + +Actar::~Actar(){ +} +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +void Actar::AddDetector(G4ThreeVector POS, string Shape){ + // Convert the POS value to R theta Phi as Spherical coordinate is easier in G4 + m_R.push_back(POS.mag()); + m_Theta.push_back(POS.theta()); + m_Phi.push_back(POS.phi()); + m_Shape.push_back(Shape); +} + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +void Actar::AddDetector(double R, double Theta, double Phi, string Shape){ + m_R.push_back(R); + m_Theta.push_back(Theta); + m_Phi.push_back(Phi); + m_Shape.push_back(Shape); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +G4LogicalVolume* Actar::BuildDetector(){ + + G4Material* Cu= MaterialManager::getInstance()->GetMaterialFromLibrary("Cu"); + G4Material* Si= MaterialManager::getInstance()->GetMaterialFromLibrary("Si"); + G4Material* Al= MaterialManager::getInstance()->GetMaterialFromLibrary("Al"); + G4Material* Mylar= MaterialManager::getInstance()->GetMaterialFromLibrary("Mylar"); + G4Material* MaterialCsI = MaterialManager::getInstance()->GetMaterialFromLibrary("CsI"); + + if(!m_SquareDetector){ + // Main volume + G4Box* sChamber = new G4Box("Actar_Box",Actar_NS::ChamberWidth*0.5, + Actar_NS::ChamberHeight*0.5,Actar_NS::ChamberThickness*0.5); + + //Nose volume + G4Tubs* sNose = new G4Tubs("Actar_Nose",Actar_NS::Nose_Rmin, Actar_NS::Nose_Rmax,Actar_NS::Nose_Length*0.5, + 0*deg, 360*deg); + + //Mylar volume + G4Tubs* sWindows = new G4Tubs("Actar_Windows",0, Actar_NS::Mylar_Rmax,Actar_NS::Mylar_Thickness*0.5, + 0*deg, 360*deg); + + // Cage volume + G4Box* sCage = new G4Box("Actar_Gas",Actar_NS::XGazVolume*0.5, + Actar_NS::YGazVolume*0.5,Actar_NS::ZGazVolume*0.5); + + // Pad + G4Box* sPad = new G4Box("Actar_Pad",2*mm*0.5, + 1*um*0.5,2*mm*0.5); + + // Cathode + G4Box* sCathode = new G4Box("Actar_Cathode",26.5*cm*0.5, + 1*um*0.5,25.6*cm*0.5); + + + + unsigned const int NumberOfGasMix = m_GasMaterial.size(); + + double density=0; + double density_sum=0; + vector<G4Material*> GasComponent; + vector<double> FractionMass; + + for(unsigned int i=0; i<NumberOfGasMix; i++){ + GasComponent.push_back(MaterialManager::getInstance()->GetGasFromLibrary(m_GasMaterial[i],m_Pressure,m_Temperature) ); + } + for(unsigned int i=0; i<NumberOfGasMix; i++){ + density += ((double)m_GasFraction[i]/100)*GasComponent[i]->GetDensity(); + density_sum += GasComponent[i]->GetDensity(); + } + //cout << "density = " << density*cm3/g << endl; + + for(unsigned int i=0; i<NumberOfGasMix; i++){ + FractionMass.push_back(GasComponent[i]->GetDensity()/density_sum); + } + + G4Material* GasMaterial = new G4Material("GasMix", density, NumberOfGasMix, kStateGas, m_Temperature, m_Pressure); + G4Material* DriftGasMaterial = new G4Material("GasMix", density, NumberOfGasMix, kStateGas, m_Temperature, m_Pressure); + + for(unsigned int i=0; i<NumberOfGasMix; i++){ + GasMaterial->AddMaterial(GasComponent[i], FractionMass[i]); + DriftGasMaterial->AddMaterial(GasComponent[i], FractionMass[i]); + } + + G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable(); + MPT->AddConstProperty("DE_PAIRENERGY",20*eV); + MPT->AddConstProperty("DE_YIELD",1e-2); + //MPT->AddConstProperty("DE_YIELD",1e-1); + //MPT->AddConstProperty("DE_AMPLIFICATION",5); + MPT->AddConstProperty("DE_ABSLENGTH",1*pc); + MPT->AddConstProperty("DE_DRIFTSPEED",5.*cm/microsecond); + //MPT->AddConstProperty("DE_TRANSVERSALSPREAD",5e-5*mm2/ns); + //MPT->AddConstProperty("DE_LONGITUDINALSPREAD",5e-5*mm2/ns); + MPT->AddConstProperty("DE_TRANSVERSALSPREAD",7e-6*mm2/ns); + MPT->AddConstProperty("DE_LONGITUDINALSPREAD",7e-6*mm2/ns); + + DriftGasMaterial->SetMaterialPropertiesTable(MPT); + + G4MaterialPropertiesTable* MPT2 = new G4MaterialPropertiesTable(); + MPT2->AddConstProperty("DE_YIELD",1); + MPT2->AddConstProperty("DE_AMPLIFICATION",2); + MPT2->AddConstProperty("DE_ABSLENGTH",1*pc); + + Al->SetMaterialPropertiesTable(MPT2); + + m_SquareDetector = new G4LogicalVolume(sChamber,GasMaterial,"logic_Actar_Box",0,0,0); + G4LogicalVolume* logicGas = new G4LogicalVolume(sCage,DriftGasMaterial,"logic_Gas",0,0,0); + G4LogicalVolume* logicNose = new G4LogicalVolume(sNose,Al,"logic_Nose",0,0,0); + G4LogicalVolume* logicPad = new G4LogicalVolume(sPad,Cu,"logic_Pad",0,0,0); + G4LogicalVolume* logicCathode = new G4LogicalVolume(sCathode,Cu,"logic_Cathode",0,0,0); + G4LogicalVolume* logicWindows = new G4LogicalVolume(sWindows,Mylar,"logic_Windows",0,0,0); + + G4RotationMatrix* Rot = new G4RotationMatrix(); + new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,0,-Actar_NS::ChamberThickness*0.5+Actar_NS::Nose_Length*0.5)), + logicNose, + "ActarNose",m_SquareDetector,false,0); + + new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,0,-Actar_NS::ChamberThickness*0.5+Actar_NS::Mylar_Thickness*0.5+Actar_NS::Nose_Length)), + logicWindows, + "ActarEntranceWindows",m_SquareDetector,false,0); + + new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,0,0)), + logicGas, + "ActarGas",m_SquareDetector,false,0); + + int pad=0; + m_PadToXRow.clear(); + m_PadToZColumn.clear(); + for(int i=0; i<Actar_NS::PadX; i++){ + for(int j=0; j<Actar_NS::PadZ; j++){ + m_PadToXRow[pad] = i; + m_PadToZColumn[pad] = j; + double X=(i-64)*2*mm; + double Z=(j-64)*2*mm; + new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(X,Actar_NS::YGazVolume*0.5,Z)), + logicPad, + "ActarPad",logicGas,false,pad+1); + + pad++; + } + } + + /*new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(3*cm-0.5*1*um,0,0)), + logicCathode, + "ActarCathode",logicGas,false,0); + + + + new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,0,-6*cm+6*micrometer)), + logicWindows, + "ActarEntranceWindows",m_SquareDetector,false,0);*/ + + G4ElectricField* field = new G4UniformElectricField(G4ThreeVector(0.0,-100*volt/cm,0.0)); + // Create an equation of motion for this field + G4EqMagElectricField* Equation = new G4EqMagElectricField(field); + G4MagIntegratorStepper* Stepper = new G4ClassicalRK4( Equation, 8 ); + + // Get the global field manager + G4FieldManager* FieldManager= new G4FieldManager(); + // Set this field to the global field manager + FieldManager->SetDetectorField(field); + logicGas->SetFieldManager(FieldManager,true); + + G4MagInt_Driver* IntgrDriver = new G4MagInt_Driver(0.1*mm, + Stepper, + Stepper->GetNumberOfVariables() ); + + G4ChordFinder* ChordFinder = new G4ChordFinder(IntgrDriver); + FieldManager->SetChordFinder( ChordFinder ); + + + logicPad->SetSensitiveDetector(m_ActarScorer); + logicNose->SetVisAttributes(m_VisChamber); + m_SquareDetector->SetVisAttributes(m_VisChamber); + logicGas->SetVisAttributes(m_VisGas); + logicWindows->SetVisAttributes(m_VisWindows); + logicPad->SetVisAttributes(m_VisPads); + //m_SquareDetector->SetSensitiveDetector(m_ActarScorer); + } + + + /////////////////////////////////////////////////// + ///////////////////// Thin Si ///////////////////// + /////////////////////////////////////////////////// + if(m_build_Silicon){ + G4Box* solidSi = new G4Box("Si", 0.5*Actar_NS::SiliconWidth, 0.5*Actar_NS::SiliconHeight, 0.5*Actar_NS::SiliconThickness); ; + m_LogicSilicon = new G4LogicalVolume(solidSi, Si, "logicSi", 0, 0, 0); + + int SiliconNumber=0; + for(int k=0;k<4; k++){ + for(int p=0; p<5; p++){ + double PosX; + double PosY; + if(k==0) PosY= -1.5*Actar_NS::SiliconHeight-2*Actar_NS::DistInterSi; + if(k==1) PosY= -0.5*Actar_NS::SiliconHeight-1*Actar_NS::DistInterSi; + if(k==2) PosY= 0.5*Actar_NS::SiliconHeight+1*Actar_NS::DistInterSi; + if(k==3) PosY= 1.5*Actar_NS::SiliconHeight+2*Actar_NS::DistInterSi; + if(p==0) PosX= -2*Actar_NS::SiliconWidth-2*Actar_NS::DistInterSi; + if(p==1) PosX= -1*Actar_NS::SiliconWidth-1*Actar_NS::DistInterSi; + if(p==2) PosX= 0; + if(p==3) PosX= 1*Actar_NS::SiliconWidth+2*Actar_NS::DistInterSi; + if(p==4) PosX= 2*Actar_NS::SiliconWidth+1*Actar_NS::DistInterSi; + + G4ThreeVector positionSi = G4ThreeVector(PosX, PosY, Actar_NS::Si_PosZ); + new G4PVPlacement(new G4RotationMatrix(0,0,0), + positionSi, + m_LogicSilicon,"Si", + m_SquareDetector,false,SiliconNumber); + SiliconNumber++; + } + } + + // Set Si sensible + m_LogicSilicon->SetSensitiveDetector(m_SiliconScorer); + + // Visualisation of ThinSi + m_LogicSilicon->SetVisAttributes(m_SiliconVisAtt) ; + } + + /////////////////////////////////////////////// + ///////////////////// CsI ///////////////////// + /////////////////////////////////////////////// + if(m_build_CsI){ + G4Box* solidCsI = new G4Box("Si", 0.5*Actar_NS::CsIWidth, 0.5*Actar_NS::CsIHeight, 0.5*Actar_NS::CsIThickness); ; + m_LogicCsICrystal = new G4LogicalVolume(solidCsI, MaterialCsI, "logicCsI", 0, 0, 0); + + int CsINumber=0; + for(int k=0;k<8; k++){ + for(int p=0; p<10; p++){ + double PosX; + double PosY; + if(k<4) PosY= -0.5*Actar_NS::CsIHeight-0.5*Actar_NS::DistInterCsI+(k-3)*(Actar_NS::CsIHeight+Actar_NS::DistInterCsI); + if(k>3) PosY= 0.5*Actar_NS::CsIHeight+0.5*Actar_NS::DistInterCsI+(k-4)*(Actar_NS::CsIHeight+Actar_NS::DistInterCsI); + + if(p<5) PosX= 0.5*Actar_NS::CsIWidth+0.5*Actar_NS::DistInterCsI+(4-p)*(Actar_NS::CsIWidth+Actar_NS::DistInterCsI); + if(p>4) PosX= -0.5*Actar_NS::CsIWidth-0.5*Actar_NS::DistInterCsI+(5-p)*(Actar_NS::CsIWidth+Actar_NS::DistInterCsI); + + G4ThreeVector positionCsI = G4ThreeVector(PosX, PosY, Actar_NS::CsI_PosZ); + new G4PVPlacement(new G4RotationMatrix(0,0,0), + positionCsI, + m_LogicCsICrystal,"CsI", + m_SquareDetector,false,CsINumber); + CsINumber++; + } + } + + // Set Si sensible + m_LogicCsICrystal->SetSensitiveDetector(m_CsIScorer); + + // Visualisation of ThinSi + m_LogicCsICrystal->SetVisAttributes(m_CsIVisAtt) ; + } + + return m_SquareDetector; +} + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// Virtual Method of NPS::VDetector class + +// Read stream at Configfile to pick-up parameters of detector (Position,...) +// Called in DetecorConstruction::ReadDetextorConfiguration Method +void Actar::ReadConfiguration(NPL::InputParser parser){ + vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Actar"); + if(NPOptionManager::getInstance()->GetVerboseLevel()) + cout << "//// " << blocks.size() << " detectors found " << endl; + + vector<string> cart = {"POS","Shape","GasMaterial","GasFraction","Temperature","Pressure","Si","CsI"}; + vector<string> sphe = {"R","Theta","Phi","Shape","GasMaterial","GasFraction","Temperature","Pressure","Si","CsI"}; + + for(unsigned int i = 0 ; i < blocks.size() ; i++){ + if(blocks[i]->HasTokenList(cart)){ + if(NPOptionManager::getInstance()->GetVerboseLevel()) + cout << endl << "//// Actar " << i+1 << endl; + + G4ThreeVector Pos = NPS::ConvertVector(blocks[i]->GetTVector3("POS","mm")); + string Shape = blocks[i]->GetString("Shape"); + vector<string> GasName = blocks[i]->GetVectorString("GasMaterial"); + vector<int> GasFraction = blocks[i]->GetVectorInt("GasFraction"); + for(unsigned int j=0; j<GasName.size(); j++){ + m_GasMaterial.push_back(GasName[j]); + m_GasFraction.push_back(GasFraction[j]); + } + m_Temperature = blocks[i]->GetDouble("Temperature","kelvin"); + m_Pressure = blocks[i]->GetDouble("Pressure","bar"); + m_build_Silicon = blocks[i]->GetInt("Si"); + m_build_CsI = blocks[i]->GetInt("CsI"); + + AddDetector(Pos,Shape); + } + else if(blocks[i]->HasTokenList(sphe)){ + if(NPOptionManager::getInstance()->GetVerboseLevel()) + cout << endl << "//// Actar " << i+1 << endl; + double R = blocks[i]->GetDouble("R","mm"); + double Theta = blocks[i]->GetDouble("Theta","deg"); + double Phi = blocks[i]->GetDouble("Phi","deg"); + string Shape = blocks[i]->GetString("Shape"); + vector<string> GasName = blocks[i]->GetVectorString("GasMaterial"); + vector<int> GasFraction = blocks[i]->GetVectorInt("GasFraction"); + for(unsigned int j=0; j<GasName.size(); j++){ + m_GasMaterial.push_back(GasName[j]); + m_GasFraction.push_back(GasFraction[j]); + } + m_Temperature = blocks[i]->GetDouble("Temperature","kelvin"); + m_Pressure = blocks[i]->GetDouble("Pressure","bar"); + m_build_Silicon = blocks[i]->GetInt("Si"); + m_build_CsI = blocks[i]->GetInt("CsI"); + + AddDetector(R,Theta,Phi,Shape); + } + else{ + cout << "ERROR: check your input file formatting " << endl; + exit(1); + } + } +} + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +// Construct detector and inialise sensitive part. +// Called After DetecorConstruction::AddDetector Method +void Actar::ConstructDetector(G4LogicalVolume* world){ + for (unsigned short i = 0 ; i < m_R.size() ; i++) { + + G4double wX = m_R[i] * sin(m_Theta[i] ) * cos(m_Phi[i] ) ; + G4double wY = m_R[i] * sin(m_Theta[i] ) * sin(m_Phi[i] ) ; + G4double wZ = m_R[i] * cos(m_Theta[i] ) ; + G4ThreeVector Det_pos = G4ThreeVector(wX, wY, wZ) ; + // So the face of the detector is at R instead of the middle + Det_pos+=Det_pos.unit()*Actar_NS::ChamberThickness*0.5; + // Building Detector reference frame + G4double ii = cos(m_Theta[i]) * cos(m_Phi[i]); + G4double jj = cos(m_Theta[i]) * sin(m_Phi[i]); + G4double kk = -sin(m_Theta[i]); + G4ThreeVector Y(ii,jj,kk); + G4ThreeVector w = Det_pos.unit(); + G4ThreeVector u = w.cross(Y); + G4ThreeVector v = w.cross(u); + v = v.unit(); + u = u.unit(); + + G4RotationMatrix* Rot = new G4RotationMatrix(u,v,w); + + new G4PVPlacement(G4Transform3D(*Rot,Det_pos), + BuildDetector(), + "Actar",world,false,i+1); + } +} +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// Add Detector branch to the EventTree. +// Called After DetecorConstruction::AddDetector Method +void Actar::InitializeRootOutput(){ + RootOutput *pAnalysis = RootOutput::getInstance(); + TTree *pTree = pAnalysis->GetTree(); + if(!pTree->FindBranch("Actar")){ + pTree->Branch("Actar", "TActarData", &m_Event) ; + } + pTree->SetBranchAddress("Actar", &m_Event) ; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// Read sensitive part and fill the Root tree. +// Called at in the EventAction::EndOfEventAvtion +void Actar::ReadSensitive(const G4Event* event){ + m_Event->Clear(); + + ////////////// + // Pad scorer + NPS::HitsMap<G4double*>* PadHitMap; + std::map<G4int, G4double**>::iterator Pad_itr; + + G4int PadCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("ActarScorer/Actar_dig"); + PadHitMap = (NPS::HitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(PadCollectionID)); + + // Loop on the Pad map + TH1D* h = new TH1D("h","h",25000,0,25000); + for (Pad_itr = PadHitMap->GetMap()->begin() ; Pad_itr != PadHitMap->GetMap()->end() ; Pad_itr++){ + G4double* Info = *(Pad_itr->second); + // Interraction Coordinates + ms_InterCoord->SetDetectedPositionX(Info[2]) ; + ms_InterCoord->SetDetectedPositionY(Info[3]) ; + ms_InterCoord->SetDetectedPositionZ(Info[4]) ; + ms_InterCoord->SetDetectedAngleTheta(Info[5]/deg) ; + ms_InterCoord->SetDetectedAnglePhi(Info[6]/deg) ; + + double Count = Info[0]; + double Time = RandGauss::shoot(Info[1],Actar_NS::ResoTime); + //int iTime = ((int) Time*20/512)+1; + int PadNbr = Info[7]; + if(Count>Actar_NS::ChargeThreshold){ + m_Event->SetTime(Time); + m_Event->SetPadNumber(PadNbr); + m_Event->SetCharge(Count); + m_Event->SetRowNumber(m_PadToXRow[PadNbr]); + m_Event->SetColumnNumber(m_PadToZColumn[PadNbr]); + } + + if(Count){ + h->Fill(Info[1],Info[0]); + } + } + + vector<double> Q, T; + for(int i=0; i<h->GetNbinsX(); i++){ + double count = h->GetBinContent(i); + double time = h->GetBinCenter(i); + if(count){ + Q.push_back(count); + T.push_back(time+500); + + } + } + // clear map for next event + //SimulateDigitizer(Q,T,1.40*microsecond,0,8750,25,5); + delete h; + PadHitMap->clear(); + + // Silicon // + if(m_build_Silicon){ + NPS::HitsMap<G4double*>* SiHitMap; + std::map<G4int, G4double**>::iterator Si_itr; + + G4int SiCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("SiliconScorer/SiliconScorer"); + SiHitMap = (NPS::HitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(SiCollectionID)); + + // Loop on the ThinSi map + for (Si_itr = SiHitMap->GetMap()->begin() ; Si_itr != SiHitMap->GetMap()->end() ; Si_itr++){ + G4double* Info = *(Si_itr->second); + double E_Si = RandGauss::shoot(Info[0],Actar_NS::ResoSilicon); + + if(E_Si>Actar_NS::EnergyThreshold){ + m_Event->SetSiliconEnergy(E_Si); + m_Event->SetSiliconDetectorNumber(Info[7]); + m_Event->SetSiliconTime(Info[1]); + } + } + + // Clear Map for next event + SiHitMap->clear(); + } + + // CsI // + if(m_build_CsI){ + NPS::HitsMap<G4double*>* CsIHitMap; + std::map<G4int, G4double**>::iterator CsI_itr; + + G4int CsICollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("CsIScorer/CsI"); + CsIHitMap = (NPS::HitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(CsICollectionID)); + + // Loop on the CsI map + for (CsI_itr = CsIHitMap->GetMap()->begin() ; CsI_itr !=CsIHitMap->GetMap()->end() ; CsI_itr++){ + G4double* Info = *(CsI_itr->second); + double E_CsI = RandGauss::shoot(Info[0],Actar_NS::ResoCsI); + + if(E_CsI>Actar_NS::EnergyThreshold){ + m_Event->SetCsIEnergy(E_CsI); + m_Event->SetCsICrystalNumber(Info[2]); + } + } + // Clear Map for next event + CsIHitMap->clear(); + } +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +void Actar::SimulateDigitizer(vector<double> E, vector<double> T, double fallTime,double start,double stop, double step,double noise){ + + static string formula; + formula= ""; + static string Es,Ts,var,cond; + static string fall; + fall=std::to_string(fallTime); + + for(unsigned int i = 0 ; i < E.size() ; i++){ + if(E[i]!=0 && T[i]!=0){ + Es = std::to_string(E[i]); + Ts = std::to_string(T[i]); + cond = ")*(x>"+Ts+")+"; + var = "(x-"+Ts+")"; + formula += Es+"*-1*exp(-"+var+"/"+fall+cond; + } + } + formula+="0"; + //cout << formula << endl; + TF1* f = new TF1("f",formula.c_str(),start,stop); + unsigned int size = (stop-start)/step; + for(unsigned int i = 0 ; i < size ; i++){ + double time = start+i*step; + double energy = f->Eval(time)+noise*(1-2*G4UniformRand()); + m_Event->AddEnergyPoint(energy,time); + } + + delete f; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//////////////////////////////////////////////////////////////// +void Actar::InitializeScorers() { + // This check is necessary in case the geometry is reloaded + + bool already_exist = false; + vector<G4int> NestingLevel; + NestingLevel.push_back(0); + NestingLevel.push_back(2); + + m_ActarScorer = CheckScorer("ActarScorer",already_exist) ; + m_SiliconScorer = CheckScorer("SiliconScorer",already_exist); + m_CsIScorer = CheckScorer("CsIScorer",already_exist); + + if(already_exist) return; + + G4VPrimitiveScorer *SiScorer = new SILICONSCORERS::PS_Silicon_Rectangle("SiliconScorer",0,Actar_NS::SiliconHeight,Actar_NS::SiliconWidth,1,1); + m_SiliconScorer->RegisterPrimitive(SiScorer); + + G4VPrimitiveScorer* CsIScorer= new CALORIMETERSCORERS::PS_Calorimeter("CsI",NestingLevel); + m_CsIScorer->RegisterPrimitive(CsIScorer); + + vector<int> level; level.push_back(0); + G4VPrimitiveScorer* Actar_dig= new DRIFTELECTRONSCORERS::PS_DECathode("Actar_dig",0) ; + m_ActarScorer->RegisterPrimitive(Actar_dig); + + G4SDManager::GetSDMpointer()->AddNewDetector(m_ActarScorer); + G4SDManager::GetSDMpointer()->AddNewDetector(m_SiliconScorer); + G4SDManager::GetSDMpointer()->AddNewDetector(m_CsIScorer) ; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//////////////////////////////////////////////////////////////////////////////// +// Construct Method to be pass to the DetectorFactory // +//////////////////////////////////////////////////////////////////////////////// +NPS::VDetector* Actar::Construct(){ + return (NPS::VDetector*) new Actar(); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//////////////////////////////////////////////////////////////////////////////// +// Registering the construct method to the factory // +//////////////////////////////////////////////////////////////////////////////// +extern"C" { + class proxy_nps_Actar{ + public: + proxy_nps_Actar(){ + NPS::DetectorFactory::getInstance()->AddToken("Actar","Actar"); + NPS::DetectorFactory::getInstance()->AddDetector("Actar",Actar::Construct); + } + }; + + proxy_nps_Actar p_nps_Actar; +} diff --git a/NPSimulation/Detectors/Actar/Actar.hh b/NPSimulation/Detectors/Actar/Actar.hh new file mode 100644 index 0000000000000000000000000000000000000000..43aecfd54d1a928425bb60f9c4be53edc3825a04 --- /dev/null +++ b/NPSimulation/Detectors/Actar/Actar.hh @@ -0,0 +1,136 @@ +#ifndef Actar_h +#define Actar_h 1 +/***************************************************************************** + * Copyright (C) 2009-2017 this file is part of the NPTool Project * + * * + * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * + * For the list of contributors see $NPTOOL/Licence/Contributors * + *****************************************************************************/ + +/***************************************************************************** + * Original Author: Pierre Morfouace contact address: morfouac@nscl.msu.edu * + * * + * Creation Date : September 2017 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class describe Actar simulation * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + *****************************************************************************/ + +// C++ header +#include <string> +#include <vector> +#include <map> +using namespace std; + +// G4 headers +#include "G4ThreeVector.hh" +#include "G4RotationMatrix.hh" +#include "G4LogicalVolume.hh" +#include "G4MultiFunctionalDetector.hh" + +// NPTool header +#include "NPSVDetector.hh" +#include "TActarData.h" +#include "NPInputParser.h" + +class Actar : public NPS::VDetector{ + //////////////////////////////////////////////////// + /////// Default Constructor and Destructor ///////// + //////////////////////////////////////////////////// +public: + Actar() ; + virtual ~Actar() ; + + //////////////////////////////////////////////////// + /////// Specific Function of this Class /////////// + //////////////////////////////////////////////////// +public: + // Cartesian + void AddDetector(G4ThreeVector POS, string Shape); + // Spherical + void AddDetector(double R,double Theta,double Phi,string Shape); + + + G4LogicalVolume* BuildDetector(); + +private: + G4LogicalVolume* m_SquareDetector; + bool m_build_Silicon; + bool m_build_CsI; + G4LogicalVolume* m_LogicSilicon; + G4LogicalVolume* m_LogicCsICrystal; + + //////////////////////////////////////////////////// + ////// Inherite from NPS::VDetector class ///////// + //////////////////////////////////////////////////// +public: + // Read stream at Configfile to pick-up parameters of detector (Position,...) + // Called in DetecorConstruction::ReadDetextorConfiguration Method + void ReadConfiguration(NPL::InputParser) ; + + // Construct detector and inialise sensitive part. + // Called After DetecorConstruction::AddDetector Method + void ConstructDetector(G4LogicalVolume* world) ; + + // Add Detector branch to the EventTree. + // Called After DetecorConstruction::AddDetector Method + void InitializeRootOutput() ; + + // Read sensitive part and fill the Root tree. + // Called at in the EventAction::EndOfEventAvtion + void ReadSensitive(const G4Event* event) ; + +public: // Scorer + // Initialize all Scorer + void InitializeScorers() ; + void SimulateDigitizer(vector<double> E, vector<double> T, double fallTime,double start,double stop,double step,double noise); + + // Associated Scorer + G4MultiFunctionalDetector* m_ActarScorer ; + G4MultiFunctionalDetector* m_CsIScorer ; + G4MultiFunctionalDetector* m_SiliconScorer ; + //////////////////////////////////////////////////// + ///////////Event class to store Data//////////////// + //////////////////////////////////////////////////// +private: + TActarData* m_Event ; + + //////////////////////////////////////////////////// + ///////////////Private intern Data////////////////// + //////////////////////////////////////////////////// +private: // Geometry + // Detector Coordinate + vector<double> m_R; + vector<double> m_Theta; + vector<double> m_Phi; + + // Shape type + vector<string> m_Shape ; + + // map + map<int, int> m_PadToXRow; + map<int, int> m_PadToZColumn; + + // token + vector<string> m_GasMaterial; + vector<int> m_GasFraction; + double m_Pressure; // bar + double m_Temperature; // kelvin + + // Visualisation Attribute + G4VisAttributes* m_VisChamber; + G4VisAttributes* m_VisWindows; + G4VisAttributes* m_VisGas; + G4VisAttributes* m_VisPads; + G4VisAttributes* m_SiliconVisAtt; + G4VisAttributes* m_CsIVisAtt; + // Needed for dynamic loading of the library +public: + static NPS::VDetector* Construct(); +}; +#endif diff --git a/NPSimulation/Detectors/Actar/CMakeLists.txt b/NPSimulation/Detectors/Actar/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..2e1163d4fb99272a3e21fc19e495963c0a482ff6 --- /dev/null +++ b/NPSimulation/Detectors/Actar/CMakeLists.txt @@ -0,0 +1,2 @@ +add_library(NPSActar SHARED Actar.cc) +target_link_libraries(NPSActar NPSCore ${ROOT_LIBRARIES} ${Geant4_LIBRARIES} ${NPLib_LIBRARIES} -lNPActar) diff --git a/NPSimulation/Detectors/ForwardArray/ForwardArray.cc b/NPSimulation/Detectors/ForwardArray/ForwardArray.cc index 8bb623e61991d72a8c221f88bc9dd21974d91dae..71385a978517517895c473667e28b38887d53430 100644 --- a/NPSimulation/Detectors/ForwardArray/ForwardArray.cc +++ b/NPSimulation/Detectors/ForwardArray/ForwardArray.cc @@ -1,14 +1,14 @@ /***************************************************************************** - * Copyright (C) 2009-2017 this file is part of the NPTool Project * + * Copyright (C) 2009-2017 this file is part of the NPTool Project * * * * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * * For the list of contributors see $NPTOOL/Licence/Contributors * *****************************************************************************/ /***************************************************************************** - * Original Author: Pierre Morfouace contact address: morfouac@nscl.msu.edu * + * Original Author: Pierre Morfouace contact address: morfouac@nscl.msu.edu * * * - * Creation Date : May 2017 * + * Creation Date : May 2017 * * Last update : * *---------------------------------------------------------------------------* * Decription: * @@ -57,7 +57,7 @@ using namespace CLHEP; namespace ForwardArray_NS{ // Energy and time Resolution const double EnergyThreshold = 0.1*MeV; - const double ResoTime = 4.5*ns ; + const double ResoTime = 0.3*ns ; const double ResoEnergy = 1.0*MeV ; //const double Width = 100*mm ; diff --git a/NPSimulation/Detectors/Tiara/Tiara.cc b/NPSimulation/Detectors/Tiara/Tiara.cc index 235a63a3a7e80ed97e415e573565a25198099ae7..f2d240f862e70f08687aab6e913f3d4f33f940fa 100644 --- a/NPSimulation/Detectors/Tiara/Tiara.cc +++ b/NPSimulation/Detectors/Tiara/Tiara.cc @@ -246,24 +246,50 @@ void Tiara::ReadSensitive(const G4Event* event){ double EF = RandGauss::shoot(Info[0],ResoEnergyHyball); if(EF>EnergyThreshold){ int RingNumber=Info[8]; - RingNumber=abs(RingNumber-17); + //by Shuya 171009 + RingNumber=abs(RingNumber); + //RingNumber=abs(RingNumber-17); Info[8]=RingNumber; m_EventHyball->SetRingE(Info[7],Info[8],EF); m_EventHyball->SetRingT(Info[7],Info[8],Info[1]); } // Back Energy - double EB = RandGauss::shoot(Info[1]+Info[0],ResoEnergyHyball); +//by Shuya 171009. Infor[1] is Timing data... no make sense. + //double EB = RandGauss::shoot(Info[1]+Info[0],ResoEnergyHyball); + double EB = RandGauss::shoot(Info[0],ResoEnergyHyball); if(EB>EnergyThreshold){ - m_EventHyball->SetSectorE(Info[7],Info[9],EF); +//by Shuya 171012. + double m_axis = -100.0; + double m_phi = Info[6]; + if(Info[6]<0) m_phi = Info[6]+2.0*M_PI; + if(Info[7]==1) m_axis = 210.0/180.0*M_PI; + else if(Info[7]==2) m_axis = 150.0/180.0*M_PI; + else if(Info[7]==3) m_axis = 90.0/180.0*M_PI; + else if(Info[7]==4) m_axis = 30.0/180.0*M_PI; + else if(Info[7]==5) m_axis = 330.0/180.0*M_PI; + else if(Info[7]==6) m_axis = 270.0/180.0*M_PI; + double m_StripPitchSector_Tiara = (0.5*HYBALL_ActiveWafer_Angle-(-0.5*HYBALL_ActiveWafer_Angle))/HYBALL_NumberOfRadialStrip; + Info[9] = (int)((m_phi - (m_axis - 0.5*HYBALL_ActiveWafer_Angle)) / m_StripPitchSector_Tiara ) + 1 ; +//by Shuya 171009 + //m_EventHyball->SetSectorE(Info[7],Info[9],EF); + m_EventHyball->SetSectorE(Info[7],Info[9],EB); m_EventHyball->SetSectorT(Info[7],Info[9],Info[1]); } // Interaction Coordinates +//by Shuya 171009 +/* ms_InterCoord->SetDetectedPositionX(Info[5]) ; ms_InterCoord->SetDetectedPositionY(Info[6]) ; ms_InterCoord->SetDetectedPositionZ(Info[7]) ; ms_InterCoord->SetDetectedAngleTheta(Info[8]/deg) ; ms_InterCoord->SetDetectedAnglePhi(Info[9]/deg) ; +*/ + ms_InterCoord->SetDetectedPositionX(Info[2]) ; + ms_InterCoord->SetDetectedPositionY(Info[3]) ; + ms_InterCoord->SetDetectedPositionZ(Info[4]) ; + ms_InterCoord->SetDetectedAngleTheta(Info[5]/deg) ; + ms_InterCoord->SetDetectedAnglePhi(Info[6]/deg) ; } // Clear Map for next event HyballHitMap->clear(); diff --git a/NPSimulation/Scorers/DriftElectronScorers.cc b/NPSimulation/Scorers/DriftElectronScorers.cc index 4c0548f425673d72e9cd3bbfe81496dfd4905aa3..2f66eaeece89e23770406a042565ed8307199fdf 100644 --- a/NPSimulation/Scorers/DriftElectronScorers.cc +++ b/NPSimulation/Scorers/DriftElectronScorers.cc @@ -40,11 +40,11 @@ G4bool PS_DECathode::ProcessHits(G4Step* aStep, G4TouchableHistory*){ // contain Energy Time, DetNbr, StripFront and StripBack G4double* Infos = new G4double[9]; Infos[0] = 0; - Infos[1] = aStep->GetPreStepPoint()->GetGlobalTime(); - + Infos[1] = aStep->GetPreStepPoint()->GetProperTime(); + m_DetectorNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(m_Level); m_Position = aStep->GetPreStepPoint()->GetPosition(); - + // Interaction coordinates (used to fill the InteractionCoordinates branch) Infos[2] = m_Position.x(); Infos[3] = m_Position.y(); @@ -52,10 +52,10 @@ G4bool PS_DECathode::ProcessHits(G4Step* aStep, G4TouchableHistory*){ Infos[5] = m_Position.theta(); Infos[6] = m_Position.phi(); Infos[7] = m_DetectorNumber; - + m_Index = m_DetectorNumber * 1e3 ; G4String PID = aStep->GetTrack()->GetDefinition()->GetParticleName(); - + if(PID=="driftelectron"){ Infos[0] = 1; } @@ -68,7 +68,7 @@ G4bool PS_DECathode::ProcessHits(G4Step* aStep, G4TouchableHistory*){ Infos[0]+=dummy[0]; Infos[1]=dummy[1]; } - + EvtMap->set(m_Index, Infos); return TRUE; } @@ -92,13 +92,13 @@ void PS_DECathode::clear(){ for (MapIterator = EvtMap->GetMap()->begin() ; MapIterator != EvtMap->GetMap()->end() ; MapIterator++){ delete *(MapIterator->second); } - + EvtMap->clear(); } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void PS_DECathode::DrawAll(){ - + } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... @@ -126,11 +126,11 @@ G4bool PS_DEDigitizer::ProcessHits(G4Step* aStep, G4TouchableHistory*){ // contain Energy Time, DetNbr, StripFront and StripBack G4double* Infos = new G4double[9]; Infos[0] = 0; - Infos[1] = aStep->GetPreStepPoint()->GetGlobalTime(); - + Infos[1] = aStep->GetPreStepPoint()->GetProperTime(); + m_DetectorNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(m_Level); m_Position = aStep->GetPreStepPoint()->GetPosition(); - + // Interaction coordinates (used to fill the InteractionCoordinates branch) Infos[2] = m_Position.x(); Infos[3] = m_Position.y(); @@ -138,10 +138,10 @@ G4bool PS_DEDigitizer::ProcessHits(G4Step* aStep, G4TouchableHistory*){ Infos[5] = m_Position.theta(); Infos[6] = m_Position.phi(); Infos[7] = m_DetectorNumber; - + m_Index = aStep->GetTrack()->GetTrackID() + m_DetectorNumber*1e6 ; G4String PID = aStep->GetTrack()->GetDefinition()->GetParticleName(); - + if(PID=="driftelectron"){ Infos[0] = 1; } @@ -154,7 +154,7 @@ G4bool PS_DEDigitizer::ProcessHits(G4Step* aStep, G4TouchableHistory*){ Infos[0]+=dummy[0]; Infos[1]=dummy[1]; } - + EvtMap->set(m_Index, Infos); return TRUE; } @@ -178,13 +178,13 @@ void PS_DEDigitizer::clear(){ for (MapIterator = EvtMap->GetMap()->begin() ; MapIterator != EvtMap->GetMap()->end() ; MapIterator++){ delete *(MapIterator->second); } - + EvtMap->clear(); } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void PS_DEDigitizer::DrawAll(){ - + } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... @@ -194,5 +194,3 @@ void PS_DEDigitizer::PrintAll(){ G4cout << " Number of entries " << EvtMap->entries() << G4endl ; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... - - diff --git a/NPSimulation/Scorers/SiliconScorers.cc b/NPSimulation/Scorers/SiliconScorers.cc index 56fce8e262f14f9ba846e21245a2f3335653460b..b9674a66c3e751dcee5679795697b50d444ed4c7 100644 --- a/NPSimulation/Scorers/SiliconScorers.cc +++ b/NPSimulation/Scorers/SiliconScorers.cc @@ -338,8 +338,7 @@ G4bool PS_Silicon_Annular::ProcessHits(G4Step* aStep, G4TouchableHistory*){ m_StripRingNumber = (int) ((m_Position.rho() - m_StripPlaneInnerRadius) / m_StripPitchRing ) + 1 ; double phi = m_Position.phi(); - if(phi<0) - phi+=2*M_PI; + m_StripSectorNumber = (int) ((phi - m_PhiStart) / m_StripPitchSector ) + 1 ; m_StripQuadrantNumber = (int) ((phi - m_PhiStart) /m_StripPitchQuadrant) + 1 ; diff --git a/Projects/Actar/Analysis.cxx b/Projects/Actar/Analysis.cxx new file mode 100644 index 0000000000000000000000000000000000000000..23c45664f868629cba52cf075ca4231dab73cd85 --- /dev/null +++ b/Projects/Actar/Analysis.cxx @@ -0,0 +1,71 @@ +/***************************************************************************** + * Copyright (C) 2009-2016 this file is part of the NPTool Project * + * * + * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * + * For the list of contributors see $NPTOOL/Licence/Contributors * + *****************************************************************************/ + +/***************************************************************************** + * Original Author: XAUTHORX contact address: XMAILX * + * * + * Creation Date : XMONTHX XYEARX * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class describe Actar analysis project * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + *****************************************************************************/ + +#include<iostream> +using namespace std; +#include"Analysis.h" +#include"NPAnalysisFactory.h" +#include"NPDetectorManager.h" +//////////////////////////////////////////////////////////////////////////////// +Analysis::Analysis(){ +} +//////////////////////////////////////////////////////////////////////////////// +Analysis::~Analysis(){ +} + +//////////////////////////////////////////////////////////////////////////////// +void Analysis::Init(){ + Actar= (TActarPhysics*) m_DetectorManager->GetDetector("Actar"); +} + +//////////////////////////////////////////////////////////////////////////////// +void Analysis::TreatEvent(){ + //cout << "/// Size=" << Actar->PadRow.size() << endl; + for(unsigned int i=0; i<Actar->PadRow.size(); i++){ + //cout << "Row= " << Actar->PadRow[i] << endl; + } +} + +//////////////////////////////////////////////////////////////////////////////// +void Analysis::End(){ +} + + +//////////////////////////////////////////////////////////////////////////////// +// Construct Method to be pass to the DetectorFactory // +//////////////////////////////////////////////////////////////////////////////// +NPL::VAnalysis* Analysis::Construct(){ + return (NPL::VAnalysis*) new Analysis(); +} + +//////////////////////////////////////////////////////////////////////////////// +// Registering the construct method to the factory // +//////////////////////////////////////////////////////////////////////////////// +extern "C"{ +class proxy{ + public: + proxy(){ + NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct); + } +}; + +proxy p; +} diff --git a/Projects/Actar/Analysis.h b/Projects/Actar/Analysis.h new file mode 100644 index 0000000000000000000000000000000000000000..075b9abe68fd6c3255f58bd01f22acbddee18ffd --- /dev/null +++ b/Projects/Actar/Analysis.h @@ -0,0 +1,42 @@ +#ifndef Analysis_h +#define Analysis_h +/***************************************************************************** + * Copyright (C) 2009-2016 this file is part of the NPTool Project * + * * + * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * + * For the list of contributors see $NPTOOL/Licence/Contributors * + *****************************************************************************/ + +/***************************************************************************** + * Original Author: XAUTHORX contact address: XMAILX * + * * + * Creation Date : XMONTHX XYEARX * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class describe Actar analysis project * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + *****************************************************************************/ + +#include"NPVAnalysis.h" +#include"TActarPhysics.h" +class Analysis: public NPL::VAnalysis{ + public: + Analysis(); + ~Analysis(); + + public: + void Init(); + void TreatEvent(); + void End(); + + static NPL::VAnalysis* Construct(); + + private: + TActarPhysics* Actar; + +}; +#endif diff --git a/Projects/Actar/CMakeLists.txt b/Projects/Actar/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..22c74affdfc45019bdda2594f8439c52d4ab97ec --- /dev/null +++ b/Projects/Actar/CMakeLists.txt @@ -0,0 +1,5 @@ +cmake_minimum_required (VERSION 2.8) +# Setting the policy to match Cmake version +cmake_policy(VERSION ${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION}) +# include the default NPAnalysis cmake file +include("../../NPLib/ressources/CMake/NPAnalysis.cmake") diff --git a/Projects/HiraU_experiment/.ls_return b/Projects/HiraU_experiment/.ls_return index 6eb02b2613bfaf2946d0742a93989fef82426bc5..c40ba7451e09136b0f2db068aa7b0e9fb4d8552e 100644 --- a/Projects/HiraU_experiment/.ls_return +++ b/Projects/HiraU_experiment/.ls_return @@ -1 +1 @@ -../../Outputs/Simulation/ca40ni58_e140_b2.root +../../Outputs/Simulation/ca40ni58_e50_b2.root diff --git a/Projects/HiraU_experiment/Analysis.cxx b/Projects/HiraU_experiment/Analysis.cxx index 8b49c4d2cfb6fdcc84c99ee7a5b01c087ba663fe..a09b52fa18bab8b31e52ae506000ffe16e5bfe2f 100644 --- a/Projects/HiraU_experiment/Analysis.cxx +++ b/Projects/HiraU_experiment/Analysis.cxx @@ -130,6 +130,7 @@ void Analysis::TreatEvent(){ for(int i=0; i<FA->Energy.size(); i++){ ThetaLabFA = InteractionCoordinates->GetDetectedAngleTheta(i); PhiLabFA = InteractionCoordinates->GetDetectedAnglePhi(i); + StartTimeFA = FA->Time[0]; } FAMultiplicity = FA->DetectorNumber.size(); @@ -163,6 +164,7 @@ void Analysis::InitOutputBranch() { RootOutput::getInstance()->GetTree()->Branch("ThetaLabFA",&ThetaLabFA,"ThetaLabFA/D"); RootOutput::getInstance()->GetTree()->Branch("PhiLabFA",&PhiLabFA,"PhiLabFA/D"); + RootOutput::getInstance()->GetTree()->Branch("StartTimeFA",&StartTimeFA,"StartTimeFA/D"); RootOutput::getInstance()->GetTree()->Branch("MBMultiplicity",&MBMultiplicity,"MBMultiplicity/I"); RootOutput::getInstance()->GetTree()->Branch("NWMultiplicity",&NWMultiplicity,"NWMultiplicity/I"); @@ -200,6 +202,7 @@ void Analysis::ReInitValue(){ EF = -100; HiraELab = -100; E_CsI = -100; + StartTimeFA = -100; FAMultiplicity = -1; NWMultiplicity = -1; diff --git a/Projects/HiraU_experiment/Analysis.h b/Projects/HiraU_experiment/Analysis.h index 545dc165f32a138a4d1472e3904e615377db45b2..67544bd735501b35944dcc7362c2756a1a78342b 100644 --- a/Projects/HiraU_experiment/Analysis.h +++ b/Projects/HiraU_experiment/Analysis.h @@ -59,6 +59,7 @@ private: int MBMultiplicity; int NWMultiplicity; int FAMultiplicity; + double StartTimeFA; double X_Hira; double Y_Hira; double Z_Hira; diff --git a/Projects/HiraU_experiment/RunToTreat.txt b/Projects/HiraU_experiment/RunToTreat.txt index 9a91e98a01774ba6e72dc27229314702cf7f85b8..40acab93914e12c670b25c79f80104d3d55b4be3 100755 --- a/Projects/HiraU_experiment/RunToTreat.txt +++ b/Projects/HiraU_experiment/RunToTreat.txt @@ -1,6 +1,6 @@ TTreeName SimulatedTree RootFileName - ../../Outputs/Simulation/ca40ni58_e140_b2.root + ../../Outputs/Simulation/ca40ni58_e50_b2.root %../../Outputs/Simulation/hiraU_12Cpp.root %../../Outputs/Simulation/hiraU_1Hpp.root diff --git a/Projects/T40/Analysis.cxx b/Projects/T40/Analysis.cxx index 9078e40695ffd876e78dd4245fcbb8758b7d4922..8d0e271847461b58cbc46c7a206614576f941030 100644 --- a/Projects/T40/Analysis.cxx +++ b/Projects/T40/Analysis.cxx @@ -36,6 +36,38 @@ using namespace std; #include "TFitResult.h" #include "TFitResultPtr.h" +////// INTERNAL FUNCTIONS ////////// +namespace { +double calculate_fit_slope(int len, double* Aw_X, double* Aw_Z, double& R2) +{ + vector<double> X, Z; + for(int i=0; i< len; ++i) { + if(Aw_X[i] != -1000) { X.push_back(Aw_X[i]); } + if(Aw_Z[i] != -1000) { Z.push_back(Aw_Z[i]); } + } + + Long64_t N = X.size(); + double meanZ = TMath::Mean(N, &Z[0]); + double meanX = TMath::Mean(N, &X[0]); + double meanZ2 = 0, meanXZ = 0, meanX2 = 0; + for(size_t i=0; i< N; ++i) { + meanZ2 += Z[i]*Z[i]; + meanX2 += X[i]*X[i]; + meanXZ += Z[i]*X[i]; + } + meanZ2 /= N; + meanXZ /= N; + + double slope = (meanXZ - meanX*meanZ) / (meanZ2 - meanZ*meanZ); + R2 = pow(meanXZ - meanX*meanZ, 2) / + ((meanZ2 - meanZ*meanZ) * (meanX2 - meanX*meanX)); + + /// TODO::: R2 doesn't seem to make sense... look into it! + + return slope; +} } + + //////////////////////////////////////////////////////////////////////////////// Analysis::Analysis(){ } @@ -59,7 +91,7 @@ void Analysis::Init(){ cout << "Original Beam energy (entrance of target): " << OriginalBeamEnergy << endl ; // target thickness - TargetThickness = 0*m_DetectorManager->GetTargetThickness(); + TargetThickness = m_DetectorManager->GetTargetThickness(); string TargetMaterial = m_DetectorManager->GetTargetMaterial(); // energy losses @@ -86,6 +118,10 @@ void Analysis::Init(){ //LightAl = NPL::EnergyLoss("He4_Al.SRIM","SRIM",10); LightSi = NPL::EnergyLoss(light+"_Si.SRIM","SRIM",10); //LightSi = NPL::EnergyLoss("He4_Si.SRIM","SRIM",10); + +//by Shuya 170530 + //LightCBacking = NPL::EnergyLoss(light+"_C.SRIM","SRIM",10); + BeamTarget = NPL::EnergyLoss(beam+"_"+TargetMaterial+".SRIM","SRIM",10); FinalBeamEnergy = BeamTarget.Slow(OriginalBeamEnergy, TargetThickness*0.5, 0); myReaction->SetBeamEnergy(FinalBeamEnergy); @@ -120,7 +156,9 @@ void Analysis::Init(){ Micro1_E_row1_2 = 0; // Energy from micromega rows 1 & 2 ("delta E in stopping mode") Micro2_E_row1_2 = 0; // Energy from micromega rows 1-2 ("E in stopping mode") Micro1_E_row1 = 0 ;// Energy from micromega row 1 - Micro1_E_col4 = 0 ;// energy from micromega col 4 +//by Shuya 170912 + //Micro1_E_col4 = 0 ;// energy from micromega col 4 + Micro1_E_col4_sum = 0 ;// energy from micromega col 4 Plast_E = 0; // Energy Plastic for(int i=0; i< kNumAw; ++i) { Aw_X[i] = -1000; @@ -130,19 +168,35 @@ void Analysis::Init(){ Aw_ThetaFit = -1000; Aw_ThetaFit_R2 = -1000; //by Shuya 170516 - Micro1_E_col1 = 0. ;// energy from micromega col 1 - Micro1_E_col2 = 0. ;// energy from micromega col 2 - Micro1_E_col3 = 0. ;// energy from micromega col 3 - Micro1_E_col5 = 0. ;// energy from micromega col 5 - Micro1_E_col6 = 0. ;// energy from micromega col 6 - Micro1_E_col7 = 0. ;// energy from micromega col 7 - Micro2_E_col1 = 0. ;// energy from micromega2 col 1 - Micro2_E_col2 = 0. ;// energy from micromega2 col 2 - Micro2_E_col3 = 0. ;// energy from micromega2 col 3 - Micro2_E_col4 = 0. ;// energy from micromega2 col 3 - Micro2_E_col5 = 0. ;// energy from micromega2 col 5 - Micro2_E_col6 = 0. ;// energy from micromega2 col 6 - Micro2_E_col7 = 0. ;// energy from micromega2 col 7 +//by Shuya 170912 + Micro1_E_col1_sum = 0. ;// energy from micromega col 1 + Micro1_E_col2_sum = 0. ;// energy from micromega col 2 + Micro1_E_col3_sum = 0. ;// energy from micromega col 3 + Micro1_E_col5_sum = 0. ;// energy from micromega col 5 + Micro1_E_col6_sum = 0. ;// energy from micromega col 6 + Micro1_E_col7_sum = 0. ;// energy from micromega col 7 + Micro2_E_col1_sum = 0. ;// energy from micromega2 col 1 + Micro2_E_col2_sum = 0. ;// energy from micromega2 col 2 + Micro2_E_col3_sum = 0. ;// energy from micromega2 col 3 + Micro2_E_col4_sum = 0. ;// energy from micromega2 col 3 + Micro2_E_col5_sum = 0. ;// energy from micromega2 col 5 + Micro2_E_col6_sum = 0. ;// energy from micromega2 col 6 + Micro2_E_col7_sum = 0. ;// energy from micromega2 col 7 +//by Shuya 170912 + Micro1_E_col1_mult = 0. ;// energy from micromega col 1 + Micro1_E_col2_mult = 0. ;// energy from micromega col 2 + Micro1_E_col3_mult = 0. ;// energy from micromega col 3 + Micro1_E_col4_mult = 0. ;// energy from micromega col 3 + Micro1_E_col5_mult = 0. ;// energy from micromega col 5 + Micro1_E_col6_mult = 0. ;// energy from micromega col 6 + Micro1_E_col7_mult = 0. ;// energy from micromega col 7 + Micro2_E_col1_mult = 0. ;// energy from micromega2 col 1 + Micro2_E_col2_mult = 0. ;// energy from micromega2 col 2 + Micro2_E_col3_mult = 0. ;// energy from micromega2 col 3 + Micro2_E_col4_mult = 0. ;// energy from micromega2 col 3 + Micro2_E_col5_mult = 0. ;// energy from micromega2 col 5 + Micro2_E_col6_mult = 0. ;// energy from micromega2 col 6 + Micro2_E_col7_mult = 0. ;// energy from micromega2 col 7 //TAC TacSiGeOR = -1000; @@ -175,11 +229,16 @@ void Analysis::TreatEvent(){ ThetaNormalTarget = 0; if(XTarget>-1000 && YTarget>-1000){ TVector3 BeamImpact(XTarget,YTarget,0); + TVector3 HitDirection = TH -> GetRandomisedPositionOfInteraction(countTiaraHyball) - BeamImpact ; ThetaLab = HitDirection.Angle( BeamDirection ); ThetaTHSurface = HitDirection.Angle(TVector3(0,0,-1)); // vector Normal on Hyball ThetaNormalTarget = HitDirection.Angle( TVector3(0,0,1) ) ; + + //by Shuya 171019 + PhiLab = HitDirection.Phi(); + PhiLab = PhiLab/(TMath::Pi())*180.0; } else{ BeamDirection = TVector3(-1000,-1000,-1000); @@ -193,7 +252,9 @@ void Analysis::TreatEvent(){ Si_E_TH = TH->Strip_E[countTiaraHyball]; Energy = Si_E_TH; // calibration for hyball is in MeV // Correct for energy loss using the thickness of the target and the dead layer - ELab = LightSi.EvaluateInitialEnergy( Energy ,0.61*micrometer , ThetaTHSurface); // 0.1 um of Aluminum + ELab = LightSi.EvaluateInitialEnergy( Energy ,0.61*micrometer , ThetaTHSurface); // equivalent to 0.1 um of Aluminum +//by Shuya 170530 + //if(ThetaNormalTarget < halfpi) ELab = LightCBacking.EvaluateInitialEnergy( ELab ,0.044*micrometer , ThetaNormalTarget); //10 ug/cm2 carbon ELab = LightTarget.EvaluateInitialEnergy( ELab ,TargetThickness/2., ThetaNormalTarget); ///////////////////////////// @@ -203,6 +264,9 @@ void Analysis::TreatEvent(){ ThetaCM = myReaction -> EnergyLabToThetaCM( ELab , ThetaLab)/deg; ThetaLab=ThetaLab/deg; +//by Shuya 170703 + Ex_Hyball = Ex; + ///////////////////////////// // Part 5 : Implementing randomised position impact matrix for the Hyball TVector3 HyballRandomImpactPosition = TH -> GetRandomisedPositionOfInteraction(countTiaraHyball); @@ -220,6 +284,7 @@ void Analysis::TreatEvent(){ ThetaNormalTarget = 0; if(XTarget>-1000 && YTarget>-1000){ TVector3 BeamImpact(XTarget,YTarget,0); + TVector3 HitDirection = TB -> GetRandomisedPositionOfInteraction(countTiaraBarrel) - BeamImpact ; //Angle of emission wrt to beam ThetaLab = HitDirection.Angle( BeamDirection ); @@ -229,6 +294,10 @@ void Analysis::TreatEvent(){ int det = TB->Detector_N[countTiaraBarrel]; NormalOnBarrel.RotateZ((3-det)*45*deg); ThetaTBSurface = HitDirection.Angle(NormalOnBarrel); + + //by Shuya 171019 + PhiLab = HitDirection.Phi(); + PhiLab = PhiLab/(TMath::Pi())*180.0; } else{ BeamDirection = TVector3(-1000,-1000,-1000); @@ -251,11 +320,13 @@ void Analysis::TreatEvent(){ // Evaluate energy using the thickness, Target and Si dead layer Correction ELab = LightSi.EvaluateInitialEnergy( Energy ,0.3*micrometer, ThetaTBSurface); - //ELab = LightTarget.EvaluateInitialEnergy( ELab ,TargetThickness/2., ThetaNormalTarget); + ELab = LightTarget.EvaluateInitialEnergy( ELab ,TargetThickness/2., ThetaNormalTarget); ///////////////////////////// // Part 3 : Excitation Energy Calculation Ex = myReaction -> ReconstructRelativistic( ELab , ThetaLab ); +//by Shuya 170703 + Ex_Barrel = Ex; ////////////////////////////// // Part 4 : Theta CM Calculation @@ -301,7 +372,8 @@ void Analysis::TreatEvent(){ if(TF->MicroRowNumber.size()) { Micro1_E_row1 = TF->GetMicroGroupEnergy(1,1,1,1,7); // energy sum from the row 1 - Micro1_E_col4 = TF->GetMicroGroupEnergy(1,1,4,4,4); // energy sum from the col 4 + //by Shuya 170912 + Micro1_E_col4_sum = TF->GetMicroGroupEnergy(1,1,4,4,4); // energy sum from the col 4 //by Shuya 170516. Since Micro2_E is dE detector, we always expect the particles penetrate through the whole rows. That is why you should use GetMicroRowGeomEnergy() instead of GetMicroGroupEnergy(). //Micro2_E = TF->GetMicroGroupEnergy(2,1,4,1,7); // energy sum from all the pads Micro2_E = TF->GetMicroRowGeomEnergy(2,1,4,0); // energy sum from all the pads @@ -309,42 +381,76 @@ void Analysis::TreatEvent(){ Micro2_E_row1_2 = TF->GetMicroGroupEnergy(2,1,2,1,7); // energy sum from row 3-6 //by Shuya 170516. For Micro1's energy sum, you need to choose which of GetMicroGroupEnergy() and GetMicroRowGeomEnergy(). If you're using the particles penetrate the Micro1, better to use GetMicroRowGeomEnergy(). - Micro1_E_col1 = TF->GetMicroGroupEnergy(1,1,4,1,1); // energy sum from the col 1 - Micro1_E_col2 = TF->GetMicroGroupEnergy(1,1,4,2,2); // energy sum from the col 2 - Micro1_E_col3 = TF->GetMicroGroupEnergy(1,1,4,3,3); // energy sum from the col 3 - Micro1_E_col5 = TF->GetMicroGroupEnergy(1,1,4,5,5); // energy sum from the col 5 - Micro1_E_col6 = TF->GetMicroGroupEnergy(1,1,4,6,6); // energy sum from the col 6 - Micro1_E_col7 = TF->GetMicroGroupEnergy(1,1,4,7,7); // energy sum from the col 7 - - Micro2_E_col1 = TF->GetMicroRowGeomEnergy(2,1,4,1); // energy sum from the col1. - Micro2_E_col2 = TF->GetMicroRowGeomEnergy(2,1,4,2); // energy sum from the col2. - Micro2_E_col3 = TF->GetMicroRowGeomEnergy(2,1,4,3); // energy sum from the col3. - Micro2_E_col4 = TF->GetMicroRowGeomEnergy(2,1,4,4); // energy sum from the col4. - Micro2_E_col5 = TF->GetMicroRowGeomEnergy(2,1,4,5); // energy sum from the col5. - Micro2_E_col6 = TF->GetMicroRowGeomEnergy(2,1,4,6); // energy sum from the col6. - Micro2_E_col7 = TF->GetMicroRowGeomEnergy(2,1,4,7); // energy sum from the col7. + //by Shuya 170912 + Micro1_E_col1_sum = TF->GetMicroGroupEnergy(1,1,4,1,1); // energy sum from the col 1 + Micro1_E_col2_sum = TF->GetMicroGroupEnergy(1,1,4,2,2); // energy sum from the col 2 + Micro1_E_col3_sum = TF->GetMicroGroupEnergy(1,1,4,3,3); // energy sum from the col 3 + Micro1_E_col5_sum = TF->GetMicroGroupEnergy(1,1,4,5,5); // energy sum from the col 5 + Micro1_E_col6_sum = TF->GetMicroGroupEnergy(1,1,4,6,6); // energy sum from the col 6 + Micro1_E_col7_sum = TF->GetMicroGroupEnergy(1,1,4,7,7); // energy sum from the col 7 + Micro2_E_col1_mult = TF->GetMicroRowGeomEnergy(2,1,4,1); // energy sum from the col1. + Micro2_E_col2_mult = TF->GetMicroRowGeomEnergy(2,1,4,2); // energy sum from the col2. + Micro2_E_col3_mult = TF->GetMicroRowGeomEnergy(2,1,4,3); // energy sum from the col3. + Micro2_E_col4_mult = TF->GetMicroRowGeomEnergy(2,1,4,4); // energy sum from the col4. + Micro2_E_col5_mult = TF->GetMicroRowGeomEnergy(2,1,4,5); // energy sum from the col5. + Micro2_E_col6_mult = TF->GetMicroRowGeomEnergy(2,1,4,6); // energy sum from the col6. + Micro2_E_col7_mult = TF->GetMicroRowGeomEnergy(2,1,4,7); // energy sum from the col7. + + //by Shuya 170912 + Micro1_E_col1_mult = TF->GetMicroRowGeomEnergy(1,1,4,1); // energy sum from the col1. + Micro1_E_col2_mult = TF->GetMicroRowGeomEnergy(1,1,4,2); // energy sum from the col2. + Micro1_E_col3_mult = TF->GetMicroRowGeomEnergy(1,1,4,3); // energy sum from the col3. + Micro1_E_col4_mult = TF->GetMicroRowGeomEnergy(1,1,4,4); // energy sum from the col4. + Micro1_E_col5_mult = TF->GetMicroRowGeomEnergy(1,1,4,5); // energy sum from the col5. + Micro1_E_col6_mult = TF->GetMicroRowGeomEnergy(1,1,4,6); // energy sum from the col6. + Micro1_E_col7_mult = TF->GetMicroRowGeomEnergy(1,1,4,7); // energy sum from the col7. + Micro2_E_col1_sum = TF->GetMicroGroupEnergy(2,1,4,1,1); // energy sum from the col 1 + Micro2_E_col2_sum = TF->GetMicroGroupEnergy(2,1,4,2,2); // energy sum from the col 2 + Micro2_E_col3_sum = TF->GetMicroGroupEnergy(2,1,4,3,3); // energy sum from the col 3 + Micro2_E_col4_sum = TF->GetMicroGroupEnergy(2,1,4,4,4); // energy sum from the col 3 + Micro2_E_col5_sum = TF->GetMicroGroupEnergy(2,1,4,5,5); // energy sum from the col 5 + Micro2_E_col6_sum = TF->GetMicroGroupEnergy(2,1,4,6,6); // energy sum from the col 6 + Micro2_E_col7_sum = TF->GetMicroGroupEnergy(2,1,4,7,7); // energy sum from the col 7 } else { Micro1_E_row1 = -1000; - Micro1_E_col4 = -1000; + //by Shuya 170912 + Micro1_E_col4_sum = -1000; Micro1_E_row1_2 = -1000; Micro2_E_row1_2 = -1000; Micro2_E = -1000; //by Shuya 170516 - Micro1_E_col1 = -1000; - Micro1_E_col2 = -1000; - Micro1_E_col3 = -1000; - Micro1_E_col5 = -1000; - Micro1_E_col6 = -1000; - Micro1_E_col7 = -1000; - Micro2_E_col1 = -1000; - Micro2_E_col2 = -1000; - Micro2_E_col3 = -1000; - Micro2_E_col4 = -1000; - Micro2_E_col5 = -1000; - Micro2_E_col6 = -1000; - Micro2_E_col7 = -1000; + //by Shuya 170912 + Micro1_E_col1_sum = -1000; + Micro1_E_col2_sum = -1000; + Micro1_E_col3_sum = -1000; + Micro1_E_col5_sum = -1000; + Micro1_E_col6_sum = -1000; + Micro1_E_col7_sum = -1000; + Micro2_E_col1_sum = -1000; + Micro2_E_col2_sum = -1000; + Micro2_E_col3_sum = -1000; + Micro2_E_col4_sum = -1000; + Micro2_E_col5_sum = -1000; + Micro2_E_col6_sum = -1000; + Micro2_E_col7_sum = -1000; + + //by Shuya 170912 + Micro1_E_col1_mult = -1000; + Micro1_E_col2_mult = -1000; + Micro1_E_col3_mult = -1000; + Micro1_E_col4_mult = -1000; + Micro1_E_col5_mult = -1000; + Micro1_E_col6_mult = -1000; + Micro1_E_col7_mult = -1000; + Micro2_E_col1_mult = -1000; + Micro2_E_col2_mult = -1000; + Micro2_E_col3_mult = -1000; + Micro2_E_col4_mult = -1000; + Micro2_E_col5_mult = -1000; + Micro2_E_col6_mult = -1000; + Micro2_E_col7_mult = -1000; } // Delta E ion chamber @@ -376,9 +482,10 @@ void Analysis::TreatEvent(){ for(int i=0; i< kNumAw; ++i) { if(Aw_X[i] != -1000) { ++numValid; } if(numValid == 2) { // at least 2 points to calculate an angle - Aw_ThetaFit = TF->AWireAngle*(180/TMath::Pi()); - Aw_ThetaFit_R2 = TF->AWireFitR2; + Aw_ThetaFit = TF->AWireAngle*(180/TMath::Pi()); + Aw_ThetaFit_R2 = TF->AWireFitR2; break; + } } @@ -393,11 +500,16 @@ void Analysis::TreatEvent(){ TacSiMicro = TF->MicroTimeOR[0]; for(UInt_t ti = 0; ti< TF->MicroTimeOR.size(); ++ti) { + //by Shuya 170905 - uncomment the second one and comment out the first one if you want to have a MicroMegas_dE Timing data in an appropriate Tree (not PlastLeftTime). + //Note with this, you can't get timing data in neither TacSiMicro_dE or TacSiMicro_E. Only TacSiMicro take the data equivalent to TacSiMicro_E, and TacSiMicro_dE goes to PlastLeftTime due to the setting in t40.txt configuration file. switch(TF->MicroTimeRowNumber[ti]) { + //switch(TF->MicroTimeDetNumber[ti]) { case 0: + //case 2: TacSiMicro_dE = TF->MicroTimeOR[ti]; break; case 5: + //case 1: TacSiMicro_E = TF->MicroTimeOR[ti]; break; default: @@ -406,6 +518,10 @@ void Analysis::TreatEvent(){ } } +//by Shuya 170912 + if(TacSiMicro_dE<-1000 || TacSiMicro_dE>5000.0) TacSiMicro_dE=-500.0; + if(TacSiMicro_E<-1000 || TacSiMicro_E>5000.0) TacSiMicro_E=-500.0; + // For the plastic there's two ways to calculate the times, both ar OR. // The two available time channels i.e. Plast Right and Plast left are used in this case if(TF->PlastRightTime.size()==1) @@ -413,11 +529,17 @@ void Analysis::TreatEvent(){ else TacSiPlastRight = -999; +//by Shuya 170906 + if(TacSiPlastRight<-999.0 || TacSiPlastRight>5000.0) TacSiPlastRight=-500.0; + if(TF->PlastLeftTime.size()==1) TacSiPlastLeft = TF->PlastLeftTime[0]; else TacSiPlastLeft = -999; +//by Shuya 170906 + if(TacSiPlastLeft<-999.0 || TacSiPlastLeft>5000.0) TacSiPlastLeft=-500.0; + if(TG->GeTime.size()==1) TacSiGeOR = TG->GeTime[0]; @@ -451,6 +573,11 @@ void Analysis::ReInitValue(){ ThetaLab = -1000; ThetaCM = -1000; LightParticleDetected = false ; +//by Shuya 170703 + Ex_Hyball = -1000 ; + Ex_Barrel = -1000 ; +//by Shuya 171019 + PhiLab = -1000; //Simu //Original_ELab = -1000; @@ -459,25 +586,43 @@ void Analysis::ReInitValue(){ //FPD Delta_E = -1000; Micro1_E_row1 = -1000; - Micro1_E_col4 = -1000; +//by Shuya 170912 + Micro1_E_col4_sum = -1000; Micro1_E_row1_2 = -1000; Micro2_E_row1_2 = -1000; Micro2_E = -1000; Plast_E = -1000; //by Shuya 170516 - Micro1_E_col1 = -1000; - Micro1_E_col2 = -1000; - Micro1_E_col3 = -1000; - Micro1_E_col5 = -1000; - Micro1_E_col6 = -1000; - Micro1_E_col7 = -1000; - Micro2_E_col1 = -1000; - Micro2_E_col2 = -1000; - Micro2_E_col3 = -1000; - Micro2_E_col4 = -1000; - Micro2_E_col5 = -1000; - Micro2_E_col6 = -1000; - Micro2_E_col7 = -1000; +//by Shuya 170912 + Micro1_E_col1_sum = -1000; + Micro1_E_col2_sum = -1000; + Micro1_E_col3_sum = -1000; + Micro1_E_col5_sum = -1000; + Micro1_E_col6_sum = -1000; + Micro1_E_col7_sum = -1000; + Micro2_E_col1_sum = -1000; + Micro2_E_col2_sum = -1000; + Micro2_E_col3_sum = -1000; + Micro2_E_col4_sum = -1000; + Micro2_E_col5_sum = -1000; + Micro2_E_col6_sum = -1000; + Micro2_E_col7_sum = -1000; + + Micro1_E_col1_mult = -1000; + Micro1_E_col2_mult = -1000; + Micro1_E_col3_mult = -1000; + Micro1_E_col4_mult = -1000; + Micro1_E_col5_mult = -1000; + Micro1_E_col6_mult = -1000; + Micro1_E_col7_mult = -1000; + Micro2_E_col1_mult = -1000; + Micro2_E_col2_mult = -1000; + Micro2_E_col3_mult = -1000; + Micro2_E_col4_mult = -1000; + Micro2_E_col5_mult = -1000; + Micro2_E_col6_mult = -1000; + Micro2_E_col7_mult = -1000; + for(int i=0; i< kNumAw; ++i) { Aw_X[i] = -1000; @@ -505,9 +650,16 @@ void Analysis::ReInitValue(){ void Analysis::InitOutputBranch() { //Tiara RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex,"Ex/D"); +//by Shuya 170703 + RootOutput::getInstance()->GetTree()->Branch("Ex_Hyball",&Ex_Hyball,"Ex_Hyball/D"); + RootOutput::getInstance()->GetTree()->Branch("Ex_Barrel",&Ex_Barrel,"Ex_Barrel/D"); + RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab,"ELab/D"); RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab,"ThetaLab/D"); RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM,"ThetaCM/D"); +//by Shuya 171019 + RootOutput::getInstance()->GetTree()->Branch("PhiLab",&PhiLab,"PhiLab/D"); + RootOutput::getInstance()->GetTree()->Branch("TiaraImpactMatrixX",&TiaraIMX,"TiaraImpactMatrixX/D"); RootOutput::getInstance()->GetTree()->Branch("TiaraImpactMatrixY",&TiaraIMY,"TiaraImpactMatrixY/D"); RootOutput::getInstance()->GetTree()->Branch("TiaraImpactMatrixZ",&TiaraIMZ,"TiaraImpactMatrixZ/D"); @@ -520,7 +672,6 @@ void Analysis::InitOutputBranch() { //FPD RootOutput::getInstance()->GetTree()->Branch("Delta_E",&Delta_E,"Delta_E/D"); RootOutput::getInstance()->GetTree()->Branch("Micro1_E_row1",&Micro1_E_row1,"Micro1_E_row1/D"); - RootOutput::getInstance()->GetTree()->Branch("Micro1_E_col4",&Micro1_E_col4,"Micro1_E_col4/D"); RootOutput::getInstance()->GetTree()->Branch("Micro1_E_row1_2", &Micro1_E_row1_2, "Micro1_E_row1_2/D"); RootOutput::getInstance()->GetTree()->Branch("Micro2_E_row1_2", &Micro2_E_row1_2, "Micro2_E_row1_2/D"); RootOutput::getInstance()->GetTree()->Branch("Micro2_E",&Micro2_E,"Micro2_E/D"); @@ -531,32 +682,54 @@ void Analysis::InitOutputBranch() { RootOutput::getInstance()->GetTree()->Branch("Aw_ThetaFit",&Aw_ThetaFit,"Aw_ThetaFit/D"); RootOutput::getInstance()->GetTree()->Branch("Aw_ThetaFit_R2",&Aw_ThetaFit_R2,"Aw_ThetaFit_R2/D"); //by Shuya 170516 - RootOutput::getInstance()->GetTree()->Branch("Micro1_E_col1",&Micro1_E_col1,"Micro1_E_col1/D"); - RootOutput::getInstance()->GetTree()->Branch("Micro1_E_col2",&Micro1_E_col2,"Micro1_E_col2/D"); - RootOutput::getInstance()->GetTree()->Branch("Micro1_E_col3",&Micro1_E_col3,"Micro1_E_col3/D"); - RootOutput::getInstance()->GetTree()->Branch("Micro1_E_col5",&Micro1_E_col5,"Micro1_E_col5/D"); - RootOutput::getInstance()->GetTree()->Branch("Micro1_E_col6",&Micro1_E_col6,"Micro1_E_col6/D"); - RootOutput::getInstance()->GetTree()->Branch("Micro1_E_col7",&Micro1_E_col7,"Micro1_E_col7/D"); - RootOutput::getInstance()->GetTree()->Branch("Micro2_E_col1",&Micro2_E_col1,"Micro2_E_col1/D"); - RootOutput::getInstance()->GetTree()->Branch("Micro2_E_col2",&Micro2_E_col2,"Micro2_E_col2/D"); - RootOutput::getInstance()->GetTree()->Branch("Micro2_E_col3",&Micro2_E_col3,"Micro2_E_col3/D"); - RootOutput::getInstance()->GetTree()->Branch("Micro2_E_col4",&Micro2_E_col4,"Micro2_E_col4/D"); - RootOutput::getInstance()->GetTree()->Branch("Micro2_E_col5",&Micro2_E_col5,"Micro2_E_col5/D"); - RootOutput::getInstance()->GetTree()->Branch("Micro2_E_col6",&Micro2_E_col6,"Micro2_E_col6/D"); - RootOutput::getInstance()->GetTree()->Branch("Micro2_E_col7",&Micro2_E_col7,"Micro2_E_col7/D"); +//by Shuya 170912 + RootOutput::getInstance()->GetTree()->Branch("Micro1_E_col1_sum",&Micro1_E_col1_sum,"Micro1_E_col1_sum/D"); + RootOutput::getInstance()->GetTree()->Branch("Micro1_E_col2_sum",&Micro1_E_col2_sum,"Micro1_E_col2_sum/D"); + RootOutput::getInstance()->GetTree()->Branch("Micro1_E_col3_sum",&Micro1_E_col3_sum,"Micro1_E_col3_sum/D"); + RootOutput::getInstance()->GetTree()->Branch("Micro1_E_col4_sum",&Micro1_E_col4_sum,"Micro1_E_col4_sum/D"); + RootOutput::getInstance()->GetTree()->Branch("Micro1_E_col5_sum",&Micro1_E_col5_sum,"Micro1_E_col5_sum/D"); + RootOutput::getInstance()->GetTree()->Branch("Micro1_E_col6_sum",&Micro1_E_col6_sum,"Micro1_E_col6_sum/D"); + RootOutput::getInstance()->GetTree()->Branch("Micro1_E_col7_sum",&Micro1_E_col7_sum,"Micro1_E_col7_sum/D"); + RootOutput::getInstance()->GetTree()->Branch("Micro2_E_col1_sum",&Micro2_E_col1_sum,"Micro2_E_col1_sum/D"); + RootOutput::getInstance()->GetTree()->Branch("Micro2_E_col2_sum",&Micro2_E_col2_sum,"Micro2_E_col2_sum/D"); + RootOutput::getInstance()->GetTree()->Branch("Micro2_E_col3_sum",&Micro2_E_col3_sum,"Micro2_E_col3_sum/D"); + RootOutput::getInstance()->GetTree()->Branch("Micro2_E_col4_sum",&Micro2_E_col4_sum,"Micro2_E_col4_sum/D"); + RootOutput::getInstance()->GetTree()->Branch("Micro2_E_col5_sum",&Micro2_E_col5_sum,"Micro2_E_col5_sum/D"); + RootOutput::getInstance()->GetTree()->Branch("Micro2_E_col6_sum",&Micro2_E_col6_sum,"Micro2_E_col6_sum/D"); + RootOutput::getInstance()->GetTree()->Branch("Micro2_E_col7_sum",&Micro2_E_col7_sum,"Micro2_E_col7_sum/D"); + + //by Shuya 170912 + RootOutput::getInstance()->GetTree()->Branch("Micro1_E_col1_mult",&Micro1_E_col1_mult,"Micro1_E_col1_mult/D"); + RootOutput::getInstance()->GetTree()->Branch("Micro1_E_col2_mult",&Micro1_E_col2_mult,"Micro1_E_col2_mult/D"); + RootOutput::getInstance()->GetTree()->Branch("Micro1_E_col3_mult",&Micro1_E_col3_mult,"Micro1_E_col3_mult/D"); + RootOutput::getInstance()->GetTree()->Branch("Micro1_E_col4_mult",&Micro1_E_col4_mult,"Micro1_E_col4_mult/D"); + RootOutput::getInstance()->GetTree()->Branch("Micro1_E_col5_mult",&Micro1_E_col5_mult,"Micro1_E_col5_mult/D"); + RootOutput::getInstance()->GetTree()->Branch("Micro1_E_col6_mult",&Micro1_E_col6_mult,"Micro1_E_col6_mult/D"); + RootOutput::getInstance()->GetTree()->Branch("Micro1_E_col7_mult",&Micro1_E_col7_mult,"Micro1_E_col7_mult/D"); + RootOutput::getInstance()->GetTree()->Branch("Micro2_E_col1_mult",&Micro2_E_col1_mult,"Micro2_E_col1_mult/D"); + RootOutput::getInstance()->GetTree()->Branch("Micro2_E_col2_mult",&Micro2_E_col2_mult,"Micro2_E_col2_mult/D"); + RootOutput::getInstance()->GetTree()->Branch("Micro2_E_col3_mult",&Micro2_E_col3_mult,"Micro2_E_col3_mult/D"); + RootOutput::getInstance()->GetTree()->Branch("Micro2_E_col4_mult",&Micro2_E_col4_mult,"Micro2_E_col4_mult/D"); + RootOutput::getInstance()->GetTree()->Branch("Micro2_E_col5_mult",&Micro2_E_col5_mult,"Micro2_E_col5_mult/D"); + RootOutput::getInstance()->GetTree()->Branch("Micro2_E_col6_mult",&Micro2_E_col6_mult,"Micro2_E_col6_mult/D"); + RootOutput::getInstance()->GetTree()->Branch("Micro2_E_col7_mult",&Micro2_E_col7_mult,"Micro2_E_col7_mult/D"); //TACS RootOutput::getInstance()->GetTree()->Branch("TacSiGeOR",&TacSiGeOR,"TacSiGeOR/D"); RootOutput::getInstance()->GetTree()->Branch("TacSiMicro",&TacSiMicro,"TacSiMicro/D"); + + RootOutput::getInstance()->GetTree()->Branch("TacSiMicro_E",&TacSiMicro_E,"TacSiMicro_E/D"); RootOutput::getInstance()->GetTree()->Branch("TacSiMicro_dE",&TacSiMicro_dE,"TacSiMicro_dE/D"); - RootOutput::getInstance()->GetTree()->Branch("TacSiPlastLeft",& TacSiPlastLeft," TacSiPlastLeft/D"); - RootOutput::getInstance()->GetTree()->Branch("TacSiPlastRight",& TacSiPlastRight," TacSiPlastRight/D"); + + RootOutput::getInstance()->GetTree()->Branch("TacSiPlastLeft",&TacSiPlastLeft,"TacSiPlastLeft/D"); + RootOutput::getInstance()->GetTree()->Branch("TacSiPlastRight",&TacSiPlastRight,"TacSiPlastRight/D"); + // Other - RootOutput::getInstance()->GetTree()->Branch("RunNumber", &RunNumber," RunNumber/I"); + RootOutput::getInstance()->GetTree()->Branch("RunNumber", &RunNumber,"RunNumber/I"); // by Shuya 170524. - RootOutput::getInstance()->GetTree()->Branch("EntryNumber", &EntryNumber," EntryNumber/I"); + RootOutput::getInstance()->GetTree()->Branch("EntryNumber", &EntryNumber,"EntryNumber/I"); //Simulation //RootOutput::getInstance()->GetTree()->Branch("Original_ELab",&Original_ELab,"Original_ELab/D"); diff --git a/Projects/T40/Analysis.h b/Projects/T40/Analysis.h index ab8f31df7936426cb323cdb6276426818706a3b5..5d31cd9146bff1fe524cb17e2c0ffd814be3510e 100644 --- a/Projects/T40/Analysis.h +++ b/Projects/T40/Analysis.h @@ -53,6 +53,10 @@ class Analysis: public NPL::VAnalysis{ private: double Ex; +//by Shuya 170703 + double Ex_Hyball; + double Ex_Barrel; + double ELab; double ThetaLab; double ThetaCM; @@ -62,6 +66,8 @@ class Analysis: public NPL::VAnalysis{ TInitialConditions* Initial; NPL::Reaction* myReaction; bool LightParticleDetected; +//by Shuya 171019 + double PhiLab; // Energy loss table: the G4Table are generated by the simulation @@ -69,6 +75,8 @@ class Analysis: public NPL::VAnalysis{ NPL::EnergyLoss LightAl ; NPL::EnergyLoss LightSi; NPL::EnergyLoss BeamTarget ; +//by Shuya 170530 + NPL::EnergyLoss LightCBacking; double OriginalBeamEnergy; double FinalBeamEnergy; @@ -103,7 +111,8 @@ class Analysis: public NPL::VAnalysis{ static const Int_t kNumAw = 4; // number of wires double Delta_E ; double Micro1_E_row1 ; - double Micro1_E_col4 ; + //by Shuya 170912 + double Micro1_E_col4_sum ; double Micro1_E_row1_2; double Micro2_E ; double Micro2_E_row1_2; @@ -114,19 +123,35 @@ class Analysis: public NPL::VAnalysis{ double Aw_ThetaFit ; // Theta calculated from FITTING all wires double Aw_ThetaFit_R2; // Goodness of fit value of theta from fitting //by Shuya 170516 - double Micro1_E_col1 ; - double Micro1_E_col2 ; - double Micro1_E_col3 ; - double Micro1_E_col5 ; - double Micro1_E_col6 ; - double Micro1_E_col7 ; - double Micro2_E_col1 ; - double Micro2_E_col2 ; - double Micro2_E_col3 ; - double Micro2_E_col4 ; - double Micro2_E_col5 ; - double Micro2_E_col6 ; - double Micro2_E_col7 ; + //by Shuya 170912 + double Micro1_E_col1_sum ; + double Micro1_E_col2_sum ; + double Micro1_E_col3_sum ; + double Micro1_E_col5_sum ; + double Micro1_E_col6_sum ; + double Micro1_E_col7_sum ; + double Micro2_E_col1_sum ; + double Micro2_E_col2_sum ; + double Micro2_E_col3_sum ; + double Micro2_E_col4_sum ; + double Micro2_E_col5_sum ; + double Micro2_E_col6_sum ; + double Micro2_E_col7_sum ; + //by Shuya 170912 + double Micro1_E_col1_mult ; + double Micro1_E_col2_mult ; + double Micro1_E_col3_mult ; + double Micro1_E_col4_mult ; + double Micro1_E_col5_mult ; + double Micro1_E_col6_mult ; + double Micro1_E_col7_mult ; + double Micro2_E_col1_mult ; + double Micro2_E_col2_mult ; + double Micro2_E_col3_mult ; + double Micro2_E_col4_mult ; + double Micro2_E_col5_mult ; + double Micro2_E_col6_mult ; + double Micro2_E_col7_mult ; //TACS double TacSiGeOR ; diff --git a/Projects/T40/T40.detector b/Projects/T40/T40.detector index 96c5a2c1c284344019ea0c4701eee16d4f63ca16..f0a35d8eae29cf8edab72fa0f062b7d445c739dc 100644 --- a/Projects/T40/T40.detector +++ b/Projects/T40/T40.detector @@ -3,9 +3,9 @@ GeneralTarget %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %0.5mg/cm2 Target - THICKNESS= 0.0 micrometer + THICKNESS= 0.1132 micrometer RADIUS= 5 mm - MATERIAL= CD2 + MATERIAL= 6LiF ANGLE= 0 deg X= 0 mm Y= 0 mm diff --git a/Projects/T40/configs/ConfigGeTamu.dat b/Projects/T40/configs/ConfigGeTamu.dat index 2cdacb530b4094c1c27808dadb3325a589422039..52af7baae72f2da9bb2bc74debade14a49165141 100644 --- a/Projects/T40/configs/ConfigGeTamu.dat +++ b/Projects/T40/configs/ConfigGeTamu.dat @@ -6,10 +6,12 @@ ConfigGeTamu %DISABLE_ALL CLOVER04 %DISABLE_CHANNEL CLOVER01_SEG01 %%%%%%%% Thresh. in "channels" for Raw, in keV for Calibrated - CRY_E_RAW_THRESHOLD 300 - SEG_E_RAW_THRESHOLD 300 - CRY_E_THRESHOLD 100 - SEG_E_THRESHOLD 100 + %CRY_E_RAW_THRESHOLD 300 + %SEG_E_RAW_THRESHOLD 300 + CRY_E_RAW_THRESHOLD 50 + SEG_E_RAW_THRESHOLD 50 + CRY_E_THRESHOLD 50 + SEG_E_THRESHOLD 50 %%%%%%%% ADD Back scheme for the analysis ADD_BACK_CLOVER %ADD_BACK_ARRAY diff --git a/Projects/T40/configs/ConfigTiaraHyball.dat b/Projects/T40/configs/ConfigTiaraHyball.dat new file mode 100644 index 0000000000000000000000000000000000000000..2424c7e306f8b9fc4d2696e5829cf76a618d404a --- /dev/null +++ b/Projects/T40/configs/ConfigTiaraHyball.dat @@ -0,0 +1,21 @@ +ConfigTiaraHyball + %MAX_STRIP_MULTIPLICITY + %Energy match up is made by if abs((RingE-SectorE)/2) < STRIP_ENERGY_MATCHING_SIGMA*STRIP_ENERGY_MATCHING_NUMBER_OF_SIGMA + %in unit of MeV + STRIP_ENERGY_MATCHING_SIGMA 0.01 + STRIP_ENERGY_MATCHING_NUMBER_OF_SIGMA 5 + %If you want to use Ring signal instead of Ring-Sector average, comment this out + %To get weighted averaged energy, you need to input these below (default; RingE 18 keV, SectorE=22 keV) + TAKE_E_RING_SECTOR_Average + SIGMA_RING_E 0.018 + SIGMA_SECTOR_E 0.022 + TAKE_E_RING + %TAKE_E_SECTOR + %TAKE_T_RING + %TAKE_T_SECTOR + %STRIP_RING_E_RAW_THRESHOLD + %STRIP_SECTOR_E_RAW_THRESHOLD + %STRIP_RING_E_THRESHOLD + %STRIP_SECTOR_THRESHOLD + %DISABLE_ALL + %DISABLE_CHANNEL