diff --git a/NPLib/Detectors/PISTA/TFPMWPhysics.cxx b/NPLib/Detectors/PISTA/TFPMWPhysics.cxx
index f422030e2c22a7daf854a232d37c12bd9d3f85fc..4da55fb1b7ee62fb1737a23a5672f06f7a45700d 100644
--- a/NPLib/Detectors/PISTA/TFPMWPhysics.cxx
+++ b/NPLib/Detectors/PISTA/TFPMWPhysics.cxx
@@ -46,6 +46,8 @@ TFPMWPhysics::TFPMWPhysics()
: m_EventData(new TFPMWData),
m_PreTreatedData(new TFPMWData),
m_EventPhysics(this),
+ m_Zf(7600),
+ m_GapSize(0),
m_NumberOfDetectors(0){
}
@@ -104,6 +106,8 @@ void TFPMWPhysics::BuildPhysicalEvent() {
QSumX[NX] += QX;
if(MaxQX.find(NX)==MaxQX.end() || MaxQX[NX].second<QX){
MaxQX[NX] = make_pair(StrX,QX);
+ //QXmax[NX] = QX;
+ //StripXmax[NX] = StrX;
}
}
}
@@ -116,13 +120,23 @@ void TFPMWPhysics::BuildPhysicalEvent() {
QSumY[NY] += QY;
if(MaxQY.find(NY)==MaxQY.end() || MaxQY[NY].second<QY){
MaxQY[NY] = make_pair(StrY,QY);
+ //QYmax[NY] = QY;
+ //StripYmax[NY] = StrY;
}
}
}
for(auto &DetN : DetectorHit){
+ //double PosX = WeightedAverage(MapX[DetN]);
+ //double PosY = WeightedAverage(MapY[DetN]);
double PosX = AnalyticHyperbolicSecant(MaxQX[DetN],MapX[DetN]);
double PosY = AnalyticHyperbolicSecant(MaxQY[DetN],MapY[DetN]);
+ //double PosX = FittedHyperbolicSecant(MaxQX[DetN],MapX[DetN]);
+ //double PosY = FittedHyperbolicSecant(MaxQY[DetN],MapY[DetN]);
+
+
+ PosX = -PosX - DetPosX[DetN];
+ PosY = PosY - DetPosY[DetN];
int sx0 = (int) PosX;
int sx1 = sx0+1;
@@ -136,8 +150,7 @@ void TFPMWPhysics::BuildPhysicalEvent() {
PositionY.push_back(PosY);
}
- double Zf = 9000;
- CalculateFocalPlanePosition(Zf);
+ CalculateFocalPlanePosition(m_Zf);
CalculateTargetPosition();
}
@@ -145,14 +158,19 @@ void TFPMWPhysics::BuildPhysicalEvent() {
void TFPMWPhysics::CalculateFocalPlanePosition(double Zf){
if(PositionX.size()==4){
- double Z2 = DetPosZ[2];
- double Z3 = DetPosZ[3];
+ double Z2 = DetPosZ[2]-m_GapSize;
+ double Z3 = DetPosZ[3]-m_GapSize;
double X2 = PositionX[2];
double X3 = PositionX[3];
+ double Y2 = PositionY[2];
+ double Y3 = PositionY[3];
Xf = X2 + (X3-X2)/(Z3-Z2)*(Zf-Z2);
Thetaf = atan((X3-X2)/(Z3-Z2));
+ Z2 = DetPosZ[2]+m_GapSize;
+ Z3 = DetPosZ[3]+m_GapSize;
+ Yf = Y2 + (Y3-Y2)/(Z3-Z2)*(Zf-Z2);
}
}
@@ -161,25 +179,46 @@ void TFPMWPhysics::CalculateFocalPlanePosition(double Zf){
void TFPMWPhysics::CalculateTargetPosition(){
if(PositionX.size()>1){
- double Z0 = DetPosZ[0];
- double Z1 = DetPosZ[1];
+ double Z0 = DetPosZ[0]-m_GapSize;
+ double Z1 = DetPosZ[1]-m_GapSize;
double X0 = PositionX[0];
double X1 = PositionX[1];
double Y0 = PositionY[0];
double Y1 = PositionY[1];
+ Xt = X1 + (X0-X1)*Z1/Z0;
- TVector3 vFF = TVector3(X1-X0,Y1-Y0,Z1-Z0);
+ Z0 = DetPosZ[0]+m_GapSize;
+ Z1 = DetPosZ[1]+m_GapSize;
+ Yt = Y1 + (Y0-Y1)*Z1/Z0;
+ Z0 = DetPosZ[0];
+ Z1 = DetPosZ[1];
+ TVector3 vFF = TVector3(X1-X0,Y1-Y0,Z1-Z0);
Theta_in = vFF.Theta();
Phi_in = vFF.Phi();
+ }
- Xt = X1 + (X0-X1)*Z1/Z0;
- Yt = Y1 + (Y0-Y1)*Z1/Z0;
+}
+/////////////////////////////////////////////////////////////////////
+double TFPMWPhysics::WeightedAverage(std::vector<std::pair<int,double>>& Map){
+
+ unsigned int sizeQ = Map.size();
+
+ double num=0;
+ double denum=0;
+ for(unsigned int i = 0 ; i < sizeQ ; i++){
+ num += Map[i].first * Map[i].second;
+ denum += Map[i].second;
}
+ double weighted_average = num/denum;
+
+ return weighted_average;
+
+
}
-
+
/////////////////////////////////////////////////////////////////////
double TFPMWPhysics::AnalyticHyperbolicSecant(std::pair<int,double>& MaxQ,std::vector<std::pair<int,double>>& Map){
double sech = -1000 ;
@@ -224,6 +263,45 @@ double TFPMWPhysics::AnalyticHyperbolicSecant(std::pair<int,double>& MaxQ,std::v
return sech ;
}
+///////////////////////////////////////////////////////////////////////////
+double TFPMWPhysics::FittedHyperbolicSecant(std::pair<int,double>& MaxQ,std::vector<std::pair<int,double>>& Map){
+ // Warning: should not delete static variable
+ static TF1* f = new TF1("sechs","[0]/(cosh(TMath::Pi()*(x-[1])/[2])*cosh(TMath::Pi()*(x-[1])/[2]))",1,1000);
+
+ // Help the fit by computing the position of the maximum by analytic method
+ double StartingPoint = AnalyticHyperbolicSecant(MaxQ,Map);
+ // if analytic method fails then the starting point in strip max
+ if(StartingPoint==-1000) StartingPoint = MaxQ.first;
+
+ // Maximum is close to charge max, Mean value is close to Analytic one, typical width is 3.8 strip
+ f->SetParameters(MaxQ.first,StartingPoint,1);
+
+ static vector<double> y ;
+ static vector<double> q ;
+ y.clear(); q.clear();
+ double final_size = 0 ;
+ unsigned int sizeQ = Map.size();
+
+ for(unsigned int i = 0 ; i < sizeQ ; i++){
+ if(Map[i].second > (MaxQ.second)*0.2){
+ q.push_back(Map[i].second);
+ y.push_back(Map[i].first);
+ final_size++;
+ }
+ }
+
+ // requiered at least 3 point to perfom a fit
+ if(final_size<3){
+ return -1000 ;
+ }
+
+ TGraph* g = new TGraph(q.size(),&y[0],&q[0]);
+ g->Fit(f,"QN0");
+ delete g;
+ return f->GetParameter(1) ;
+}
+
+
///////////////////////////////////////////////////////////////////////////
void TFPMWPhysics::PreTreat() {
@@ -336,11 +414,19 @@ void TFPMWPhysics::Clear() {
MaxQY.clear();
Xf = -1000;
+ Yf = -1000;
Thetaf = -1000;
Theta_in = -1000;
Phi_in = -1000;
Xt = -1000;
Yt = -1000;
+
+ /*for(int i=0; i<4; i++){
+ QXmax[i] = -1;
+ //QYmax[i] = -1;
+ StripXmax[i] = -1;
+ //StripYmax[i] = -1;
+ }*/
}
diff --git a/NPLib/Detectors/PISTA/TFPMWPhysics.h b/NPLib/Detectors/PISTA/TFPMWPhysics.h
index dc49aa89a6e1aa456d52504b2861cd10e470696d..b55f47a58eb629b03326b0d50445fca1d2f344e6 100644
--- a/NPLib/Detectors/PISTA/TFPMWPhysics.h
+++ b/NPLib/Detectors/PISTA/TFPMWPhysics.h
@@ -32,6 +32,7 @@ using namespace std;
// ROOT headers
#include "TObject.h"
#include "TH1.h"
+#include "TF1.h"
#include "TVector3.h"
#include "TSpline.h"
// NPTool headers
@@ -66,7 +67,12 @@ class TFPMWPhysics : public TObject, public NPL::VDetector {
vector<double> PositionY;
vector<double> ChargeX;
vector<double> ChargeY;
+ //double QXmax[4];
+ //int StripXmax[4];
+ //double QYmax[4];
+ //int StripYmax[4];
double Xf;
+ double Yf;
double Thetaf;
double Xt;
double Yt;
@@ -132,7 +138,12 @@ class TFPMWPhysics : public TObject, public NPL::VDetector {
// clear the raw and physical data objects event by event
void ClearEventPhysics() {Clear();}
- void ClearEventData() {m_EventData->Clear();}
+ void ClearEventData() {m_EventData->Clear();}
+
+ // Get Detector Position
+ double GetDetectorPositionX(int Det) {return DetPosX[Det];}
+ double GetDetectorPositionY(int Det) {return DetPosY[Det];}
+ double GetDetectorPositionZ(int Det) {return DetPosZ[Det];}
//////////////////////////////////////////////////////////////
@@ -154,7 +165,9 @@ class TFPMWPhysics : public TObject, public NPL::VDetector {
double fThresholdX;
double fThresholdY;
+ double WeightedAverage(std::vector<std::pair<int,double>>& Map);
double AnalyticHyperbolicSecant(std::pair<int,double>& MaxQ, std::vector<std::pair<int, double>>& Map);
+ double FittedHyperbolicSecant(std::pair<int,double>& MaxQ, std::vector<std::pair<int, double>>& Map);
// objects are not written in the TTree
@@ -171,6 +184,8 @@ class TFPMWPhysics : public TObject, public NPL::VDetector {
// parameters used in the analysis
private:
double m_E_Threshold; //!
+ double m_Zf; //! Z focal plane position
+ double m_GapSize;
// number of detectors
private:
diff --git a/NPLib/Detectors/PISTA/TICPhysics.cxx b/NPLib/Detectors/PISTA/TICPhysics.cxx
index 34219350aa7e020e94fe566801df939efcd95f33..2a5d69665176af5d43dd01b7a7d2f2bed47270e7 100644
--- a/NPLib/Detectors/PISTA/TICPhysics.cxx
+++ b/NPLib/Detectors/PISTA/TICPhysics.cxx
@@ -86,10 +86,12 @@ void TICPhysics::BuildPhysicalEvent() {
if(fIC[1]>0 && fIC[5]>0){
DE = 0.5*(fIC[0] + fIC[1] + fIC[2] + fIC[3]) + fIC[4];
Eres = fIC[5] + fIC[6] + fIC[7] + fIC[8] + fIC[9];
+ Etot =0.02411*(0.8686*fIC[0]+0.7199*fIC[1]+0.6233*fIC[2]+0.4697*fIC[3]+0.9787*fIC[4]+0.9892*fIC[5]+2.1038*fIC[6]+1.9429*fIC[7]+1.754*fIC[8]+2.5*fIC[9]);
}
else{
DE = -100;
Eres = -100;
+ Etot = -100;
}
}
@@ -175,6 +177,7 @@ void TICPhysics::ReadAnalysisConfig() {
void TICPhysics::Clear() {
DE = -100;
Eres = -100;
+ Etot = -100;
}
diff --git a/NPLib/Detectors/PISTA/TICPhysics.h b/NPLib/Detectors/PISTA/TICPhysics.h
index beac07d38f693281f85f4015b023b6ae33521c2a..255bf831b7ef28db00468da71be66e7d8f26597e 100644
--- a/NPLib/Detectors/PISTA/TICPhysics.h
+++ b/NPLib/Detectors/PISTA/TICPhysics.h
@@ -63,7 +63,7 @@ class TICPhysics : public TObject, public NPL::VDetector {
public:
double DE;
double Eres;
-
+ double Etot;
private:
/// A usefull method to bundle all operation to add a detector