diff --git a/NPLib/Detectors/Sofia/GladFieldMap.cxx b/NPLib/Detectors/Sofia/GladFieldMap.cxx
index 0ffcc2ca3a486a1d2d7d029da33a91e9646d5187..8c8719def0ef8f882e33a641c574f885ed257b83 100644
--- a/NPLib/Detectors/Sofia/GladFieldMap.cxx
+++ b/NPLib/Detectors/Sofia/GladFieldMap.cxx
@@ -155,56 +155,62 @@ TGraph* GladFieldMap::BrhoScan(double Brho_min, double Brho_max, double Brho_ste
//////////////////////////////////////////////////////////////////////
TVector3 GladFieldMap::PropagateToMWPC(TVector3 pos, TVector3 dir){
// go to MWPC3 frame reference
- //pos.RotateY(-m_Angle_MWPC3);
- //dir.RotateY(-m_Angle_MWPC3);
+ pos.RotateY(-m_CentralTheta);
+ dir.RotateY(-m_CentralTheta);
double deltaZ = m_MWPC3_POS.Z() - pos.Z();
dir*=deltaZ/dir.Z();
pos+=dir;
pos.SetX(pos.X());
- //pos.RotateY(m_Angle_MWPC3);
+ pos.RotateY(m_CentralTheta);
return pos;
}
//////////////////////////////////////////////////////////////////////
vector<TVector3> GladFieldMap::Propagate(double Brho, TVector3 Pos, TVector3 Dir, bool store){
- //Pos.RotateY(m_Tilt);
- //Dir.RotateY(m_Tilt);
+ Pos.RotateY(m_Tilt);
+ Dir.RotateY(m_Tilt);
+ Dir = Dir.Unit();
+
static NPL::Particle N("90Zr");
N.SetBrho(Brho);
-
+
// track result
static std::vector<TVector3> track;
track.clear();
// starting point of the track
if(store){
- //Pos.RotateY(-m_Tilt);
+ Pos.RotateY(-m_Tilt);
track.push_back(Pos);
- //Pos.RotateY(-m_Tilt);
+ Pos.RotateY(m_Tilt);
}
- //static double r;
- //r = sqrt(Pos.X()*Pos.X() + Pos.Z()*Pos.Z());
+ Dir = Dir.Unit();
+ static double r;
+ r = sqrt(Pos.X()*Pos.X() + Pos.Z()*Pos.Z());
+ // number of step taken
+ static unsigned int count;
+ count = 0;
// Propagate to m_R_max with one line
- /*while(r>m_R_max && count<m_Limit){
+ while(r>m_R_max && count<m_Limit){
Pos+=(r-m_R_max)/cos(Dir.Theta())*Dir.Unit();
r = 1.01*sqrt(Pos.X()*Pos.X() + Pos.Z()*Pos.Z());
- }
+ }
- if(r<=m_R_max){ // success
+ if(r<=m_R_max){ // success
if(store){
- //Pos.RotateY(-m_Tilt);
- track.push_back(Pos);
- //Pos.RotateY(m_Tilt);
- }
+ Pos.RotateY(-m_Tilt);
+ track.push_back(Pos);
+ Pos.RotateY(m_Tilt);
+ }
}
else{
- cout << "Fail" << endl;
- return track;
- }*/
+ cout << "Fail" << endl;
+ return track;
+ }
static TVector3 xk1, xk2, xk3, xk4;
static TVector3 pk1, pk2, pk3, pk4;
@@ -217,15 +223,13 @@ vector<TVector3> GladFieldMap::Propagate(double Brho, TVector3 Pos, TVector3 Dir
E = T + M;
P = sqrt(T*T + 2*M*T)/c_light;
- Dir = Dir.Unit();
px = P*Dir.X();
py = P*Dir.Y();
pz = P*Dir.Z();
Imp = TVector3(px,py,pz);
- unsigned int count=0;
- while(Pos.Z()<=m_Zmax && count<m_Limit){
+ while(Pos.Z()<=m_R_max && count<m_Limit){
func(N, Pos, Imp, xk1, pk1);
func(N, Pos + (m_dt/2)*xk1, Imp + (m_dt/2)*pk1, xk2, pk2);
func(N, Pos + (m_dt/2)*xk2, Imp + (m_dt/2)*pk2, xk3, pk3);
@@ -235,17 +239,17 @@ vector<TVector3> GladFieldMap::Propagate(double Brho, TVector3 Pos, TVector3 Dir
Imp += (m_dt/6)*(pk1 + 2*pk2 + 2*pk3 + pk4);
if(store){
- //Pos.RotateY(-m_Tilt);
+ Pos.RotateY(-m_Tilt);
track.push_back(Pos);
- //Pos.RotateY(m_Tilt);
+ Pos.RotateY(m_Tilt);
}
- //r = sqrt(Pos.X()*Pos.X() + Pos.Z()*Pos.Z());
+ r = sqrt(Pos.X()*Pos.X() + Pos.Z()*Pos.Z());
count++;
}
Imp=Imp.Unit();
Pos = PropagateToMWPC(Pos,Imp);
- //Pos.RotateY(-m_Tilt);
+ Pos.RotateY(-m_Tilt);
track.push_back(Pos);
return track;
@@ -272,10 +276,10 @@ void GladFieldMap::func(NPL::Particle& N, TVector3 Pos, TVector3 Imp, TVector3&
xk.SetY(vy);
xk.SetZ(vz);
- /*B = InterpolateB(Pos);
- Bx = B[0];
- By = B[1];
- Bz = B[2];*/
+ //B = InterpolateB(Pos);
+ //Bx = B[0];
+ //By = B[1];
+ //Bz = B[2];
Bx = GetB(Pos,"X");
By = GetB(Pos,"Y");
@@ -331,12 +335,12 @@ void GladFieldMap::LoadMap(string filename) {
ifile >> m_y_min >> m_y_max >> m_Ny;
ifile >> m_z_min >> m_z_max >> m_Nz;
- m_x_min = m_x_min*10;
- m_x_max = m_x_max*10;
- m_y_min = m_y_min*10;
- m_y_max = m_y_max*10;
- m_z_min = m_z_min*10;
- m_z_max = m_z_max*10;
+ m_x_min = m_x_min;
+ m_x_max = m_x_max;
+ m_y_min = m_y_min;
+ m_y_max = m_y_max;
+ m_z_min = m_z_min;
+ m_z_max = m_z_max;
unsigned int count=0;
int index = 0;
@@ -352,18 +356,14 @@ void GladFieldMap::LoadMap(string filename) {
//cout << x << " " << y << " " << z << " " << Bx << " " << By << " " << Bz << endl;
- x = x*10;
- y = y*10;
- z = z*10;
-
- //m_Leff[ix][iy] += abs(By)*m_bin;
+ //m_Leff[ix][iy] += By*m_bin;
// Need to fill this TGraph before scaling the field to get the proper Leff //
gBy->SetPoint(iz,z,By);
x += m_Glad_Entrance.X();
y += m_Glad_Entrance.Y();
-
- z = z + x*sin(m_Tilt);
+
+ //z = z + x*sin(m_Tilt);
z += m_Glad_Entrance.Z();
Bx *= -m_Scale;
@@ -374,18 +374,18 @@ void GladFieldMap::LoadMap(string filename) {
m_By.push_back(By*tesla);
m_Bz.push_back(Bz*tesla);
- /*vector<double> p = {x,y,z};
- Bx*=tesla;
- By*=tesla;
- Bz*=tesla;
- vector<double> B = {Bx,By,Bz};
- m_field[p] = B;*/
+ vector<double> p = {x,y,z};
+ Bx*=tesla;
+ By*=tesla;
+ Bz*=tesla;
+ vector<double> B = {Bx,By,Bz};
+ m_field[p] = B;
}
m_Leff[ix][iy] = gBy->Integral()/m_Bmax;
}
}
- m_R_max = m_z_max;
+ m_R_max = m_x_max;
cout << endl;
cout << "///////// ASCII file loaded"<< endl;
cout << "m_field size= " << m_By.size() << endl;
diff --git a/NPLib/Detectors/Sofia/GladFieldMap.h b/NPLib/Detectors/Sofia/GladFieldMap.h
index 0071bee2c2075f156a1a370f0c0975b441948dce..d9a956537e75ba270b0bbe714be286451810d2fe 100644
--- a/NPLib/Detectors/Sofia/GladFieldMap.h
+++ b/NPLib/Detectors/Sofia/GladFieldMap.h
@@ -98,6 +98,7 @@ class GladFieldMap{
m_Scale = m_Current/3583.81;
m_B = 2.2*m_Scale;
}
+ void SetBin(double val) {m_bin = val;}
void SetCentralTheta(double val) {m_CentralTheta = val;}
void Set_MWPC3_Position(double x, double y, double z) {m_MWPC3_POS = TVector3(x,y,z);}
@@ -122,8 +123,9 @@ class GladFieldMap{
double GetZmin() {return m_z_min;}
double GetZmax() {return m_z_max;}
double GetCentralTheta() {return m_CentralTheta;}
+ double GetBin() {return m_bin;}
TVector3 Get_MWPC3_Position() {return m_MWPC3_POS;}
-
+
public:
void LoadMap(string filename);
vector<double> InterpolateB(const vector<double>& pos);
diff --git a/Projects/Sofia/macro/ShowResultSimu.C b/Projects/Sofia/macro/ShowResultSimu.C
index a2ccbf897a539d1772c46eee422e03bcd2d018e9..bd3dd82b159aa5802bb68d9c2095292a50bab9b6 100644
--- a/Projects/Sofia/macro/ShowResultSimu.C
+++ b/Projects/Sofia/macro/ShowResultSimu.C
@@ -4,7 +4,8 @@ TChain* chain=NULL;
void LoadRootFile(string nucleus){
chain = new TChain("SimulatedTree");
- string rootfilename = "../../../Outputs/Simulation/sofia_simu_" + nucleus + "_*";
+ //string rootfilename = "../../../Outputs/Simulation/sofia_simu_" + nucleus + "_*";
+ string rootfilename = "../../../Outputs/Simulation/sofia_simu_1.root";
chain->Add(rootfilename.c_str());
}
diff --git a/Projects/Vendeta/macro/FillTOFHisto.C b/Projects/Vendeta/macro/FillTOFHisto.C
index 558383cc70048659cedcb3ecf36b9062859766d7..4e944ac5aebc1eacb38636866ef5e760fb99a708 100644
--- a/Projects/Vendeta/macro/FillTOFHisto.C
+++ b/Projects/Vendeta/macro/FillTOFHisto.C
@@ -6,99 +6,99 @@ int nentries=1e6;
int run=86;
/////////////////////////////////////
void LoadRootFile(){
- chain = new TChain("PhysicsTree");
- chain->Add(Form("/home/faster/nptool/Outputs/Analysis/run%i.root",run));
- //chain->Add("/home/faster/nptool/Outputs/Analysis/test_sampler_qdc_cf_1.root");
- //chain->Add("/home/faster/nptool/Outputs/Analysis/test_sampler_qdc_cf_2.root");
- //chain->Add("/home/faster/nptool/Outputs/Analysis/test_sampler_qdc_cf_3.root");
- //chain->Add("/home/faster/nptool/Outputs/Analysis/test_sampler_qdc_cf_4.root");
+ chain = new TChain("PhysicsTree");
+ chain->Add(Form("/home/faster/nptool/Outputs/Analysis/run%i.root",run));
+ //chain->Add("/home/faster/nptool/Outputs/Analysis/test_sampler_qdc_cf_1.root");
+ //chain->Add("/home/faster/nptool/Outputs/Analysis/test_sampler_qdc_cf_2.root");
+ //chain->Add("/home/faster/nptool/Outputs/Analysis/test_sampler_qdc_cf_3.root");
+ //chain->Add("/home/faster/nptool/Outputs/Analysis/test_sampler_qdc_cf_4.root");
}
/////////////////////////////////////
void FillTOFHisto(){
- LoadRootFile();
- nentries = chain->GetEntries();
- cout << "Number of entries: " << nentries << endl;
-
- TFile* ofile = new TFile(Form("histo_tof_file_run%i.root",run),"recreate");
- TH1F* hLG[791];
- TH1F* hHG[791];
-
- vector<double>* FC_Q1 = new vector<double>();
- vector<int>* FC_Anode_ID = new vector<int>();
- vector<bool>* FC_FakeFission = new vector<bool>();
-
- vector<double>* LG_Tof = new vector<double>();
- vector<int>* LG_ID = new vector<int>();
- vector<double>* LG_Q1 = new vector<double>();
- vector<double>* LG_Q2 = new vector<double>();
-
- vector<double>* HG_Tof = new vector<double>();
- vector<int>* HG_ID = new vector<int>();
- vector<double>* HG_Q1 = new vector<double>();
- vector<double>* HG_Q2 = new vector<double>();
-
- TFissionChamberPhysics* FC = new TFissionChamberPhysics();
- chain->SetBranchAddress("FissionChamber",&FC);
-
- chain->SetBranchAddress("FC_Q1",&FC_Q1);
- chain->SetBranchAddress("FC_Anode_ID",&FC_Anode_ID);
- chain->SetBranchAddress("FC_FakeFission",&FC_FakeFission);
- chain->SetBranchAddress("LG_Tof",&LG_Tof);
- chain->SetBranchAddress("LG_ID",&LG_ID);
- chain->SetBranchAddress("LG_Q1",&LG_Q1);
- chain->SetBranchAddress("LG_Q2",&LG_Q2);
- chain->SetBranchAddress("HG_Tof",&HG_Tof);
- chain->SetBranchAddress("HG_ID",&HG_ID);
- chain->SetBranchAddress("HG_Q1",&HG_Q1);
- chain->SetBranchAddress("HG_Q2",&HG_Q2);
-
- for(int i=0; i<NumberOfDetectors; i++){
- for(int j=0; j<NumberOfAnodes; j++){
- int index = j + i*NumberOfAnodes;
- TString histo_name = Form("hLG_Det%i_Anode%i",i+1,j+1);
- hLG[index] = new TH1F(histo_name,histo_name,2000,-100,300);
-
- histo_name = Form("hHG_Det%i_Anode%i",i+1,j+1);
- hHG[index] = new TH1F(histo_name,histo_name,2500,0,500);
- }
- }
- for(int i=0; i<nentries; i++){
- chain->GetEntry(i);
-
- if(i%100000==0){
- cout << setprecision(2) << "\033[34m\r Processing tree..." << (double)i/nentries*100 << "\% done" << flush;
- }
-
- if(FC_Anode_ID->size()>0){
- bool Fake = FC_FakeFission->at(0);
- int anode = FC_Anode_ID->at(0);
-
-
- int mysize = LG_Tof->size();
- for(int j=0; j<mysize; j++){
- // LG //
- int index_LG = (anode-1) + (LG_ID->at(j)-1)*NumberOfAnodes;
- double PSD = LG_Q2->at(j)/LG_Q1->at(j);
- if(LG_ID->at(j)>0 && anode>0 && Fake==0 && LG_Q1->at(j)>2000){
- hLG[index_LG]->Fill(LG_Tof->at(j));
- }
- }
-
- mysize = HG_Tof->size();
- for(int j=0; j<mysize; j++){
- // HG //
- int index_HG = (anode-1) + (HG_ID->at(j)-1)*NumberOfAnodes;
- double PSD = HG_Q2->at(j)/HG_Q1->at(j);
- if(HG_ID->at(j)>0 && anode>0 && Fake==0 && HG_ID->size()==1){
- hHG[index_HG]->Fill(HG_Tof->at(j));
- }
- }
- }
- }
-
- ofile->Write();
- ofile->Close();
+ LoadRootFile();
+ nentries = chain->GetEntries();
+ cout << "Number of entries: " << nentries << endl;
+
+ TFile* ofile = new TFile(Form("histo_tof_file_run%i.root",run),"recreate");
+ TH1F* hLG[791];
+ TH1F* hHG[791];
+
+ vector<double>* FC_Q1 = new vector<double>();
+ vector<int>* FC_Anode_ID = new vector<int>();
+ vector<bool>* FC_FakeFission = new vector<bool>();
+
+ vector<double>* LG_Tof = new vector<double>();
+ vector<int>* LG_ID = new vector<int>();
+ vector<double>* LG_Q1 = new vector<double>();
+ vector<double>* LG_Q2 = new vector<double>();
+
+ vector<double>* HG_Tof = new vector<double>();
+ vector<int>* HG_ID = new vector<int>();
+ vector<double>* HG_Q1 = new vector<double>();
+ vector<double>* HG_Q2 = new vector<double>();
+
+ TFissionChamberPhysics* FC = new TFissionChamberPhysics();
+ chain->SetBranchAddress("FissionChamber",&FC);
+
+ chain->SetBranchAddress("FC_Q1",&FC_Q1);
+ chain->SetBranchAddress("FC_Anode_ID",&FC_Anode_ID);
+ chain->SetBranchAddress("FC_FakeFission",&FC_FakeFission);
+ chain->SetBranchAddress("LG_Tof",&LG_Tof);
+ chain->SetBranchAddress("LG_ID",&LG_ID);
+ chain->SetBranchAddress("LG_Q1",&LG_Q1);
+ chain->SetBranchAddress("LG_Q2",&LG_Q2);
+ chain->SetBranchAddress("HG_Tof",&HG_Tof);
+ chain->SetBranchAddress("HG_ID",&HG_ID);
+ chain->SetBranchAddress("HG_Q1",&HG_Q1);
+ chain->SetBranchAddress("HG_Q2",&HG_Q2);
+
+ for(int i=0; i<NumberOfDetectors; i++){
+ for(int j=0; j<NumberOfAnodes; j++){
+ int index = j + i*NumberOfAnodes;
+ TString histo_name = Form("hLG_Det%i_Anode%i",i+1,j+1);
+ hLG[index] = new TH1F(histo_name,histo_name,2000,-100,300);
+
+ histo_name = Form("hHG_Det%i_Anode%i",i+1,j+1);
+ hHG[index] = new TH1F(histo_name,histo_name,2500,0,500);
+ }
+ }
+ for(int i=0; i<nentries; i++){
+ chain->GetEntry(i);
+
+ if(i%100000==0){
+ cout << setprecision(2) << "\033[34m\r Processing tree..." << (double)i/nentries*100 << "\% done" << flush;
+ }
+
+ if(FC_Anode_ID->size()>0){
+ bool Fake = FC_FakeFission->at(0);
+ int anode = FC_Anode_ID->at(0);
+
+
+ int mysize = LG_Tof->size();
+ for(int j=0; j<mysize; j++){
+ // LG //
+ int index_LG = (anode-1) + (LG_ID->at(j)-1)*NumberOfAnodes;
+ double PSD = LG_Q2->at(j)/LG_Q1->at(j);
+ if(LG_ID->at(j)>0 && anode>0 && Fake==0 && LG_Q1->at(j)>2000){
+ hLG[index_LG]->Fill(LG_Tof->at(j));
+ }
+ }
+
+ mysize = HG_Tof->size();
+ for(int j=0; j<mysize; j++){
+ // HG //
+ int index_HG = (anode-1) + (HG_ID->at(j)-1)*NumberOfAnodes;
+ double PSD = HG_Q2->at(j)/HG_Q1->at(j);
+ if(HG_ID->at(j)>0 && anode>0 && Fake==0 && HG_ID->size()==1){
+ hHG[index_HG]->Fill(HG_Tof->at(j));
+ }
+ }
+ }
+ }
+
+ ofile->Write();
+ ofile->Close();
}
diff --git a/Projects/Vendeta/macro/FitTofGammaPeak.C b/Projects/Vendeta/macro/FitTofGammaPeak.C
index 62d4967acbb704150e0ec3efccbc886ea0751966..0d499792bec8687f968e674f50e14d36131f080f 100644
--- a/Projects/Vendeta/macro/FitTofGammaPeak.C
+++ b/Projects/Vendeta/macro/FitTofGammaPeak.C
@@ -25,19 +25,19 @@ void FitTofGammaPeak(){
ofstream ofile;
ofile.open(Form("Vendeta_Time_run%i.cal",run));
- TFile* orootfile = new TFile(Form("histo_tof_fitted_run%i.root",run),"recreate");
+ TFile* orootfile = new TFile(Form("histo_tof_fitted_run%i.root",run),"recreate");
for(int i=0; i<NumberOfAnodes; i++){
- sigma_anode[i]=0;
+ sigma_anode[i]=0;
}
Double_t* mean = new Double_t[1];
Double_t* sigma = new Double_t[1];
- TGraph* gSigma_LG = new TGraph();
- TGraph* gSigma_HG = new TGraph();
+ TGraph* gSigma_LG = new TGraph();
+ TGraph* gSigma_HG = new TGraph();
- int countLG=0;
- int countHG=0;
+ int countLG=0;
+ int countHG=0;
for(int i=0; i<NumberOfDetectors; i++){
for(int j=0; j<NumberOfAnodes; j++){
// LG //
@@ -52,12 +52,12 @@ void FitTofGammaPeak(){
TString LG_token = Form("Vendeta_DET%i_LG_ANODE%i_TIMEOFFSET",i+1,j+1);
if(Finder(h, mean, sigma)){
ofile << LG_token << " " << -mean[0]+PosGammaPeak << endl;
- gSigma_LG->SetPoint(countLG,countLG+1,sigma[0]);
+ gSigma_LG->SetPoint(countLG,countLG+1,sigma[0]);
- sigma_anode[j] += sigma[0];
+ sigma_anode[j] += sigma[0];
- countLG++;
- h->Write();
+ countLG++;
+ h->Write();
}
else{
ofile << LG_token << " 0" << endl;
@@ -75,9 +75,9 @@ void FitTofGammaPeak(){
TString HG_token = Form("Vendeta_DET%i_HG_ANODE%i_TIMEOFFSET",i+1,j+1);
if(Finder(h, mean, sigma)){
ofile << HG_token << " " << -mean[0]+PosGammaPeak << endl;
- gSigma_HG->SetPoint(countHG,countHG+1,sigma[0]);
- countHG++;
- h->Write();
+ gSigma_HG->SetPoint(countHG,countHG+1,sigma[0]);
+ countHG++;
+ h->Write();
}
else{
ofile << HG_token << " 0" << endl;
@@ -85,29 +85,29 @@ void FitTofGammaPeak(){
}
}
- TCanvas* c1 = new TCanvas("Sigma","Sigma",1200,600);
- c1->Divide(2,1);
+ TCanvas* c1 = new TCanvas("Sigma","Sigma",1200,600);
+ c1->Divide(2,1);
- gSigma_LG->SetMarkerStyle(8);
- gSigma_HG->SetMarkerStyle(8);
+ gSigma_LG->SetMarkerStyle(8);
+ gSigma_HG->SetMarkerStyle(8);
- c1->cd(1);
- gSigma_LG->Draw();
- c1->cd(2);
- gSigma_HG->Draw();
+ c1->cd(1);
+ gSigma_LG->Draw();
+ c1->cd(2);
+ gSigma_HG->Draw();
- gSigma_LG->SetName("sigma_LG");
- gSigma_HG->SetName("sigma_HG");
-
- gSigma_LG->Write();
- gSigma_HG->Write();
+ gSigma_LG->SetName("sigma_LG");
+ gSigma_HG->SetName("sigma_HG");
- orootfile->Write();
- orootfile->Close();
+ gSigma_LG->Write();
+ gSigma_HG->Write();
- for(int i=0; i<NumberOfAnodes; i++){
- cout << "Anode= " << i+1 << " / " << sigma_anode[i]/NumberOfDetectors << endl;
- }
+ orootfile->Write();
+ orootfile->Close();
+
+ for(int i=0; i<NumberOfAnodes; i++){
+ cout << "Anode= " << i+1 << " / " << sigma_anode[i]/NumberOfDetectors << endl;
+ }
}
@@ -139,7 +139,7 @@ bool Finder(TH1F* h, Double_t *mean, Double_t *sigma){
return true;
}
- else if(nfound != m_NumberOfGammaPeak){
+ else if(nfound != m_NumberOfGammaPeak){
cout << "Warning. Number of peak different of " << m_NumberOfGammaPeak << " !! / nfound = " << nfound << endl;
return false;
}