Commit 5a7988a4 authored by adrien-matta's avatar adrien-matta
Browse files

* Progress on FDC2

parent 323a41fd
/*****************************************************************************
* Copyright (C) 2009-2016 this file is part of the NPTool Project *
* Copyright (C) 2009-2020 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 *
......@@ -58,7 +58,7 @@ void TSamuraiFDC2Physics::BuildSimplePhysicalEvent(){
///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC2Physics::BuildPhysicalEvent(){
PreTreat();
//RemoveNoise();
RemoveNoise();
// Map[detector & plane angle, vector of spatial information]
static map<std::pair<unsigned int,double>, vector<double> > X ;
......@@ -68,7 +68,7 @@ void TSamuraiFDC2Physics::BuildPhysicalEvent(){
X.clear();Z.clear();R.clear();
unsigned int size = Detector.size();
for(unsigned int i = 0 ; i < size ; i++){
if( DriftLength[i] > DriftLowThreshold && DriftLength[i] < DriftUpThreshold){
if(DriftLength[i] > DriftLowThreshold && DriftLength[i] < DriftUpThreshold){
det = Detector[i];
layer = Layer[i];
wire = Wire[i];
......@@ -88,9 +88,9 @@ void TSamuraiFDC2Physics::BuildPhysicalEvent(){
VX0.clear();VX100.clear(),D.clear();
for(auto it = X.begin();it!=X.end();++it){
D[it->first]=m_reconstruction.BuildTrack2D(X[it->first],Z[it->first],R[it->first],X0,X100,a,b);
std::ofstream f("distance.txt", std::ios::app);
/* std::ofstream f("distance.txt", std::ios::app);
f<< D[it->first] << endl;
f.close();
f.close();*/
// very large a means track perpendicular to the chamber, what happen when there is pile up
if(abs(a)>1000)
PileUp++;
......@@ -341,9 +341,9 @@ void TSamuraiFDC2Physics::AddDC(string name, NPL::XmlParser& xml){
else if(sDir=="Y")
T=90*deg;
else if(sDir=="U")
T=30*deg;
else if(sDir=="V")
T=-30*deg;
else if(sDir=="V")
T=+30*deg;
else{
cout << "ERROR: Unknown layer orientation for Samurai FDC2"<< endl;
exit(1);
......
#ifndef TSAMURAIFDC2PHYSICS_H
#define TSAMURAIFDC2PHYSICS_H
/*****************************************************************************
* Copyright (C) 2009-2016 this file is part of the NPTool Project *
* Copyright (C) 2009-2020 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 *
......
/*****************************************************************************
* Copyright (C) 2009-2020 this file is part of the NPTool Project *
* Copyright (C) 2009-2020 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: Adrien Matta contact address: matta@lpccaen.in2p3.fr *
* Original Author: Adrien Matta contact address: matta@lpccaen.in2p3.fr *
* *
* Creation Date : October 2020 *
* Creation Date : October 2020 *
* Last update : *
*---------------------------------------------------------------------------*
* Decription: *
* This class hold eAGanil Treated data *
* This class hold eAGanil Treated data *
* *
*---------------------------------------------------------------------------*
* Comment: *
......@@ -203,30 +203,18 @@ void TeAGanilPhysics::Clear() {
///////////////////////////////////////////////////////////////////////////
void TeAGanilPhysics::ReadConfiguration(NPL::InputParser parser) {
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("eAGanil");
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithTokenAndValue("eAGanil","Spectrometer");
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << "//// " << blocks.size() << " detectors found " << endl;
vector<string> cart = {"POS","Shape"};
vector<string> sphe = {"R","Theta","Phi","Shape"};
vector<string> sphe = {"R","Theta","Phi","EntranceWidth","EntranceHeight","MomentumResolution"};
for(unsigned int i = 0 ; i < blocks.size() ; i++){
if(blocks[i]->HasTokenList(cart)){
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << endl << "//// eAGanil " << 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(blocks[i]->HasTokenList(sphe)){
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << endl << "//// eAGanil " << 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;
......
......@@ -44,7 +44,7 @@ R__LOAD_LIBRARY(libNPTrackReconstruction.dylib)
x.push_back(X);
DrawCircle(X,Z,R);
}
// add second track for pile up
/* // add second track for pile up
a = rand.Uniform(-10,10);
b = rand.Uniform(-100,100);
DrawTrack(a,b,3);
......@@ -59,52 +59,41 @@ R__LOAD_LIBRARY(libNPTrackReconstruction.dylib)
x.push_back(X);
DrawCircle(X,Z,R);
}
*/
// add noise
// Generate n points from the line
/* n = (int) rand.Uniform(1,3);
for(unsigned int i = 0 ; i < n ; i++){
double Z = i*stepZ-100;
z.push_back(Z);
double R = rand.Uniform(-maxR,maxR);
r.push_back(R);
double X = rand.Uniform(-100,100);
x.push_back(X);
DrawCircle(X,Z,R);
}
*/
double x0,x100,af,bf;
double res = dcr.BuildTrack2D(x,z,r,x0,x100,af,bf);
cout << res << endl;
DrawTrack(af,bf,1);
// Check 2D track crossing point
// Check 2D Minimize Func
c->cd(2);
bck->Draw();
double ThetaL=-60*deg;
double ThetaH= 60*deg;
TVector3 L(rand.Uniform(-100,100),0,0);
TVector3 H(rand.Uniform(-100,100),0,0);
L.RotateZ(ThetaL);
H.RotateZ(ThetaH);
// direction of U and V wire
TVector3 u(0,1,0);
u.RotateZ(ThetaL);
TVector3 v(0,1,0);
v.RotateZ(ThetaH);
// Compute the coeff of the two line of vecotr u (v) going through H (L)
// dv : y = av*x+bv
double av = v.Y()/v.X();
double bv = H.Y() - av*H.X();
DrawTrack(av,bv,1);
// du : y = au*x+bu
double au = u.Y()/u.X();
double bu = L.Y() - au*L.X();
DrawTrack(au,bu,1);
TVector3 P;
dcr.ResolvePlane(L,ThetaL,H,ThetaH,P);
auto Lmk = new TMarker(L.X(),L.Y(),23);
auto Hmk = new TMarker(H.X(),H.Y(),27);
auto Pmk = new TMarker(P.X(),P.Y(),29);
Lmk->Draw();
Hmk->Draw();
Pmk->Draw();
TGraph* ga = dcr.Scan(af,bf,0,af*0.5,af*2);
TGraph* gb = dcr.Scan(af,bf,1,bf*0.5,bf*2);
ga->SetLineColor(kAzure+7);
ga->Draw("ac");
ga->GetYaxis()->SetRangeUser(0,1000);
gb->Draw("c");
c->Update();
gPad->WaitPrimitive();
}
}
// draw a track based on z=ax+b
void DrawTrack(double a, double b, int style){
double x[2]={-100,100};
......
#include"NPDCReconstruction.h"
#include "Math/Factory.h"
#include "TError.h"
/*****************************************************************************
* Copyright (C) 2009-2020 this file is part of the NPTool Project *
* *
......@@ -19,6 +16,10 @@
* This class have all the method needed to analyse Drift Chambers *
*****************************************************************************/
#include"NPDCReconstruction.h"
#include "Math/Factory.h"
#include "TError.h"
#include "TGraph.h"
using namespace std;
using namespace NPL;
......@@ -30,16 +31,16 @@ DCReconstruction::DCReconstruction(){
m_min->SetPrintLevel(0);
// this avoid error
gErrorIgnoreLevel = kError;
gErrorIgnoreLevel = kError;
//m_min->SetMaxFunctionCalls(1000);
//m_min->SetMaxIterations(1000);
//m_min->SetTolerance(1000);
//m_min->SetTolerance(1);
//m_min->SetPrecision(1e-10);
}
}
////////////////////////////////////////////////////////////////////////////////
DCReconstruction::~DCReconstruction(){
delete m_min;
}
}
////////////////////////////////////////////////////////////////////////////////
double DCReconstruction::BuildTrack2D(const vector<double>& X,const vector<double>& Z,const vector<double>& R,double& X0,double& X100,double& a, double& b ){
......@@ -47,23 +48,22 @@ double DCReconstruction::BuildTrack2D(const vector<double>& X,const vector<doubl
fitZ=&Z;
fitR=&R;
// assume all X,Z,R of same size
size = X.size();
sizeX = X.size();
// Define the starting point of the fit: a straight line passing through the
// the first and last wire
// z = ax+b -> x=(z-b)/a
// ai = (Z[size-1]-Z[0])/(X[size-1]-R[size-1]-X[0]-R[0]);
// bi = Z[0]-ai*(X[0]+R[0]);
ai = (Z[size-1]-Z[0])/(X[size-1]-X[0]);
ai = (Z[sizeX-1]-Z[0])/(X[sizeX-1]-X[0]);
bi = Z[0]-ai*(X[0]);
// parameter[0]=ai;
// parameter[1]=bi;
m_min->Clear();
m_min->SetLimitedVariable(0,"a",ai,1,ai*0.1,ai*10);
m_min->SetLimitedVariable(1,"b",bi,1,bi*0.1,bi*10);
m_min->SetVariable(0,"a",ai,1);
m_min->SetVariable(1,"b",bi,1);
// Perform minimisation
m_min->Minimize();
//std::cout << "EDM:" << m_min->Edm() << std::endl;
// access set of parameter that produce the minimum
const double *xs = m_min->X();
a=xs[0];
......@@ -109,34 +109,54 @@ void DCReconstruction::ResolvePlane(const TVector3& L,const double& ThetaU ,cons
xM=-10000;
yM=-10000;
}
PosXY.SetXYZ(xM,yM,0);
}
////////////////////////////////////////////////////////////////////////////////
double DCReconstruction::SumD(const double* parameter ){
unsigned int sizeX = fitX->size();
// Compute the sum P of the distance between the circle and the track
P = 0;
a = parameter[0];
b = parameter[1];
ab= a*b;
a2=a*a;
for(unsigned int i = 0 ; i < sizeX ; i++){
c = (*fitX)[i];
d = (*fitZ)[i];
r = (*fitR)[i];
x = (a*d-ab+c)/(1+a2);
z = a*x+b;
//P+= sqrt(abs( (x-c)*(x-c)+(z-d)*(z-d)-r*r));
P+= abs( (x-c)*(x-c)+(z-d)*(z-d)-r*r);
p= (x-c)*(x-c)+(z-d)*(z-d)-r*r;
// numerical trick to have a smooth derivative instead of using abs
P+= sqrt(p*p+0.1);
}
// return normalized power
return P/size;
return P/sizeX;
}
////////////////////////////////////////////////////////////////////////////////
TGraph* DCReconstruction::Scan(double a, double b, int tovary, double minV, double maxV){
vector<double> x,y;
unsigned int sizeT=1000;
double step = (maxV-minV)/sizeT;
double p[2]={a,b};
for(unsigned int i = 0 ; i < sizeT ; i++){
if(!tovary){
p[0]=minV+step*i;
x.push_back(p[0]);
y.push_back(SumD(p));
}
else{
p[1]=minV+step*i;
x.push_back(p[1]);
y.push_back(SumD(p));
}
}
TGraph* g = new TGraph(x.size(),&x[0],&y[0]);
return g;
}
......@@ -34,10 +34,11 @@
* intersaction in any given plane. This is done using ResolvePlane. *
* - Resolving plane for two Z plane will provide a point and a direction *
*****************************************************************************/
#include<vector>
#include"TVector3.h"
#include <vector>
#include "TVector3.h"
#include "Math/Minimizer.h"
#include "Math/Functor.h"
class TGraph;
namespace NPL{
class DCReconstruction{
......@@ -57,6 +58,13 @@ namespace NPL{
// Function used by the minimizer in BuildTrack2D
double SumD(const double* parameter );
// For debugging/optimisation
// Scan Sumd versus parameter a or b (tovary =0 for a, 1 for b)
// return a TGraph for display
TGraph* Scan(double a, double b, int tovary, double minV, double maxV);
private: // private member used by SumD
ROOT::Math::Minimizer* m_min;
ROOT::Math::Functor m_func;
......@@ -64,8 +72,8 @@ namespace NPL{
const std::vector<double>* fitZ;
const std::vector<double>* fitR;
// used by SumD
unsigned int size ;
double P,a,b,ab,a2,c,d,x,z,r;
unsigned int sizeX ;
double P,p,a,b,ab,a2,c,d,x,z,r;
// used by BuildTrack
double ai,bi;
double parameter[2];
......@@ -73,7 +81,6 @@ namespace NPL{
long double av,bv,au,bu;
double xM,yM;
};
}
......
Markdown is supported
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