Skip to content
Snippets Groups Projects
Commit aae16765 authored by Pierre Morfouace's avatar Pierre Morfouace
Browse files

* Updating s455 project

parent 8bff9f17
No related branches found
No related tags found
No related merge requests found
Pipeline #144004 passed
Stage: build-NPLib
Stage: build-NPSimulation
Stage: test
Hide whitespace changes
Inline Side-by-side
Projects/s455/Analysis.cxx
+ 280
− 223
View file @ aae16765
......@@ -56,6 +56,7 @@ void Analysis::Init(){
void Analysis::TreatEvent(){
ReInitValue();
//cout << "************" << endl;
RunID = fRunID;
BeamAnalysis();
unsigned int sofsci_size = SofSci->DetectorNbr.size();
......@@ -67,8 +68,6 @@ void Analysis::TreatEvent(){
FissionFragmentAnalysis();
}
}
////////////////////////////////////////////////////////////////////////////////
......@@ -127,7 +126,6 @@ void Analysis::FissionFragmentAnalysis(){
int multR = mult3 + mult4;
if(softwim_size>1){
//if( (multL>15 && multR>15) || (multR==32) || (multL==32)){
if( (mult1>1 && mult1<17) || (mult2>1 && mult2<17) || (mult3>1 && mult3<17) || (mult4>1 && mult4<17)){
for(unsigned int i=0; i< softwim_size; i++){
int sec = SofTwim->SectionNbr[i];
......@@ -247,9 +245,6 @@ void Analysis::FissionFragmentAnalysis(){
// Z calibration //
double p0 = 2.80063;
double p1 = 6.91985e-2;
double p2 = 1.01598e-7;
double Z1=-1;
double Z2=-1;
double Zsum=-1;
......@@ -281,9 +276,9 @@ void Analysis::FissionFragmentAnalysis(){
}
if(Z1>0 && Z2>0){
Z1 = p0 + p1*Z1 + p2*Z1*Z1;
Z2 = p0 + p1*Z2 + p2*Z2*Z2;
Z1 = fZff_p0 + fZff_p1*Z1 + fZff_p2*Z1*Z1;
Z2 = fZff_p0 + fZff_p1*Z2 + fZff_p2*Z2*Z2;
Z1 = sqrt(Z1);
Z2 = sqrt(Z2);
......@@ -301,15 +296,6 @@ void Analysis::FissionFragmentAnalysis(){
SofFF->SetBeta(Beta3);
SofFF->SetBeta(Beta4);
/*SofFF->SetZ(E1);
SofFF->SetZ(E2);
SofFF->SetZ(E3);
SofFF->SetZ(E4);*/
/*SofFF->SetZ( sqrt(p0+p1*E1+p2*E1*E1) );
SofFF->SetZ( sqrt(p0+p1*E2+p2*E2*E2) );
SofFF->SetZ( sqrt(p0+p1*E3+p2*E3*E3) );
SofFF->SetZ( sqrt(p0+p1*E4+p2*E4*E4) );*/
SofFF->SetZ(Z1);
SofFF->SetZ(Z2);
......@@ -322,234 +308,305 @@ void Analysis::FissionFragmentAnalysis(){
SofFF->SetIntZsum(iZsum);
}
}
}
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::BeamAnalysis(){
unsigned int sofsci_size = SofSci->DetectorNbr.size();
double Y_p0 = 23943.8;
double Y_p1 = 12.362;
double alpha = 2.09862;
double beta = 0.972689;
double Z_p0 = alpha+beta*-9.31873;
double Z_p1 = beta*0.564459;
if(sofsci_size==2){
double beta = SofSci->Beta[0];
//cout << "Set beta to " << beta << endl;
SofTrim->SetBeta(beta);
SofTrim->BuildSimplePhysicalEvent();
double Zbeam,Qmax,Theta;
if(SofTrim->EnergySection.size()>0){
Zbeam = SofTrim->GetMaxEnergySection();
Qmax = DetermineQmax();
Theta = SofTrim->Theta[0];
double TofFromS2 = SofSci->CalTof[0];
double velocity_mns = SofSci->VelocityMNs[0];
double Beta = SofSci->Beta[0];
double XS2 = SofSci->PosMm[0];
//double XCC = SofSci->PosMm[1];
double XCC=0;
double YCC=0;
for(unsigned int i=0; i<SofMwpc->DetectorNbr.size(); i++){
if(SofMwpc->DetectorNbr[i]==1){
XCC = SofMwpc->PositionX1[i];
YCC = SofMwpc->PositionY[i];
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::BeamAnalysis(){
unsigned int sofsci_size = SofSci->DetectorNbr.size();
if(sofsci_size==2){
double beta = SofSci->Beta[0];
//cout << "Set beta to " << beta << endl;
SofTrim->SetBeta(beta);
SofTrim->BuildSimplePhysicalEvent();
double Zbeam,Qmax,Theta;
if(SofTrim->EnergySection.size()>0){
Zbeam = SofTrim->GetMaxEnergySection();
Qmax = DetermineQmax();
Theta = SofTrim->Theta[0];
double TofFromS2 = SofSci->CalTof[0];
double velocity_mns = SofSci->VelocityMNs[0];
double Beta = SofSci->Beta[0];
double XS2 = SofSci->PosMm[0];
//double XCC = SofSci->PosMm[1];
double XCC=0;
double YCC=0;
for(unsigned int i=0; i<SofMwpc->DetectorNbr.size(); i++){
if(SofMwpc->DetectorNbr[i]==1){
XCC = SofMwpc->PositionX1[i];
YCC = SofMwpc->PositionY[i];
}
double LS2;
LS2 = fLS2_0;//*(1 + fK_LS2*Theta);
velocity_mns = LS2/TofFromS2;
Beta = velocity_mns * m/ns / NPUNITS::c_light;
double Gamma = 1./(TMath::Sqrt(1 - TMath::Power(Beta,2)));
double Brho = fBrho0 * (1 - XS2/fDS2 - XCC/fDCC);
double AoQ = Brho / (3.10716*Gamma*Beta);
double A = AoQ * Qmax;
// Y dependence correction //
Zbeam = Zbeam/(Y_p0 + Y_p1*YCC)*Y_p0;
// Z calibration //
Zbeam = Z_p0 + Z_p1*sqrt(Zbeam);
//Zbeam = 2.09862 + 0.972689*Zbeam;
// Filling Beam tree
SofBeamID->SetZbeam(Zbeam);
SofBeamID->SetQmax(rand.Gaus(Qmax,0.15));
SofBeamID->SetAoQ(AoQ);
SofBeamID->SetAbeam(A);
SofBeamID->SetBeta(Beta);
SofBeamID->SetGamma(Gamma);
SofBeamID->SetBrho(Brho);
SofBeamID->SetXS2(XS2);
SofBeamID->SetXCC(XCC);
SofBeamID->SetYCC(YCC);
}
}
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::LoadCut(){
TString input_path = "./calibration/SofTrim/cut/";
TString rootfile;
TString cutfile;
TFile* file;
for(int i=0; i<3; i++){
// Q=77
rootfile = Form("cutsec%iQ77.root",i+1);
cutfile = input_path + rootfile;
file = new TFile(cutfile);
cutQ77[i] = (TCutG*) file->Get(Form("cutsec%iQ77",i+1));
// Q=78
rootfile = Form("cutsec%iQ78.root",i+1);
cutfile = input_path + rootfile;
file = new TFile(cutfile);
cutQ78[i] = (TCutG*) file->Get(Form("cutsec%iQ78",i+1));
// Q=79
rootfile = Form("cutsec%iQ79.root",i+1);
cutfile = input_path + rootfile;
file = new TFile(cutfile);
cutQ79[i] = (TCutG*) file->Get(Form("cutsec%iQ79",i+1));
// Q=80
rootfile = Form("cutsec%iQ80.root",i+1);
cutfile = input_path + rootfile;
file = new TFile(cutfile);
cutQ80[i] = (TCutG*) file->Get(Form("cutsec%iQ80",i+1));
double LS2;
LS2 = fLS2_0;//*(1 + fK_LS2*Theta);
velocity_mns = LS2/TofFromS2;
Beta = velocity_mns * m/ns / NPUNITS::c_light;
double Gamma = 1./(TMath::Sqrt(1 - TMath::Power(Beta,2)));
double Brho = fBrho0 * (1 - XS2/fDS2 - XCC/fDCC);
double AoQ = Brho / (3.10716*Gamma*Beta);
double A = AoQ * Qmax;
// Y dependence correction //
double Y_p0 = 23943.8;
double Y_p1 = 12.362;
Zbeam = Zbeam/(Y_p0 + Y_p1*YCC)*Y_p0;
// Z calibration //
Zbeam = fZbeam_p0 + fZbeam_p1*Zbeam + fZbeam_p2*Zbeam*Zbeam;
Zbeam = sqrt(Zbeam);
// Last beta correction //
double Beta_norm = 0.8355;
Zbeam = Zbeam/(fZBeta_p0 + fZBeta_p1*Beta)*(fZBeta_p0 + fZBeta_p1*Beta_norm);
// Filling Beam tree
SofBeamID->SetZbeam(Zbeam);
SofBeamID->SetQmax(rand.Gaus(Qmax,0.15));
SofBeamID->SetAoQ(AoQ);
SofBeamID->SetAbeam(A);
SofBeamID->SetBeta(Beta);
SofBeamID->SetGamma(Gamma);
SofBeamID->SetBrho(Brho);
SofBeamID->SetXS2(XS2);
SofBeamID->SetXCC(XCC);
SofBeamID->SetYCC(YCC);
}
}
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::LoadCut(){
TString input_path = "./calibration/SofTrim/cut/";
TString rootfile;
TString cutfile;
TFile* file;
for(int i=0; i<3; i++){
// Q=77
rootfile = Form("cutsec%iQ77.root",i+1);
cutfile = input_path + rootfile;
file = new TFile(cutfile);
cutQ77[i] = (TCutG*) file->Get(Form("cutsec%iQ77",i+1));
// Q=78
rootfile = Form("cutsec%iQ78.root",i+1);
cutfile = input_path + rootfile;
file = new TFile(cutfile);
cutQ78[i] = (TCutG*) file->Get(Form("cutsec%iQ78",i+1));
// Q=79
rootfile = Form("cutsec%iQ79.root",i+1);
cutfile = input_path + rootfile;
file = new TFile(cutfile);
cutQ79[i] = (TCutG*) file->Get(Form("cutsec%iQ79",i+1));
// Q=80
rootfile = Form("cutsec%iQ80.root",i+1);
cutfile = input_path + rootfile;
file = new TFile(cutfile);
cutQ80[i] = (TCutG*) file->Get(Form("cutsec%iQ80",i+1));
}
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::LoadSpline(){
TString input_path = "./calibration/SofTwim/spline/";
TString rootfile = input_path + "spline_beta.root";
TFile* ifile = new TFile(rootfile,"read");
TString splinename;
if(ifile->IsOpen()){
cout << "Loading Beta spline for fission fragment analysis..." << endl;
for(int i=0; i<4; i++){
splinename = Form("spline_beta_sec%i",i+1);
fcorr_z_beta[i] = (TSpline3*) ifile->FindObjectAny(splinename);
}
ifile->Close();
////////////////////////////////////////////////////////////////////////////////
void Analysis::LoadSpline(){
TString input_path = "./calibration/SofTwim/spline/";
TString rootfile = input_path + "spline_beta.root";
TFile* ifile = new TFile(rootfile,"read");
TString splinename;
if(ifile->IsOpen()){
cout << "Loading Beta spline for fission fragment analysis..." << endl;
for(int i=0; i<4; i++){
splinename = Form("spline_beta_sec%i",i+1);
fcorr_z_beta[i] = (TSpline3*) ifile->FindObjectAny(splinename);
}
else
cout << "File " << rootfile << " not found!" << endl;
//*** ***//
rootfile = input_path + "spline_dt.root";
ifile = new TFile(rootfile,"read");
if(ifile->IsOpen()){
cout << "Loading DT spline for fission fragment analysis..." << endl;
for(int i=0; i<4; i++){
splinename = Form("spline_dt_sec%i",i+1);
fcorr_z_dt[i] = (TSpline3*) ifile->FindObjectAny(splinename);
}
ifile->Close();
ifile->Close();
}
else
cout << "File " << rootfile << " not found!" << endl;
//*** ***//
rootfile = input_path + "spline_dt.root";
ifile = new TFile(rootfile,"read");
if(ifile->IsOpen()){
cout << "Loading DT spline for fission fragment analysis..." << endl;
for(int i=0; i<4; i++){
splinename = Form("spline_dt_sec%i",i+1);
fcorr_z_dt[i] = (TSpline3*) ifile->FindObjectAny(splinename);
}
else
cout << "File " << rootfile << " not found!" << endl;
ifile->Close();
}
else
cout << "File " << rootfile << " not found!" << endl;
}
////////////////////////////////////////////////////////////////////////////////
int Analysis::DetermineQmax(){
int Qmax;
int Qsec[3];
unsigned int size = SofTrim->EnergySection.size();
for(unsigned int i=0; i<size; i++){
int SecNbr = SofTrim->SectionNbr[i];
double Theta = SofTrim->Theta[i];
double Esec = SofTrim->EnergySection[i];
if(cutQ77[SecNbr-1]->IsInside(Theta,Esec))
Qsec[SecNbr-1] = 77;
else if(cutQ78[SecNbr-1]->IsInside(Theta,Esec))
Qsec[SecNbr-1] = 78;
else if(cutQ79[SecNbr-1]->IsInside(Theta,Esec))
Qsec[SecNbr-1] = 79;
else if(cutQ80[SecNbr-1]->IsInside(Theta,Esec))
Qsec[SecNbr-1] = 80;
//else if(cutQ81[SecNbr-1]->IsInside(Theta,Esec))
//Qsec[SecNbr-1] = 81;
}
////////////////////////////////////////////////////////////////////////////////
int Analysis::DetermineQmax(){
int Qmax;
int Qsec[3];
unsigned int size = SofTrim->EnergySection.size();
for(unsigned int i=0; i<size; i++){
int SecNbr = SofTrim->SectionNbr[i];
double Theta = SofTrim->Theta[i];
double Esec = SofTrim->EnergySection[i];
if(cutQ77[SecNbr-1]->IsInside(Theta,Esec))
Qsec[SecNbr-1] = 77;
else if(cutQ78[SecNbr-1]->IsInside(Theta,Esec))
Qsec[SecNbr-1] = 78;
else if(cutQ79[SecNbr-1]->IsInside(Theta,Esec))
Qsec[SecNbr-1] = 79;
else if(cutQ80[SecNbr-1]->IsInside(Theta,Esec))
Qsec[SecNbr-1] = 80;
//else if(cutQ81[SecNbr-1]->IsInside(Theta,Esec))
//Qsec[SecNbr-1] = 81;
}
Qmax = max(Qsec[0],Qsec[1]);
Qmax = max(Qsec[2],Qmax);
Qmax = max(Qsec[0],Qsec[1]);
Qmax = max(Qsec[2],Qmax);
return Qmax;
}
return Qmax;
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::End(){
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::End(){
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitParameter(){
fLS2_0 = 135.614;
fDS2 = 8000;
fDCC = -10000;
fK_LS2 = -30e-8;
fRunID = 8;
// Beam parameter //
fZBeta_p0 = 1;
fZBeta_p1 = 0;
// FF paramter //
fZff_p0 = 2.80063;
fZff_p1 = 6.91985e-2;
fZff_p2 = 1.01598e-7;
if(fRunID==1 || fRunID==2){
//fBrho0 = 10.6813; // 180Hg
fBrho0 = 10.6955; // 180Hg
fZbeam_p0 = -5303.06;
fZbeam_p1 = 0.674945;
fZbeam_p2 = -8.32085e-6;
fZBeta_p0 = 2.73;
fZBeta_p1 = 74.7153;
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitParameter(){
fLS2_0 = 135.614;
fDS2 = 8000;
fDCC = -10000;
fK_LS2 = -30e-8;
//fBrho0 = 14.1008; // 238U run 367
//fBrho0 = 12.8719; // 238U run 368
//fBrho0 = 12.3255; // 238U run 369
//fBrho0 = 12.3352; // 238U run 422
//fBrho0 = 10.9476; // 189Pb
//fBrho0 = 10.9558; // 184Hg
//fBrho0 = 10.8183; // 182Hg
//fBrho0 = 10.6814; // 180Hg
//fBrho0 = 10.8138; // 187Pb
//fBrho0 = 11.3418; // 216Th
if(fRunID==3){
fBrho0 = 10.8183; // 182Hg
fZbeam_p0 = -2737.25;
fZbeam_p1 = 0.452017;
fZbeam_p2 = -3.48831e-6;
}
if(fRunID==4){
fBrho0 = 10.9558; // 184Hg
fZbeam_p0 = -5044.61;
fZbeam_p1 = 0.639986;
fZbeam_p2 = -7.3077e-6;
}
if(fRunID==5){
fBrho0 = 10.8138; // 187Pb
fZbeam_p0 = -2858.72;
fZbeam_p1 = 0.454064;
fZbeam_p2 = -3.36443e-6;
}
if(fRunID==6){
fBrho0 = 10.9476; // 189Pb
fZbeam_p0 = 1590.66;
fZbeam_p1 = 0.0956455;
fZbeam_p2 = 3.84585e-6;
}
if(fRunID==7){
fBrho0 = 10.6814; // 175Pt
fZbeam_p0 = 459.68;
fZbeam_p1 = 0.162277;
fZbeam_p2 = 3.10164e-6;
}
if(fRunID==8){
fBrho0 = 11.0864; // 204Fr
//fBrho0 = 11.2712; // 207Fr
//fBrho0 = 10.6814; // 175Pt
//fBrho0 = 11.5067; // 221Pa
//fBrho0 = 11.0955; // 199At run 423 & 424
//fBrho0 = 10.9970; // 199At run 425 & 426
//fBrho0 = 10.8697; //197At
fZbeam_p0 = 4122.94;
fZbeam_p1 = -0.119867;
fZbeam_p2 = 8.29115e-6;
fZBeta_p0 = 63.9575;
fZBeta_p1 = 25.1988;
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::ReInitValue(){
SofBeamID->Clear();
SofFF->Clear();
if(fRunID==9){
fBrho0 = 11.2712; // 207Fr
fZbeam_p0 = -1752.27;
fZbeam_p1 = 0.346018;
fZbeam_p2 = -8.64673e-7;
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitOutputBranch(){
//RootOutput::getInstance()->GetTree()->Branch("Zbeam",&Zbeam,"Zbeam/D");
RootOutput::getInstance()->GetTree()->Branch("SofBeamID","TSofBeamID",&SofBeamID);
RootOutput::getInstance()->GetTree()->Branch("SofFissionFragment","TSofFissionFragment",&SofFF);
if(fRunID==10){
fBrho0 = 11.0955; // 199At run 423 & 424
fZbeam_p0 = -116.425;
fZbeam_p1 = 0.218256;
fZbeam_p2 = 1.62399e-6;
}
if(fRunID==11){
fBrho0 = 10.9970; // 199At run 425 & 426
fZbeam_p0 = -116.425;
fZbeam_p1 = 0.218256;
fZbeam_p2 = 1.62399e-6;
}
////////////////////////////////////////////////////////////////////////////////
// Construct Method to be pass to the DetectorFactory //
////////////////////////////////////////////////////////////////////////////////
NPL::VAnalysis* Analysis::Construct(){
return (NPL::VAnalysis*) new Analysis();
if(fRunID==12){
fBrho0 = 10.8697; //197At
fZbeam_p0 = -2683.52;
fZbeam_p1 = 0.422551;
fZbeam_p2 = -2.44471e-6;
}
if(fRunID==13){
fBrho0 = 11.3418; // 216Th
}
if(fRunID==14){
fBrho0 = 11.5067; // 221Pa
}
}
////////////////////////////////////////////////////////////////////////////////
// Registering the construct method to the factory //
////////////////////////////////////////////////////////////////////////////////
extern "C"{
class proxy{
public:
proxy(){
NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
}
};
////////////////////////////////////////////////////////////////////////////////
void Analysis::ReInitValue(){
SofBeamID->Clear();
SofFF->Clear();
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitOutputBranch(){
RootOutput::getInstance()->GetTree()->Branch("RunID",&RunID,"RunID/I");
RootOutput::getInstance()->GetTree()->Branch("SofBeamID","TSofBeamID",&SofBeamID);
RootOutput::getInstance()->GetTree()->Branch("SofFissionFragment","TSofFissionFragment",&SofFF);
}
////////////////////////////////////////////////////////////////////////////////
// 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;
}
proxy p;
}
Projects/s455/Analysis.h
+ 10
− 0
View file @ aae16765
......@@ -67,13 +67,23 @@ class Analysis: public NPL::VAnalysis{
TSofTofWPhysics* SofTofW;
TSofBeamID* SofBeamID;
TSofFissionFragment* SofFF;
int RunID;
private:
int fRunID;
double fLS2_0;
double fBrho0;
double fDS2;
double fDCC;
double fK_LS2;
double fZbeam_p0;
double fZbeam_p1;
double fZbeam_p2;
double fZff_p0;
double fZff_p1;
double fZff_p2;
double fZBeta_p0;
double fZBeta_p1;
TRandom3 rand;
TCutG* cutQ77[3];
......
Projects/s455/calibration/SofSci/SofSci_physics.cal
+ 2
− 2
View file @ aae16765
SofSci_TOF2INV_V -8.04267 0.00737389
SofSci_TOF2INV_V -8.0425 0.00737389
SofSci_LENGTH_S2 135.614
SofSci_DET1_POSPAR 90 54.6
SofSci_DET1_POSPAR 45.0 80.0
SofSci_DET2_POSPAR 943.287 86.652
Projects/s455/calibration/SofTrim/SofTrim_SectionAlign.cal
+ 2
− 2
View file @ aae16765
SofTrim_SEC1_ALIGN 755.97 0.987851
SofTrim_SEC1_ALIGN -348.015 1.03199
SofTrim_SEC2_ALIGN 0 1
SofTrim_SEC3_ALIGN -1674.22 1.15275
SofTrim_SEC3_ALIGN 2496.29 0.984162
Projects/s455/configs/ConfigSofSci.dat
+ 2
− 2
View file @ aae16765
......@@ -3,6 +3,6 @@ ConfigSofSci
DET1_POSNS_MAX 20
DET2_POSNS_MIN -20
DET2_POSNS_MAX 20
RAWTOF_MIN 1600
RAWTOF_MAX 1700
RAWTOF_MIN 1628
RAWTOF_MAX 1636
Projects/s455/macro/GetChargeDist.C 0 → 100644
+ 352
− 0
View file @ aae16765
TChain* chain=NULL;
TCutG* cut_Pb1[13];
TCutG* cut_Pb2[13];
TCutG* cut_C[13];
int A;
int Z;
NPL::Particle* iso;
TH2F* hbid;
TH1F* hsumpb;
TH1F* hsumc;
TH1F* hzpb;
TH1F* hzc;
TH1F* hzcscaled;
TH1F* hsumcscaled;
TH1F* hzsub;
void LoadRootFiles(string iso);
//////////////////////////////////////
void LoadCuts(){
TString element[13] = {"180Hg_1", "180Hg_2", "182Hg", "184Hg", "187Pb", "189Pb", "175Pt", "204Fr", "207Fr", "199At_1", "199At_2", "197At", "216Th"};
TFile* file;
TString filename1[13];
TString filename2[13];
TString filename3[13];
for(int i=0; i<13; i++){
filename1[i]= "cuts/"+element[i]+"/cut_Pb1.root";
filename2[i]= "cuts/"+element[i]+"/cut_Pb2.root";
filename3[i]= "cuts/"+element[i]+"/cut_C.root";
file = new TFile(filename1[i],"read");
cut_Pb1[i] = (TCutG*) file->FindObjectAny("cut_Pb1");
cut_Pb1[i]->SetName(element[i]+"_Pb1");
file = new TFile(filename2[i],"read");
cut_Pb2[i] = (TCutG*) file->FindObjectAny("cut_Pb2");
cut_Pb2[i]->SetName(element[i]+"_Pb2");
file = new TFile(filename3[i],"read");
cut_C[i] = (TCutG*) file->FindObjectAny("cut_C");
cut_C[i]->SetName(element[i]+"_C");
}
}
//////////////////////////////////////
void GetChargeDist(string nucleus_name){
iso = new NPL::Particle(nucleus_name);
A = iso->GetA();
Z = iso->GetZ();
double AoQ=(double)A/(double)Z;
cout << "*******************************" << endl;
cout << "Nucleus name = " << nucleus_name << endl;
cout << "A= " << A << endl;
cout << "Z= " << Z << endl;
cout << "AoQ= " << AoQ << endl;
LoadRootFiles(nucleus_name);
LoadCuts();
//*** opening output root file ***//
TString output_filename = "root/charge_dist_" + nucleus_name + ".root";
TFile* ofile = new TFile(output_filename, "recreate");
cout << "Creating file: " << output_filename << endl;
ofile->cd();
hbid = new TH2F("hbid","hbid",1000,2.2,2.5,1000,70,95);
hsumpb = new TH1F("hsumpb","hsumpb",1000,50,95);
hsumc = new TH1F("hsumc","hsumc",1000,50,95);
hzpb = new TH1F("hzpb","hzpb",1000,20,60);
hzc = new TH1F("hzc","hzc",1000,20,60);
hzcscaled = new TH1F("hzcscaled","hzcscaled",1000,20,60);
hsumcscaled = new TH1F("hsumcscaled","hsumcscaled",1000,50,95);
hzsub = new TH1F("hzsub","hzsub",1000,20,60);
//*** Defining beam cut ***//
TCutG* cut_beam = new TCutG("cut_beam");
cut_beam->SetName("cut_beam");
cut_beam->SetVarX("fBeam_AoQ");
cut_beam->SetVarY("fBeam_Z");
double Rx=0.005;
double Ry=0.5;
AoQ = AoQ;// + 0.001;
double Zcenter = Z;// - 0.5;
for(int i=0; i<360; i++){
double x = AoQ+Rx*cos(i*TMath::Pi()/180);
double y = Zcenter+Ry*sin(i*TMath::Pi()/180);
cut_beam->SetPoint(i,x,y);
}
cut_beam->SetLineWidth(2);
cut_beam->SetLineColor(2);
int nentries = chain->GetEntries();
//*** Init input branch ***//
int RunID;
TSofBeamID* SofBeam = new TSofBeamID();
TSofAtPhysics* SofAt = new TSofAtPhysics();
TSofFissionFragment* SofFF = new TSofFissionFragment();
chain->SetBranchAddress("RunID",&RunID);
chain->SetBranchAddress("SofBeamID",&SofBeam);
chain->SetBranchAddress("SofAt",&SofAt);
chain->SetBranchAddress("SofFissionFragment",&SofFF);
for(int i=0; i<nentries; i++){
chain->GetEntry(i);
if(i%100000==0){
cout << "\033[34m\r Processing tree... " << (double)i/nentries*100 << "\% done" << flush;
}
if(SofAt->Energy.size()==3 || SofAt->Energy.size()==4){
double Zbeam = SofBeam->GetZbeam();
double AoQ = SofBeam->GetAoQ();
double Anode1 = SofAt->Energy[0];
double Anode2 = SofAt->Energy[1];
double Anode3 = SofAt->Energy[2];
int k=RunID-1;
if(cut_Pb1[k]->IsInside(Anode1,Anode2) || cut_Pb2[k]->IsInside(Anode2,Anode3)){
hbid->Fill(AoQ,Zbeam);
}
if(cut_beam->IsInside(AoQ,Zbeam)){
//*** Fission Fragment information ***//
double Zsum = SofFF->GetZsum();
int iZsum = SofFF->GetIntZsum();
double Z1 = -1;
double Z2 = -1;
if(SofFF->GetMult()>0){
Z1 = SofFF->GetZ(0);
Z2 = SofFF->GetZ(1);
}
//*** Taking event from fission in Lead cathod ***//
if(cut_Pb1[k]->IsInside(Anode1,Anode2) || cut_Pb2[k]->IsInside(Anode2,Anode3)){
hsumpb->Fill(Zsum);
if(iZsum==Z){
hzpb->Fill(Z1);
hzpb->Fill(Z2);
hzsub->Fill(Z1);
hzsub->Fill(Z2);
}
}
//*** Taking event from fission in Carbon cathod ***//
if(cut_C[k]->IsInside(Anode2,Anode3)){
hsumc->Fill(Zsum);
hsumcscaled->Fill(Zsum);
if(iZsum==Z){
hzc->Fill(Z1);
hzc->Fill(Z2);
hzcscaled->Fill(Z1);
hzcscaled->Fill(Z2);
}
}
}
}
} // end of loop ovent nentries
double integral_pb = hsumpb->Integral(hsumpb->GetXaxis()->FindBin(65.5), hsumpb->GetXaxis()->FindBin(Z-1.5));
double integral_c = hsumc->Integral(hsumc->GetXaxis()->FindBin(65.5), hsumc->GetXaxis()->FindBin(Z-1.5));
double ratio = integral_pb/integral_c;
if(ratio>0){
cout << "*** ratio= " << ratio << endl;
hsumcscaled->Scale(ratio);
hzcscaled->Scale(ratio);
hzsub->Add(hzcscaled,-1);
}
cout << endl;
hbid->Write();
hsumpb->Write();
hsumc->Write();
hzpb->Write();
hzc->Write();
hzcscaled->Write();
hsumcscaled->Write();
hzsub->Write();
cut_beam->Write();
}
//////////////////////////////////////
void LoadRootFiles(string iso){
chain = new TChain("PhysicsTree");
//*** Th isotopes ***//
if(iso=="217Th" || iso=="216Th" || iso=="209Ra" || iso=="210Ra" || iso=="211Ra" || iso=="212Ra" || iso=="213Ra" || iso=="214Ra"){
chain->Add("../../../Outputs/Analysis/216Th_1.root");
}
//*** Fr isotopes ***//
if(iso=="203Fr" || iso=="204Fr" || iso=="205Fr" || iso=="206Fr"){
chain->Add("../../../Outputs/Analysis/204Fr.root");
}
if(iso=="205Fr" || iso=="206Fr" || iso=="207Fr" || iso=="208Fr"){
chain->Add("../../../Outputs/Analysis/207Fr_2.root");
}
if(iso=="208Fr" || iso=="209Fr" || iso=="210Fr" || iso=="211Fr"){
chain->Add("../../../Outputs/Analysis/216Th_1.root");
}
//*** Rn isotopes ***//
if(iso=="200Rn" || iso=="201Rn" || iso=="202Rn" || iso=="203Rn" || iso=="204Rn" || iso=="205Rn"){
chain->Add("../../../Outputs/Analysis/204Fr.root");
}
if(iso=="203Rn" || iso=="204Rn" || iso=="205Rn" || iso=="206Rn"){
chain->Add("../../../Outputs/Analysis/207Fr.root");
}
//*** At isotopes ***//
if(iso=="197At" || iso=="198At" || iso=="199At" || iso=="200At" || iso=="201At" || iso=="202At"){
chain->Add("../../../Outputs/Analysis/204Fr.root");
}
if(iso=="200At" || iso=="201At" || iso=="202At" || iso=="203At"){
chain->Add("../../../Outputs/Analysis/207Fr.root");
}
if(iso=="198At" || iso=="199At" || iso=="200At" || iso=="201At"){
chain->Add("../../../Outputs/Analysis/199At_1.root");
chain->Add("../../../Outputs/Analysis/199At_2.root");
}
if(iso=="196At" || iso=="197At" || iso=="198At" || iso=="199At"){
chain->Add("../../../Outputs/Analysis/197At.root");
}
//*** Po isotopes ***//
if(iso=="192Po" || iso=="193Po" || iso=="194Po" || iso=="195Po" || iso=="196Po" || iso=="197Po"){
chain->Add("../../../Outputs/Analysis/197At.root");
}
if(iso=="195Po" || iso=="196Po" || iso=="197Po" || iso=="198Po" || iso=="199Po" || iso=="200Po"){
chain->Add("../../../Outputs/Analysis/199At_1.root");
chain->Add("../../../Outputs/Analysis/199At_2.root");
}
//*** Bi isotopes ***//
if(iso=="189Bi" || iso=="190Bi" || iso=="191Bi" || iso=="192Bi" || iso=="193Bi" || iso=="194Bi" || iso=="195Bi"){
chain->Add("../../../Outputs/Analysis/197At.root");
}
if(iso=="192Bi" || iso=="193Bi" || iso=="194Bi" || iso=="195Bi" || iso=="196Bi" || iso=="197Bi"){
chain->Add("../../../Outputs/Analysis/199At_1.root");
chain->Add("../../../Outputs/Analysis/199At_2.root");
}
//*** Pb isotopes ***//
if(iso=="186Pb" || iso=="187Pb" || iso=="188Pb" || iso=="189Pb" || iso=="190Pb" || iso=="191Pb" || iso=="192Pb"){
chain->Add("../../../Outputs/Analysis/197At.root");
}
if(iso=="191Pb" || iso=="192Pb" || iso=="193Pb" || iso=="194Pb"){
chain->Add("../../../Outputs/Analysis/199At_1.root");
chain->Add("../../../Outputs/Analysis/199At_2.root");
}
if(iso=="188Pb" || iso=="189Pb" || iso=="190Pb" || iso=="191Pb"){
chain->Add("../../../Outputs/Analysis/189Pb.root");
}
if(iso=="186Pb" || iso=="187Pb" || iso=="188Pb" || iso=="189Pb"){
chain->Add("../../../Outputs/Analysis/187Pb.root");
}
//*** Tl isotopes ***//
if(iso=="186Tl" || iso=="187Tl" || iso=="188Tl" || iso=="189Tl" || iso=="190Tl"){
chain->Add("../../../Outputs/Analysis/197At.root");
}
if(iso=="185Tl" || iso=="186Tl" || iso=="187Tl" || iso=="188Tl" || iso=="189Tl"){
chain->Add("../../../Outputs/Analysis/189Pb.root");
}
if(iso=="183Tl" || iso=="184Tl" || iso=="185Tl" || iso=="186Tl" || iso=="187Tl"){
chain->Add("../../../Outputs/Analysis/187Pb.root");
}
if(iso=="185Tl" || iso=="186Tl" || iso=="187Tl"){
chain->Add("../../../Outputs/Analysis/184Hg.root");
}
if(iso=="183Tl" || iso=="184Tl" || iso=="185Tl"){
chain->Add("../../../Outputs/Analysis/182Hg.root");
}
//*** Hg isotopes ***//
if(iso=="182Hg" || iso=="183Hg" || iso=="184Hg" || iso=="185Hg" || iso=="186Hg" || iso=="187Hg"){
chain->Add("../../../Outputs/Analysis/189Pb.root");
}
if(iso=="180Hg" || iso=="181Hg" || iso=="182Hg" || iso=="183Hg" || iso=="184Hg" || iso=="185Hg"){
chain->Add("../../../Outputs/Analysis/187Pb.root");
}
if(iso=="183Hg" || iso=="184Hg" || iso=="185Hg" || iso=="186Hg"){
chain->Add("../../../Outputs/Analysis/184Hg.root");
}
if(iso=="181Hg" || iso=="182Hg" || iso=="183Hg" || iso=="184Hg"){
chain->Add("../../../Outputs/Analysis/182Hg.root");
}
if(iso=="178Hg" || iso=="179Hg" || iso=="180Hg" || iso=="181Hg"){
chain->Add("../../../Outputs/Analysis/180Hg_1.root");
chain->Add("../../../Outputs/Analysis/180Hg_2.root");
}
//*** Au isotopes ***//
if(iso=="177Au" || iso=="178Au" || iso=="179Au"){
chain->Add("../../../Outputs/Analysis/175Pt.root");
}
if(iso=="180Au" || iso=="181Au" || iso=="182Au"){
chain->Add("../../../Outputs/Analysis/189Pb.root");
}
if(iso=="179Au" || iso=="180Au" || iso=="181Au" || iso=="182Au"){
chain->Add("../../../Outputs/Analysis/187Pb.root");
}
if(iso=="180Au" || iso=="181Au" || iso=="182Au"){
chain->Add("../../../Outputs/Analysis/184Hg.root");
}
if(iso=="178Au" || iso=="179Au" || iso=="180Au" || iso=="181Au" || iso=="182Au"){
chain->Add("../../../Outputs/Analysis/182Hg.root");
}
if(iso=="177Au" || iso=="178Au" || iso=="179Au" || iso=="180Au"){
chain->Add("../../../Outputs/Analysis/180Hg_1.root");
chain->Add("../../../Outputs/Analysis/180Hg_2.root");
}
//*** Pt isotopes ***//
if(iso=="173Pt" || iso=="174Pt" || iso=="175Pt" || iso=="176Pt" || iso=="177Pt" || iso=="178Pt"){
chain->Add("../../../Outputs/Analysis/175Pt.root");
}
if(iso=="174Pt" || iso=="175Pt" || iso=="176Pt" || iso=="177Pt" || iso=="178Pt"){
chain->Add("../../../Outputs/Analysis/180Hg_1.root");
chain->Add("../../../Outputs/Analysis/180Hg_2.root");
}
//*** Ir isotopes ***//
if(iso=="172Ir" || iso=="173Ir" || iso=="174Ir" || iso=="175Ir" || iso=="176It"){
chain->Add("../../../Outputs/Analysis/175Pt.root");
chain->Add("../../../Outputs/Analysis/180Hg_1.root");
chain->Add("../../../Outputs/Analysis/180Hg_2.root");
}
cout << "Number of entries: " << chain->GetEntries() << endl;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment