Skip to content
Snippets Groups Projects
Commit 708d8785 authored by Cyril Lenain's avatar Cyril Lenain :surfer_tone3:
Browse files

merging Vendeta macros

parents 968cd667 1f5a4966
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 123 additions and 399 deletions
NPLib/Detectors/Sofia/GladFieldMap.cxx
+ 84
− 63
View file @ 708d8785
......@@ -48,9 +48,9 @@ GladFieldMap::GladFieldMap() {
m_bin = 50;
m_Current = 2135.;
m_Scale = m_Current/3583.81;
m_Glad_Entrance = TVector3(0,0,2774.);
m_Glad_TurningPoint = TVector3(0,0,2774.+1654);
m_Tilt = 14.*deg;
m_Glad_Entrance = TVector3(0,0,-1.135*m);
m_Glad_TurningPoint = TVector3(0,0,0);
m_Tilt = -14.*deg;
m_B = m_Scale*m_Bmax;
for(int i=0; i<81; i++){
for(int j=0; j<41; j++){
......@@ -71,9 +71,9 @@ GladFieldMap::GladFieldMap() {
m_Nx= 0;
m_Ny= 0;
m_Nz= 0;
m_CentralTheta = 20.*deg;
m_MWPC3_POS = TVector3(-1436.,0,7852);
m_Angle_MWPC3 = 20.*deg;
m_CentralTheta = -20.*deg;
m_MWPC3_POS = TVector3(-1436.,0,3402);
m_Angle_MWPC3 = -20.*deg;
m_R_MWPC3 = 4199.;
m_InitPos = TVector3(0,0,0);
......@@ -92,7 +92,7 @@ double GladFieldMap::Delta(const double* parameter){
static vector<TVector3> pos;
pos = Propagate(parameter[0], m_InitPos, m_InitDir, false);
//pos.back().RotateY(-m_Angle_MWPC3+m_Tilt);
//pos.back().RotateY(-m_CentralTheta+m_Tilt);
diff = pos.back() - m_FinalPos;
//cout << pos.back().X() << " " << pos.back().Z() << endl;
......@@ -103,23 +103,22 @@ double GladFieldMap::Delta(const double* parameter){
double GladFieldMap::FindBrho(TVector3 Pos_init, TVector3 Dir_init, TVector3 Pos_final){
if(!m_BrhoScan)
BrhoScan(1.,20,0.2,TVector3(0,0,0), TVector3(0,0,1));
BrhoScan(6,11,0.1,TVector3(0,0,-2500), TVector3(0,0,1));
m_InitPos = Pos_init;
m_InitDir = Dir_init;
m_FinalPos = Pos_final;
m_InitDir = m_InitDir.Unit();
//BrhoScan(2,15,1,m_InitPos, m_InitDir);
//BrhoScan(6,11,1,m_InitPos, m_InitDir);
static double param[1];
param[0] = m_BrhoScan->Eval(m_FinalPos.X());
//return param[0];
m_min->Clear();
m_min->SetPrecision(1e-6);
m_min->SetMaxFunctionCalls(1000);
m_min->SetLimitedVariable(0,"B",param[0],0.1,1,20);
m_min->SetLimitedVariable(0,"B",param[0],0.1,6,11);
m_min->Minimize();
return m_min->X()[0];
......@@ -141,7 +140,7 @@ TGraph* GladFieldMap::BrhoScan(double Brho_min, double Brho_max, double Brho_ste
//dir.RotateY(m_Tilt);
for(double Brho=Brho_min; Brho<Brho_max; Brho+=Brho_step){
vector<TVector3> vPos = Propagate(Brho, pos, dir, false);
//vPos.back().RotateY(-m_Angle_MWPC3);
//vPos.back().RotateY(m_CentralTheta);
m_BrhoScan->SetPoint(i++,vPos.back().X(),Brho);
//cout << vPos.back().X() << " " << Brho << endl;
......@@ -152,59 +151,87 @@ TGraph* GladFieldMap::BrhoScan(double Brho_min, double Brho_max, double Brho_ste
return m_BrhoScan;
}
//////////////////////////////////////////////////////////////////////
double GladFieldMap::GetFlightPath(TVector3 vStart, double Brho, TVector3 Pos, TVector3 Dir){
double FlightPath = 0;
vector<TVector3> track;
track = Propagate(Brho, Pos, Dir, true);
FlightPath += (Pos - vStart).Mag();
unsigned int vsize = track.size();
for(unsigned int i=0; i<vsize-1; i++){
TVector3 point1 = track[i];
TVector3 point2 = track[i+1];
TVector3 v12 = point2 - point1;
FlightPath += v12.Mag();
}
return FlightPath;
}
//////////////////////////////////////////////////////////////////////
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 +244,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 +260,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 +297,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 +356,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 +377,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 +395,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;
......
NPLib/Detectors/Sofia/GladFieldMap.h
+ 6
− 1
View file @ 708d8785
......@@ -98,6 +98,8 @@ class GladFieldMap{
m_Scale = m_Current/3583.81;
m_B = 2.2*m_Scale;
}
void SetBin(double val) {m_bin = val;}
void SetTimeStep(double val) {m_dt = 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 +124,10 @@ class GladFieldMap{
double GetZmin() {return m_z_min;}
double GetZmax() {return m_z_max;}
double GetCentralTheta() {return m_CentralTheta;}
double GetBin() {return m_bin;}
double GetTimeStep() {return m_dt;}
TVector3 Get_MWPC3_Position() {return m_MWPC3_POS;}
public:
void LoadMap(string filename);
vector<double> InterpolateB(const vector<double>& pos);
......@@ -135,6 +139,7 @@ class GladFieldMap{
TGraph* BrhoScan(double Brho_min, double Brho_max, double Brho_step, TVector3 pos, TVector3 dir);
TVector3 CalculateIntersectionPoint(vector<TVector3> vPos);
vector<TVector3> Propagate(double Brho, TVector3 Pos, TVector3 Dir, bool store);
double GetFlightPath(TVector3 vStart, double Brho, TVector3 Pos, TVector3 Dir);
TVector3 PropagateToMWPC(TVector3 pos, TVector3 dir);
void func(NPL::Particle& N, TVector3 Pos, TVector3 Imp, TVector3& xk, TVector3& pk);
double FindBrho(TVector3 Pos_init, TVector3 Dir_init, TVector3 Pos_final);
......
NPSimulation/Detectors/Vendeta/Vendeta.cc
+ 1
− 1
View file @ 708d8785
......@@ -60,7 +60,7 @@ using namespace CLHEP;
namespace Vendeta_NS{
// Energy and time Resolution
const double EnergyThreshold = 0.1*MeV;
const double ResoTime = 4.5*ns ;
const double ResoTime = 0.6*ns ;
const double ResoEnergy = 1.0*MeV ;
const double Thickness = 51.*mm ;
const double Radius = 127./2*mm ;
......
Projects/Sofia/macro/ShowResultSimu.C
+ 2
− 1
View file @ 708d8785
......@@ -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());
}
......
Projects/Vendeta/macro/FillTOFHisto.C deleted 100644 → 0
+ 0
− 130
View file @ 968cd667
// c++ includes
#include <vector>
#include <iostream>
#include <iomanip>
#include <stdlib.h>
// ROOT includes
#include "TChain.h"
#include "TFile.h"
#include "TObject.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TGraph.h"
using namespace std;
int NumberOfDetectors= 72;
int NumberOfAnodes= 11;
int nentries=1e6;
/////////////////////////////////////
int main(int argc, char** argv){
int run = atof(argv[1]);
auto chain = new TChain("PhysicsTree");
chain->Add(Form("/home/cyril/Workflow/nptool/Outputs/Analysis/run%i/run%i.root",run,run));
nentries = chain->GetEntries();
/* nentries = 50000000; */
cout << "Number of entries: " << nentries << endl;
TFile* ofile = new TFile(Form("histos_ToF/histo_tof_file_run%i.root",run),"recreate");
TH2F* hLG_LSci[11];
TH2F* hHG_LSci[11];
TH1F* hLG[791];
TH1F* hHG[791];
TGraph *gFlyPath = new TGraph();
gFlyPath->SetTitle(" ; Vendeta-Anode Index ; fly length (mm);");
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>* LG_FlyPath = 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>();
vector<double>* HG_FlyPath = new vector<double>();
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("LG_FlyPath",&LG_FlyPath);
chain->SetBranchAddress("HG_Tof",&HG_Tof);
chain->SetBranchAddress("HG_ID",&HG_ID);
chain->SetBranchAddress("HG_Q1",&HG_Q1);
chain->SetBranchAddress("HG_Q2",&HG_Q2);
chain->SetBranchAddress("HG_FlyPath",&HG_FlyPath);
for(int j=0; j<NumberOfAnodes; j++){
TString histo_name = Form("hLG_LSci_Anode%i",j+1);
hLG_LSci[j] = new TH2F(histo_name,histo_name,72,1,73,2000,-100,300);
histo_name = Form("hHG_LSci_Anode%i",j+1);
hHG_LSci[j] = new TH2F(histo_name,histo_name,72,1,73,2000,-100,300);
for(int i=0; i<NumberOfDetectors; i++){
int index = j + i*NumberOfAnodes;
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,2000,-100,300);
}
}
for(int i=0; i<nentries; i++){
chain->GetEntry(i);
if(i%100000==0){
cout << setprecision(2) << "\033[34m\r Processing run "<< run <<" ..." << (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));
hLG_LSci[anode-1]->Fill(LG_ID->at(j),LG_Tof->at(j));
gFlyPath->SetPoint(index_LG,index_LG,LG_FlyPath->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));
hHG_LSci[anode-1]->Fill(HG_ID->at(j),HG_Tof->at(j));
}
}
}
}
ofile->WriteObject(gFlyPath,"gFlyPath");
ofile->Write();
ofile->Close();
}
Projects/Vendeta/macro/FitTofGammaPeak.C deleted 100644 → 0
+ 0
− 186
View file @ 968cd667
// c++ includes
#include <vector>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <stdlib.h>
// ROOT includes
#include "TChain.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TObject.h"
#include "TSpectrum.h"
#include "TGraph.h"
#include "TF1.h"
#include "TH1F.h"
#include "TH2F.h"
using namespace std;
int NumberOfDetectors=72;
int NumberOfAnodes=11;
int m_NumberOfGammaPeak=1;
double PosGammaPeak = 3.2;
double c = 299.792458; // mm/ns
double sigma_anode[11], sigma_anode_HG[11];
bool Finder(TH1F* h, Double_t *mean, Double_t *sigma);
/////////////////////////////////////////////////////
int main(int argc, char** argv){
int run = atof(argv[1]);
TString Append = argv[2];
auto ifile = new TFile(Form("histos_ToF/histo_tof_file_run%d"+Append+".root",run));
auto gFlyPath = (TGraph*)ifile->Get("gFlyPath");
TCanvas *cLG[11];
TCanvas *cHG[11];
ofstream ofile;
ofile.open(Form("ToF_calibs/Vendeta_Time_run%d.cal",run));
TFile* orootfile = new TFile(Form("histos_ToF/histo_tof_fitted_run%d"+Append+".root",run),"recreate");
for(int i=0; i<NumberOfAnodes; i++){
sigma_anode[i]=0;
sigma_anode_HG[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();
int countLG=0;
int countHG=0;
for(int i=0; i<NumberOfDetectors; i++){
for(int j=0; j<NumberOfAnodes; j++){
// LG //
TString histo_name = Form("hLG_Det%i_Anode%i",i+1,j+1);
TH1F* h = (TH1F*) ifile->FindObjectAny(histo_name);
//cLG[j]->cd(i+1);
h->Draw();
int index = j + i*11;
double FlyPath = gFlyPath->GetPointY(index);
PosGammaPeak = FlyPath/c ;
mean[0] = 0;
sigma[0] = 0;
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]);
sigma_anode[j] += sigma[0];
countLG++;
h->Write();
}
else{
ofile << LG_token << " 0" << endl;
}
// HG //
histo_name = Form("hHG_Det%i_Anode%i",i+1,j+1);
h = (TH1F*) ifile->FindObjectAny(histo_name);
//cHG[j]->cd(i+1);
h->Draw();
mean[0] = 0;
sigma[0] = 0;
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]);
sigma_anode_HG[j] += sigma[0];
countHG++;
h->Write();
}
else{
ofile << HG_token << " 0" << endl;
}
}
}
TCanvas* c1 = new TCanvas("Sigma","Sigma",1200,600);
c1->Divide(2,1);
gSigma_LG->SetMarkerStyle(8);
gSigma_HG->SetMarkerStyle(8);
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();
orootfile->Write();
orootfile->Close();
double totLG = 0, totHG=0;
for(int i=0; i<NumberOfAnodes; i++){
cout << "Anode= " << i+1 << " | LG: " << sigma_anode[i]/NumberOfDetectors << " HG: " << sigma_anode_HG[i]/NumberOfDetectors << endl;
totLG += sigma_anode[i];
totHG += sigma_anode_HG[i];
}
totLG = totLG/NumberOfDetectors/11.;
totHG = totHG/NumberOfDetectors/11.;
cout << "<sigma> LG: "<< totLG << " HG: "<<totHG << endl;
}
/////////////////////////////////////////////////////
bool Finder(TH1F* h, Double_t *mean, Double_t *sigma){
Double_t resolsig = 2;
Double_t resolsigTSpec = 1;
Double_t seuil = 0.3;
TSpectrum* s = new TSpectrum(m_NumberOfGammaPeak,resolsigTSpec);
Int_t nfound = s->Search(h,resolsig,"new",seuil);
Double_t *xpeaks = s->GetPositionX();
Double_t linf=0;
Double_t lsup=0;
if(nfound == m_NumberOfGammaPeak){
cout << "Gamma Peak Found" << endl;
for(int i=0; i<nfound; i++){
linf = xpeaks[i]-1.1;
lsup = xpeaks[i]+0.8;
/* if(anode==9){ */
/* lsup =xpeaks[i]+0.8 */
/* } */
TF1* gauss = new TF1("gaus","gaus",linf,lsup);
h->Fit(gauss,"RQ");
mean[i] = gauss->GetParameter(1);
sigma[i] = gauss->GetParameter(2);
}
for(int i=0; i<nfound; i++){
linf = xpeaks[i]-1.3*sigma[i];
lsup = xpeaks[i]+1*sigma[i];
TF1* gauss = new TF1("gaus","gaus",linf,lsup);
h->Fit(gauss,"RQ");
mean[i] = gauss->GetParameter(1);
sigma[i] = gauss->GetParameter(2);
}
return true;
}
else if(nfound != m_NumberOfGammaPeak){
cout << "Warning. Number of peak different of " << m_NumberOfGammaPeak << " !! / nfound = " << nfound << endl;
return false;
}
}
Projects/s455/Analysis.cxx
+ 30
− 17
View file @ 708d8785
......@@ -103,14 +103,20 @@ void Analysis::Init(){
m_GladField = new GladFieldMap();
m_GladField->SetCurrent(2135.);
m_GladField->SetGladEntrance(0, 0.02*m, 2.774*m + 0.5405*m);
m_GladField->SetGladTurningPoint(0, 0.02*m, 2.774*m + 0.5405*m + 1.135*m);
m_GladField->SetGladTiltAngle(14.*deg);
m_GladField->LoadMap("GladFieldMap.dat");
m_GladField->SetCentralTheta(20.*deg);
//m_GladField->SetGladEntrance(0, 0.02*m, 2.774*m + 0.5405*m);
m_GladField->SetGladEntrance(0, 0, -1.135*m);
//m_GladField->SetGladTurningPoint(0, 0.02*m, 2.774*m + 0.5405*m + 1.135*m);
m_GladField->SetGladTurningPoint(0, 0, 0);
m_GladField->SetGladTiltAngle(-14.*deg);
m_GladField->LoadMap("GladFieldMap_50mm.dat");
m_GladField->SetBin(50);
m_GladField->SetTimeStep(0.8);
m_GladField->SetCentralTheta(-20.*deg);
double Z_MWPC3 = 7.852*m;
double X_MWPC3 = -(Z_MWPC3 - m_GladField->GetGladTurningPoint().Z())*tan(m_GladField->GetCentralTheta());
//double Z_MWPC3 = 7.852*m;
//double Z_MWPC3 = 3.402*m;
double Z_MWPC3 = 3.65*m;
double X_MWPC3 = (Z_MWPC3 - m_GladField->GetGladTurningPoint().Z())*tan(m_GladField->GetCentralTheta());
m_GladField->Set_MWPC3_Position(X_MWPC3,0,Z_MWPC3);
InitParameter();
......@@ -586,17 +592,21 @@ void Analysis::FissionFragmentAnalysis(){
// *** Calculation Theta_out *** //
double Theta0 = m_GladField->GetCentralTheta();
double DistanceStartToG = 2.774*m + 0.5405*m +1.135*m;
double DistanceStartToA = 2.3155*m;
double DistanceStartToMW2 = 2.651*m;
double Theta0 = 20.*deg;//m_GladField->GetCentralTheta();
double XA = 0;
double ZA = 2315.5;
double ZA = DistanceStartToA - DistanceStartToG;
double ZG = m_GladField->GetGladTurningPoint().Z();
double ZMW3 = m_GladField->Get_MWPC3_Position().Z();
double XMW3 = m_GladField->Get_MWPC3_Position().X();
double ZMW2 = 2651;
double ZMW2 = DistanceStartToMW2 - DistanceStartToG;
double X3lab = 0;
double Z3lab = 0;;
double Tilt = m_GladField->GetGladTiltAngle();;//14.*deg;
double Tilt = 14.*deg;//;
TVector3 vZ = TVector3(0,0,1);
TVector3 vStart = TVector3(0,0,-DistanceStartToG);
double XC = 0;
double ZC = 0;
double XB = 0;
......@@ -613,18 +623,19 @@ void Analysis::FissionFragmentAnalysis(){
XC = (XA+(ZG-ZA)*tan(TofHit[i].theta_in)) / (1-tan(Tilt)*tan(TofHit[i].theta_in));
ZC = ZG + XC*tan(Tilt);
TVector3 vC = TVector3(XC,0,ZC);
TofHit[i].xc = XC;
TofHit[i].yc = (ZC/8592.)*TofHit[i].y;
TofHit[i].zc = ZC;
int ix, iy;
ix = (int)(TofHit[0].xc - m_GladField->GetXmin())/50;
iy = (int)(TofHit[0].yc -20. - m_GladField->GetYmin())/50;
ix = (int)(TofHit[0].xc - m_GladField->GetXmin())/m_GladField->GetBin();
iy = (int)(TofHit[0].yc - m_GladField->GetYmin())/m_GladField->GetBin();
TofHit[i].Leff = m_GladField->GetLeff(ix,iy);
X3lab = TofHit[i].x3*cos(Theta0) + XMW3;
Z3lab = TofHit[i].x3*sin(Theta0) + ZMW3;
TVector3 vE = TVector3(X3lab, TofHit[i].y, Z3lab);
TVector3 dir = TVector3(sin(TofHit[i].theta_in), 0, cos(TofHit[i].theta_in));
TofHit[i].x3lab = X3lab;
TofHit[i].z3lab = Z3lab;
TVector3 vOut = TVector3(X3lab-XC,0,Z3lab-ZC);
......@@ -632,7 +643,8 @@ void Analysis::FissionFragmentAnalysis(){
TofHit[i].theta_out = angle;
TofHit[i].psi = TofHit[i].theta_in - TofHit[i].theta_out;
TofHit[i].rho = TofHit[i].Leff/(2*sin(0.5*TofHit[i].psi)*cos(Tilt-0.5*TofHit[i].psi));
TofHit[i].Brho = MagB*TofHit[i].rho;
//TofHit[i].Brho = MagB*TofHit[i].rho;
TofHit[i].Brho = m_GladField->FindBrho(vA,dir,vE);
// *** Extrapolate to B position *** //
double ZI = ZG - TofHit[i].Leff/(2*cos(Tilt));
......@@ -661,13 +673,14 @@ void Analysis::FissionFragmentAnalysis(){
TVector3 v3lab = TVector3(X3lab,0,Z3lab);
TVector3 v1 = TVector3(XB,0,ZB);
TVector3 v1 = TVector3(XB,0,ZB)-vStart;
TVector3 v3 = TVector3(X3lab-XD,0,Z3lab-ZD);
TofHit[i].omega = abs(2.*asin(sqrt(pow(XD-XB,2) + pow(ZD-ZB,2))/(2*TofHit[i].rho)));
double Path1 = v1.Mag();
double Path2 = TofHit[i].rho*TofHit[i].omega;
double Path3 = v3.Mag();
double PathLength = Path1 + Path2 + Path3 + 740. + 50.;
//double PathLength = Path1 + Path2 + Path3 + 740. + 50.;
double PathLength = m_GladField->GetFlightPath(vStart, TofHit[i].Brho, vA, dir) + 740. + 50;
PathLength = PathLength/1000.;
TofHit[i].flight_path = PathLength;
......
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