From b3430bfdaf02fe7ee14def8f6892e705ed0e9459 Mon Sep 17 00:00:00 2001
From: Clenain <cyril.lenain1@gmail.com>
Date: Wed, 18 Sep 2019 15:15:44 +0200
Subject: [PATCH] Adding analyse codes for MINOS
---
NPLib/Detectors/Minos/CMakeLists.txt | 28 +-
NPLib/Detectors/Minos/MTreeStructureLinkDef.h | 11 +
NPLib/Detectors/Minos/TMinosClust.cxx | 43 +
NPLib/Detectors/Minos/TMinosClust.h | 26 +
NPLib/Detectors/Minos/TMinosPhysics.cxx | 785 ++++++++++++++----
NPLib/Detectors/Minos/TMinosPhysics.h | 250 ++++--
NPLib/Detectors/Minos/TMinosResult.cxx | 37 +
NPLib/Detectors/Minos/TMinosResult.h | 25 +
NPLib/TrackReconstruction/CMakeLists.txt | 6 +-
.../NPTrackReconstructionLinkDef.h | 1 +
NPLib/TrackReconstruction/Tracking.cxx | 488 +++++++++++
NPLib/TrackReconstruction/Tracking.h | 33 +
12 files changed, 1534 insertions(+), 199 deletions(-)
create mode 100644 NPLib/Detectors/Minos/MTreeStructureLinkDef.h
create mode 100644 NPLib/Detectors/Minos/TMinosClust.cxx
create mode 100644 NPLib/Detectors/Minos/TMinosClust.h
create mode 100644 NPLib/Detectors/Minos/TMinosResult.cxx
create mode 100644 NPLib/Detectors/Minos/TMinosResult.h
create mode 100644 NPLib/TrackReconstruction/Tracking.cxx
create mode 100644 NPLib/TrackReconstruction/Tracking.h
diff --git a/NPLib/Detectors/Minos/CMakeLists.txt b/NPLib/Detectors/Minos/CMakeLists.txt
index c3874b7f6..4e7173cab 100644
--- a/NPLib/Detectors/Minos/CMakeLists.txt
+++ b/NPLib/Detectors/Minos/CMakeLists.txt
@@ -1,6 +1,28 @@
+include(CheckLanguage)
+
add_custom_command(OUTPUT TMinosPhysicsDict.cxx COMMAND ../../scripts/build_dict.sh TMinosPhysics.h TMinosPhysicsDict.cxx TMinosPhysics.rootmap libNPMinos.dylib DEPENDS TMinosPhysics.h)
add_custom_command(OUTPUT TMinosDataDict.cxx COMMAND ../../scripts/build_dict.sh TMinosData.h TMinosDataDict.cxx TMinosData.rootmap libNPMinos.dylib DEPENDS TMinosData.h)
-add_library(NPMinos SHARED TMinosSpectra.cxx TMinosData.cxx TMinosPhysics.cxx TMinosDataDict.cxx TMinosPhysicsDict.cxx )
-target_link_libraries(NPMinos ${ROOT_LIBRARIES} NPCore)
-install(FILES TMinosData.h TMinosPhysics.h TMinosSpectra.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
+add_custom_command(OUTPUT TMinosClustDict.cxx COMMAND ../../scripts/build_dict.sh TMinosClust.h TMinosClustDict.cxx TMinosClust.rootmap libNPMinos.dylib DEPENDS TMinosClust.h)
+add_custom_command(OUTPUT TMinosResultDict.cxx COMMAND ../../scripts/build_dict.sh TMinosResult.h TMinosResultDict.cxx TMinosResult.rootmap libNPMinos.dylib DEPENDS TMinosResult.h)
+
+## Check if MINUIT2 is installed along with ROOT
+find_library(libMinuit_FOUND NAMES Minuit HINTS "$ENV{ROOTSYS}/lib")
+if(libMinuit_FOUND)
+ message(STATUS "Minuit support enabled for Minos.")
+ add_definitions(-DHAVE_MINUIT)
+else()
+ message(STATUS "Minuit support disabled for Minos. Not found in $ENV{ROOTSYS}/lib")
+endif()
+
+add_library(NPMinos SHARED TMinosClust.cxx TMinosData.cxx TMinosPhysics.cxx TMinosResult TMinosDataDict.cxx TMinosPhysicsDict.cxx TMinosClustDict.cxx TMinosResultDict.cxx)
+
+if(libMinuit_FOUND)
+## Link to Minuit2 library
+ target_link_libraries(NPMinos ${ROOT_LIBRARIES} NPCore NPTrackReconstruction Minuit)
+else()
+## No Minuit2 library, skip linking
+ target_link_libraries(NPMinos ${ROOT_LIBRARIES} NPCore NPTrackReconstruction)
+endif()
+
+install(FILES TMinosClust.h TMinosData.h TMinosPhysics.h TMinosResult.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
diff --git a/NPLib/Detectors/Minos/MTreeStructureLinkDef.h b/NPLib/Detectors/Minos/MTreeStructureLinkDef.h
new file mode 100644
index 000000000..8745cb31f
--- /dev/null
+++ b/NPLib/Detectors/Minos/MTreeStructureLinkDef.h
@@ -0,0 +1,11 @@
+#ifdef __CINT__
+
+#pragma link off all globals;
+#pragma link off all classes;
+#pragma link off all functions;
+#pragma link C++ nestedclasses;
+
+#pragma link C++ class TMinosClust+;
+#pragma link C++ class TMinosResult+;
+
+#endif
diff --git a/NPLib/Detectors/Minos/TMinosClust.cxx b/NPLib/Detectors/Minos/TMinosClust.cxx
new file mode 100644
index 000000000..9f2d61d5c
--- /dev/null
+++ b/NPLib/Detectors/Minos/TMinosClust.cxx
@@ -0,0 +1,43 @@
+//---------- Written by C. Santamaria / CEA Saclay !----------
+#include "TMinosClust.h"
+
+ClassImp(TMinosClust)
+
+TMinosClust::TMinosClust()
+: TObject(), x_mm(0), y_mm(0), t_ns(0), z_mm(0), Chargemax(0), n_Cluster(0), n_Pads(0), Chi2(0){
+// Constructor
+}
+TMinosClust::TMinosClust(Double_t xmm, Double_t ymm, Double_t tns, Double_t zmm , Double_t chargemax, Double_t ncluster, Double_t npads, Double_t chi2)
+ : TObject()
+{
+ // Constructor
+ x_mm = xmm;
+ y_mm= ymm;
+ t_ns = tns;
+ z_mm = zmm;
+ Chargemax = chargemax;
+ n_Cluster = ncluster;
+ n_Pads = npads;
+ Chi2 = chi2;
+}
+TMinosClust::~TMinosClust()
+{
+// Destructor
+}
+
+void TMinosClust::Set(Double_t xmm, Double_t ymm, Double_t tns, Double_t zmm, Double_t chargemax, Double_t ncluster, Double_t npads, Double_t chi2)
+{
+ x_mm = xmm;
+ y_mm= ymm;
+ t_ns = tns;
+ z_mm = zmm;
+ Chargemax = chargemax;
+ n_Cluster = ncluster;
+ n_Pads = npads;
+ Chi2 = chi2;
+}
+
+void TMinosClust::SetZ(Double_t zmm)
+{
+ z_mm = zmm;
+}
diff --git a/NPLib/Detectors/Minos/TMinosClust.h b/NPLib/Detectors/Minos/TMinosClust.h
new file mode 100644
index 000000000..d7ee47d67
--- /dev/null
+++ b/NPLib/Detectors/Minos/TMinosClust.h
@@ -0,0 +1,26 @@
+//----------! Written by C. Santamaria / CEA Saclay !----------
+#ifndef TMINOSCLUST_H
+#define TMINOSCLUST_H
+
+#include "TObject.h"
+
+class TMinosClust : public TObject {
+ public:
+ TMinosClust();
+ TMinosClust(Double_t xmm, Double_t ymm, Double_t tns, Double_t zmm , Double_t chargemax, Double_t ncluster, Double_t npads, Double_t chi2);
+ void Set(Double_t xmm, Double_t ymm, Double_t tns, Double_t zmm , Double_t chargemax, Double_t ncluster, Double_t npads, Double_t chi2);
+ void SetZ(Double_t zmm);
+ virtual ~TMinosClust();
+ public:
+ Double_t x_mm;
+ Double_t y_mm;
+ Double_t t_ns;
+ Double_t z_mm;
+ Double_t Chargemax;
+ Double_t n_Cluster;
+ Double_t n_Pads;
+ Double_t Chi2;
+
+ ClassDef(TMinosClust,1);
+};
+#endif // end of #ifndef TMINOSCLUST_H
diff --git a/NPLib/Detectors/Minos/TMinosPhysics.cxx b/NPLib/Detectors/Minos/TMinosPhysics.cxx
index 6ea1a9931..e5302de98 100644
--- a/NPLib/Detectors/Minos/TMinosPhysics.cxx
+++ b/NPLib/Detectors/Minos/TMinosPhysics.cxx
@@ -1,18 +1,19 @@
/*****************************************************************************
- * Copyright (C) 2009-2018 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2018 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: Elidiano Tronchin contact address: tronchin@lpccaen.in2p3.fr *
+ * Original Author: Developed by C. Santamaria / CEA Saclay *
* *
- * Creation Date : October 2018 *
- * Last update : *
+ * Creation Date : 2014/11/24 *
+ * Last update : 2019/09 implemeted in NPTool by Cyril Lenain *
+ * lenain@lpccaen.in2p3.fr *
*---------------------------------------------------------------------------*
* Decription: *
- * This class hold Minos Treated data *
+ * This class hold Minos Treated data *
* *
*---------------------------------------------------------------------------*
* Comment: *
@@ -35,23 +36,40 @@ using namespace std;
#include "RootOutput.h"
#include "NPDetectorFactory.h"
#include "NPOptionManager.h"
+#include "TMinosClust.h"
+#include "TMinosResult.h"
// ROOT
#include "TChain.h"
+#include "TROOT.h"
+#include "TMath.h"
+#include "TMinuit.h"
+#include "/opt/cern/root/root_v5.34.36/math/genvector/inc/Math/Vector3D.h"
+#include "TStyle.h"
+#include "TClonesArray.h"
-ClassImp(TMinosPhysics)
-
+TMinosPhysics* current_phy = 0;
///////////////////////////////////////////////////////////////////////////
TMinosPhysics::TMinosPhysics()
: m_EventData(new TMinosData),
m_PreTreatedData(new TMinosData),
m_EventPhysics(this),
- m_Spectra(0),
+ Tracking_functions(new Tracking),
m_E_RAW_Threshold(0), // adc channels
m_E_Threshold(0), // MeV
+ hfit(new TH1F("hfit","hfit",512,0,512)),
+ grxztmp(new TGraph()),
+ gryztmp(new TGraph()),
+ npoint_temp(0),
+ allevt_2pfiltered(0),
m_NumberOfDetectors(0) {
-}
+ current_phy = this;
+ data_result.SetClass("TMinosResult");
+ fitdata.SetClass("TMinosClust");
+ }
+
+int NclusterFit;
///////////////////////////////////////////////////////////////////////////
/// A usefull method to bundle all operation to add a detector
@@ -63,70 +81,506 @@ void TMinosPhysics::AddDetector(TVector3 , string ){
}
///////////////////////////////////////////////////////////////////////////
-void TMinosPhysics::AddDetector(double R, double Theta, double Phi, string shape){
+void TMinosPhysics::AddDetector(double R, double Theta, TVector3 POS, double Phi, double TargetLenght){
// Compute the TVector3 corresponding
- TVector3 Pos(R*sin(Theta)*cos(Phi),R*sin(Theta)*sin(Phi),R*cos(Theta));
// Call the cartesian method
- AddDetector(Pos,shape);
}
///////////////////////////////////////////////////////////////////////////
void TMinosPhysics::BuildSimplePhysicalEvent() {
- /* BuildPhysicalEvent(); */
+ BuildPhysicalEvent();
}
-
-
///////////////////////////////////////////////////////////////////////////
void TMinosPhysics::BuildPhysicalEvent() {
// apply thresholds and calibration
- /* PreTreat(); */
-
- // match energy and time together
- /* unsigned int mysizeE = m_PreTreatedData->GetMultEnergy(); */
- /* unsigned int mysizeT = m_PreTreatedData->GetMultTime(); */
- /* for (UShort_t e = 0; e < mysizeE ; e++) { */
- /* for (UShort_t t = 0; t < mysizeT ; t++) { */
- /* if (m_PreTreatedData->GetE_DetectorNbr(e) == m_PreTreatedData->GetT_DetectorNbr(t)) { */
- /* DetectorNumber.push_back(m_PreTreatedData->GetE_DetectorNbr(e)); */
- /* Energy.push_back(m_PreTreatedData->Get_Energy(e)); */
- /* Time.push_back(m_PreTreatedData->Get_Time(t)); */
- /* } */
- /* } */
- /* } */
+ PreTreat();
+}
+
+// definition of the fit function for the signal Q(t)
+double TMinosPhysics::conv_fit(double *x, double *p){
+ double val;
+ if(!(x[0]<p[1] || x[0]>512.)) val = p[0] * exp(-3.*(x[0]-p[1])/p[2]) * sin((x[0]-p[1])/p[2]) * pow((x[0]-p[1])/p[2], 3) + 250;
+ //else val = p[3];
+ else val = 250;
+ return(val);
}
///////////////////////////////////////////////////////////////////////////
void TMinosPhysics::PreTreat() {
- /* // This method typically applies thresholds and calibrations */
- /* // Might test for disabled channels for more complex detector */
-
- /* // clear pre-treated object */
- /* ClearPreTreatedData(); */
-
- /* // instantiate CalibrationManager */
- /* static CalibrationManager* Cal = CalibrationManager::getInstance(); */
-
- /* // Energy */
- /* unsigned int mysize = m_EventData->GetMultEnergy(); */
- /* for (UShort_t i = 0; i < mysize ; ++i) { */
- /* if (m_EventData->Get_Energy(i) > m_E_RAW_Threshold) { */
- /* Double_t Energy = Cal->ApplyCalibration("Minos/ENERGY"+NPL::itoa(m_EventData->GetE_DetectorNbr(i)),m_EventData->Get_Energy(i)); */
- /* if (Energy > m_E_Threshold) { */
- /* m_PreTreatedData->SetEnergy(m_EventData->GetE_DetectorNbr(i), Energy); */
- /* } */
- /* } */
- /* } */
- /* // Time */
- /* mysize = m_EventData->GetMultTime(); */
- /* for (UShort_t i = 0; i < mysize; ++i) { */
- /* Double_t Time= Cal->ApplyCalibration("Minos/TIME"+NPL::itoa(m_EventData->GetT_DetectorNbr(i)),m_EventData->Get_Time(i)); */
- /* m_PreTreatedData->SetTime(m_EventData->GetT_DetectorNbr(i), Time); */
- /* } */
-}
+ // This method applies thresholds and calibrations
+ // Isolate 2D tracks with the modified HoughTransformation
+ // Calculate z for each pad
+
+ data_result.Clear();
+ fitdata.Clear();
+
+ TMinosClust *minosfitdata;
+
+ ////////////////////fit variables//////////////////////
+
+ fit_function = new TF1("fit_function",conv_fit, 0, 511, 3);
+ int fit2DStatus = 0.;
+ double fit_function_max = 0., fit_function_Tpad = 0.;
+ double Chi2 = 0.;
+ double Tshaping = 333.9;
+ const double TimeBinElec = 30.; // in ns
+ double MINOSthresh, VDrift, Z_Shift, DelayTrig;
+
+ if(SimulationBool){ // case of simulated data
+ MINOSthresh = 1000000000.; //
+ DelayTrig = 0.;
+ /* VDrift = 0.036; */
+ VDrift = 0.03;
+ Z_Shift = 0;
+ }
+ else{ // case of experiment data
+ MINOSthresh = 100000.;
+ DelayTrig = 10170;
+ VDrift = 0.03475; // in mm/ns
+ /* Z_Shift = 57.16; // */
+ Z_Shift = 0. ; //
+ }
+
+ ///////////////////////////////////////////////////////
+ ClearPreTreatedData();
+ /// Read Q(t) signal to fit and extract DriftTime
+
+ unsigned int NbOfPad = m_EventData->GetPadNumberMult();
+ filled =0;trackNbr=0;
+ trackNbr_FINAL=0;
+ ringsum=0;
+ zmax=0.;Iteration=0;
+ for (unsigned int e = 0; e < NbOfPad; e++){
+
+ x_mm = 0.; y_mm = 0.; z_mm = 0.;maxCharge = 0.; t_pad =0; q_pad=0;
+ x_mm = m_EventData->GetPadX(e);
+ y_mm = m_EventData->GetPadY(e);
+
+ if ( !(abs(x_mm)<0.0001 && abs(y_mm)<0.0001) ) {
+ for (UShort_t o = 0; o < m_EventData->GetTime(e).size() ; ++o){
+ if(m_EventData->GetCharge(e)[o]> maxCharge){
+ maxCharge = m_EventData->GetCharge(e)[o];
+ }
+ }
+
+ if((maxCharge < MINOSthresh && round((sqrt(x_mm*x_mm+y_mm*y_mm)-44.15)/2.1) < 16 && round((sqrt(x_mm*x_mm+y_mm*y_mm)-44.15)/2.1) >4 )) { // to remove rings 1,2,3,4 and 16,17,18
+ /* if((maxCharge < MINOSthresh && ( SimulationBool==1 || round((sqrt(x_mm*x_mm+y_mm*y_mm)-44.15)/2.1) < 16 && round((sqrt(x_mm*x_mm+y_mm*y_mm)-43.8)/2.1) >4 )) ) { // to remove rings 1,2,3,4 and 16,17,18 // to not apply the ring condition for simulation*/
+ Xpad.push_back(x_mm);
+ Ypad.push_back(y_mm);
+ Qpad.push_back(maxCharge);
+ filled++; //
+ }
+ }
+ } //end of NbOfPad
+
+ if(filled>0) { // Nbr of activated pads
+ while(Xpad.size()>=10 && Iteration< 20 ) { // !!!!!!!!
+ filter_result = 0; // Nbr of pads in the Track
+ Iteration++;
+ XpadTemp.clear();
+ YpadTemp.clear();
+ QpadTemp.clear();
+ clusterringboolTemp.clear();
+
+ filter_result = Tracking_functions->Hough_modified(&Xpad, &Ypad, &Qpad, &XpadTemp, &YpadTemp, &QpadTemp, &clusterringboolTemp);
+
+ if(filter_result<0) break;
+ if(filter_result>10 && clusterringboolTemp.back()==1){
+ /* if(filter_result>4 ){ // modified for the simu w 0 dispersion */
+ trackNbr++;
+ for(int ik=0; ik<filter_result; ik++) {
+ XpadNew.push_back(XpadTemp[ik]);
+ YpadNew.push_back(YpadTemp[ik]);
+ ZpadNew.push_back(-10000);
+ QpadNew.push_back(QpadTemp[ik]);
+ clusterringbool.push_back(clusterringboolTemp[ik]);
+ clusternbr.push_back(Iteration);
+ clusterpads.push_back(filter_result);
+ }
+ }
+
+
+ }
+ }
+
+ for(unsigned int il=0; il<XpadNew.size(); il++){
+ minosfitdata = (TMinosClust*)fitdata.ConstructedAt(il);
+ minosfitdata->Set(XpadNew[il], YpadNew[il],
+ -10000, -10000, QpadNew[il],
+ clusternbr[il],clusterpads[il], 0.);
+ ZpadNew.push_back(-10000.);
+
+ }
+
+ if(trackNbr>0 && trackNbr < 5 ){ ///////////
+
+ if(filled==0) cerr << "Error !!!" << endl;
+ //-------------------------------------------------------
+ // STEP 4.1: Fitting taken pads for Qmax and Ttrig info
+ //-------------------------------------------------------
+
+ for(unsigned int i=0 ; i< NbOfPad; i++ ) {
+ hfit->Reset();
+ bool fitbool = false;
+ x_mm = m_EventData->GetPadX(i);
+ y_mm = m_EventData->GetPadY(i);
+
+ for(unsigned int jj=0; jj<XpadNew.size(); jj++) {
+ if( abs(XpadNew[jj]-x_mm)<0.0001 && abs(YpadNew[jj]-y_mm)<0.0001) {
+ fitbool = true;
+ indexfill=jj;
+ break;
+ }
+ }
+
+ if( fitbool==true ) {
+
+ /* TGraph gfit( m_EventData->GetTime(i).size() ,&(m_EventData->GetTime(i)[0]),&(m_EventData->GetCharge(i)[0])); */
+ /* for (UShort_t o = 0; o < m_EventData->GetTime(i).size() ; ++o){ */
+ /* gfit.SetPoint(o,m_EventData->GetTime(i)[o],m_EventData->GetCharge(i)[o]+250); */
+ /* if(m_EventData->GetCharge(i)[o]+250> hfit_max){ */
+ /* hfit_max_T = m_EventData->GetTime(i)[o]; */
+ /* hfit_max=m_EventData->GetCharge(i)[o]+250; */
+ /* } */
+ /* } */
+
+ for(Int_t j=0; j< m_EventData->GetTime(i).size(); j++) {
+ if(m_EventData->GetCharge(i)[j]>=0){
+ hfit->SetBinContent(hfit->FindBin(m_EventData->GetTime(i)[j]), m_EventData->GetCharge(i)[j]+250);
+ }
+ }
+
+ // Fitting the hfit histogram of last channel if not empty
+
+ if(hfit->GetSumOfWeights()>0) {
+ hfit->GetXaxis()->SetRange(0,510);
+ hfit_max = hfit->GetMaximum();
+ hfit_max_T = hfit->GetMaximumBin();
+ T_min=-1;
+ T_max=-1;
+
+ // Find the T_min & T_max limits of the signal non zero
+ for(int h=hfit_max_T;h>0;h--) {
+ if(T_min == -1 && (hfit->GetBinContent(h))<=250 ) {
+ T_min = h;
+ break;
+ }
+ }
+ for(int h=hfit_max_T;h<510;h++) {
+ if(T_max == -1 && (hfit->GetBinContent(h)) < 1 ) {
+ T_max = h;
+ break;
+ }
+ }
+
+ if((hfit_max_T-3.5*(Tshaping/TimeBinElec)) > T_min) T_min = hfit_max_T-2*Tshaping/TimeBinElec;
+ if((hfit_max_T+10) < T_max || T_max==-1) T_max = hfit_max_T+10.;
+ T_min = max(T_min,0.);
+ if(T_max>510) T_max = 510;
+
+ // Set fit parameters
+ fit_function->SetParameter(0, hfit_max-250.);
+ fit_function->SetParameter(1,hfit_max_T - Tshaping/TimeBinElec);
+ fit_function->SetParameter(2, Tshaping/TimeBinElec);
+ fit_function->SetParLimits(0,0,MINOSthresh);
+ fit_function->SetParLimits(1,-20,512);
+ fit_function->SetParLimits(2,0,512);
+
+ fit2DStatus = hfit->Fit(fit_function,"QN","",T_min,T_max);
+ /* fit2DStatus = gfit.Fit(fit_function,"QN","",hfit_max-2*Tshaping,hfit_max+2*Tshaping); */
+
+ double fit_function_max = 0., fit_function_Tpad = 0.;
+ if(fit2DStatus==0) {
+ Chi2 = fit_function->GetChisquare();
+ fit_function_max = fit_function->GetMaximum();
+ fit_function_Tpad = fit_function->GetParameter(1);
+ }
+
+ //attribute q_pad and z_mm value
+ if(fit2DStatus!=0 || fit_function_max<=20. ||
+ fit_function_max >= MINOSthresh || fit_function_Tpad<=0.15 ||
+ fit_function_Tpad>=513. || fit_function->GetParameter(2)<=0.15 ||
+ fit_function->GetParameter(2)>=513.) {
+ q_pad = hfit_max-250.;
+ z_mm = -10000;
+
+ }
+
+ else {
+ t_pad = fit_function_Tpad;
+ /* z_mm = 300.-(-t_pad*TimeBinElec+DelayTrig)*VDrift[int((sqrt(x_mm*x_mm+y_mm*y_mm)-43.8)/2.1)] */
+
+ if(SimulationBool)z_mm = t_pad*TimeBinElec*VDrift; // for simu
+ else z_mm = 300.-(-t_pad*TimeBinElec+DelayTrig)*VDrift;
+
+ Q_Pad.push_back(fit_function_max-250);
+
+ T_Pad.push_back(t_pad*TimeBinElec);
+ X_Pad.push_back(x_mm);
+ Y_Pad.push_back(y_mm);
+ Z_Pad.push_back(z_mm);
+ q_pad = fit_function_max-250.;
+
+ }
+ ZpadNew[indexfill] = z_mm;
+ minosfitdata = (TMinosClust*)fitdata.ConstructedAt(indexfill);
+ minosfitdata->Set(XpadNew[indexfill], YpadNew[indexfill], t_pad*TimeBinElec, ZpadNew[indexfill], q_pad, clusternbr[indexfill], clusterpads[indexfill], Chi2);
+ QpadNew[indexfill] = q_pad ;
+
+ } // end if SumOfWeights > 0
+
+ } // end if fitbool==True
+ else continue;
+ } // end of NbOfPad
+
+ //-------------------------------------------------------
+ // STEP 3.2: Filtering the tracks off possible noise
+ // with Hough3D (3*2D planes)
+ //-------------------------------------------------------
+
+ cluster_temp = 0;
+ int ringtouch[19]={0};
+ /* for(int ko=1; ko<19; ko++) ringtouch[ko] = 0; ///// Added by Cyril */
+
+ for(unsigned int i=0;i<(XpadNew.size());i++){
+ if(xin.size()>0 && ((cluster_temp!=int(clusternbr[i]) && i!=0) || i==(XpadNew.size() - 1))){
+
+ Tracking_functions->Hough_3D(&xin, &yin, &zin, &qin, &xout, &yout, &zout, &qout);
+
+ for(unsigned int ij=0; ij<xout.size();ij++){
+ if(zout[ij]>zmax) {zmax = zout[ij];}
+ ringtouch[int(round((sqrt(xout[ij]*xout[ij]+yout[ij]*yout[ij])-44.15)/2.1))]++; // Corr by Cyril
+ }
+
+ for(int ko=0; ko<19; ko++){
+ if(ringtouch[ko]>0) ringsum++;
+ }
+
+ /* if(zmax> 290 ) ringsum = 16; //comment it if z_mm is not well reconstruct by Vdrift and DelayTrig ! */
+
+ // Tracks of interest: >10 pads and >=12 rings hit
+
+ /* if(xout.size()>10 &&(SimulationBool==1 || ringsum>=8)){ */
+ if(xout.size()>10 && ringsum>=8){
+
+ npoint=0;
+ trackNbr_FINAL++;
+ double charge_temp=0.;
+ double lenght_temp=0;
+ double zintarget=300;
+ double zouttarget=0;
+ int indexin=0; int indexout=0;
+ for(unsigned int ij=0; ij<xout.size(); ij++){
+
+ point.SetXYZ(xout[ij],yout[ij],zout[ij]);
+ point.RotateZ(0.);//15.6*TMath::DegToRad() CHANGE####
+ xoutprime.push_back(point.X());youtprime.push_back(point.Y());zoutprime.push_back(point.Z());
+ TOTxoutprime.push_back(point.X());
+ // Added line to see events used for fit of lines in tree
+ TOTyoutprime.push_back(point.Y());
+ TOTzoutprime.push_back(point.Z());
+ TOTqout.push_back(qout[ij]);
+ // save cluster_temp
+ // save numbtrack
+ trackclusternbr.push_back(clusternbr[i]);
+ tracknbr.push_back( trackNbr_FINAL ) ;
+
+ grxztmp->SetPoint(npoint,zoutprime[ij],xoutprime[ij]);
+ gryztmp->SetPoint(npoint,zoutprime[ij],youtprime[ij]);
+ charge_temp += qout[ij];
+ minosdata_result = (TMinosResult*)data_result.ConstructedAt(array_final);
+ minosdata_result->Set(xoutprime[ij], youtprime[ij], zoutprime[ij], qout[ij], trackNbr_FINAL, xout.size(), zmax);
+ array_final++;
+ npoint++;
+
+ ///////// to exclude data w/ rings from 15to18
+ if(sqrt(xout[ij]*xout[ij]+yout[ij]*yout[ij])<70 && sqrt(xout[ij]*xout[ij]+yout[ij]*yout[ij])> 53 ){
+ if(zoutprime[ij]<zintarget) {zintarget=zoutprime[ij];indexin=ij;}
+ if(zoutprime[ij]>zouttarget) {zouttarget=zoutprime[ij];indexout=ij;}
+ }
+ }
+
+ if(xout.size()>0)lenght_temp=sqrt(pow(zoutprime[indexin]-zoutprime[indexout],2)+pow(youtprime[indexin]-youtprime[indexout],2)+pow(xoutprime[indexin]-xoutprime[indexout],2));
+
+ grxz.push_back(*grxztmp);
+ gryz.push_back(*gryztmp);
+ grxztmp->Set(0);
+ gryztmp->Set(0);
+
+ chargeTot.push_back(charge_temp);
+ lenght.push_back(lenght_temp);
+ } // end of if(xout.size()>10 &&)
+
+ //X_mm->Fill();
+
+ xin.clear();
+ yin.clear();
+ zin.clear();
+ qin.clear();
+ xout.clear();
+ yout.clear();
+ zout.clear();
+ xoutprime.clear();
+ youtprime.clear();
+ zoutprime.clear();
+ qout.clear();
+ npoint_temp=0;
+ ringsum=0;// angle between 1st track and z axis in 3D in degrees
+ zmax=0.;
+ for(int ko=0; ko<18; ko++) ringtouch[ko] = 0;
+
+ } // if(xin.size()>0 && (( cluster_temp .... )
+
+ cluster_temp = clusternbr[i];
+
+ if(!(clusterpads[i]>=10 && clusterringbool[i]==1 && ZpadNew[i]>-100 && ZpadNew[i]<=310)) continue;
+ else
+ {
+ xin.push_back(XpadNew[i]);
+ yin.push_back(YpadNew[i]);
+ zin.push_back(ZpadNew[i]);
+ qin.push_back(QpadNew[i]);
+ npoint_temp++;
+ }
+
+ }//end of loop on padsNews
+
+
+ /* //------------------------------------------------------- */
+ /* // STEP 3.2: Fitting the filtered tracks in 3D */
+ /* // (weight by charge, TMinuit) */
+ /* //------------------------------------------------------- */
+
+
+ if(trackNbr_FINAL>= 1){ // !!! Just for 1 or more tracks!
+ /* if(trackNbr_FINAL== 2){ */
+
+ //////////Minimization in 3D to reconstruct track lines
+ allevt_2pfiltered++;
+
+ for(int itr= 0 ; itr<trackNbr_FINAL; itr++) {
+
+
+ pStart[0]=0; pStart[2]=0; pStart[1]=1; pStart[3]=3;
+ min = new TMinuit(4);
+ min->SetPrintLevel(-1);
+ arglist[0] = 3;
+
+ Tracking_functions->FindStart(pStart,chi,fitStatus, &grxz.at(itr), &gryz.at(itr));
+
+ NclusterFit = itr+1;
+ current_phy=this;
+ min->SetFCN(SumDistance);
+
+
+ // Set starting values and step sizes for parameters
+ min->mnparm(0,"x0",pStart[0],0.1,-500,500,iflag);
+ min->mnparm(1,"Ax",pStart[1],0.1,-10,10,iflag);
+ min->mnparm(2,"y0",pStart[2],0.1,-500,500,iflag);
+ min->mnparm(3,"Ay",pStart[3],0.1,-10,10,iflag);
+ arglist[0] = 200; // number of function calls
+ arglist[1] = 0.000001; // tolerance
+
+ min->mnexcm("MIGRAD",arglist,2,iflag); // minimization with MIGRAD
+
+ min->mnstat(amin,edm,errdef,nvpar,nparx,iflag); //returns current status of the minimization
+
+
+ // get fit parameters
+
+ for(int i = 0; i <4; i++) min->GetParameter(i,parFit_temp[i],err_temp[i]);
+
+ /* if( (parFit_temp[0] >-499 && parFit_temp[0]<499) && (parFit_temp[2] >-499 && parFit_temp[2]<499)) { */
+ parFit1.push_back(parFit_temp[0]);
+ parFit2.push_back(parFit_temp[1]);
+ parFit3.push_back(parFit_temp[2]);
+ parFit4.push_back(parFit_temp[3]);
+ /* } */
+ delete min;
+ }
+
+ if(trackNbr_FINAL==2){
+
+ double ParTrack1[4];
+ double ParTrack2[4];
+
+ ParTrack1[0] = parFit1[0];
+ ParTrack1[1] = parFit2[0];
+ ParTrack1[2] = parFit3[0];
+ ParTrack1[3] = parFit4[0];
+
+ ParTrack2[0] = parFit1[1];
+ ParTrack2[1] = parFit2[1];
+ ParTrack2[2] = parFit3[1];
+ ParTrack2[3] = parFit4[1];
+
+ Dist_min=0, Theta_tr1=0, Theta_tr2=0;
+ Tracking_functions->vertex(ParTrack1, ParTrack2, xv, yv, zv, Dist_min, Theta_tr1, Theta_tr2);
+
+ /* if(abs(xv)<20 && abs(yv)<20){ */
+ Xvertex.push_back(xv);
+ Yvertex.push_back(yv);
+ Zvertex.push_back(zv-Z_Shift);
+
+ /* if (zv <-50){ */
+ /* cout << ParTrack1[0]<< " " << ParTrack1[1]<< " " << ParTrack1[2]<<" " << ParTrack1[3]<<endl; */
+ /* cout << xv<< " " << yv<< " " << zv<<endl; */
+ /* } */
+
+
+ if(Theta_tr1>Theta_tr2){
+ Theta1.push_back(Theta_tr1);
+ Theta2.push_back(Theta_tr2);
+ }
+ else if(Theta_tr1<Theta_tr2){
+ Theta1.push_back(Theta_tr2);
+ Theta2.push_back(Theta_tr1);
+ }
+
+ sumTheta.push_back(Theta_tr1+Theta_tr2);
+ /* sumTheta.push_back(156.9-(Theta_tr1+Theta_tr2)); */
+ /* } */
+ /* else{ */
+
+ /* Xvertex.push_back(-99); */
+ /* Yvertex.push_back(-99); */
+ /* Zvertex.push_back(-99); */
+ /* sumTheta.push_back(1000); */
+ /* } */
+
+ }
+ else
+ {
+ Xvertex.push_back(-1000);
+ Yvertex.push_back(-1000);
+ Zvertex.push_back(-1000);
+ sumTheta.push_back(1000);
+ }
+
+ }// end if trackNbr_FINAL>=1
+ else
+ {
+ Xvertex.push_back(-1000);
+ Yvertex.push_back(-1000);
+ Zvertex.push_back(-1000);
+ sumTheta.push_back(1000);
+ }
+
+ } // end loop if 0<trackNbr < 5
+ else
+ {
+ Xvertex.push_back(-1000);
+ Yvertex.push_back(-1000);
+ Zvertex.push_back(-1000);
+ sumTheta.push_back(1000);
+ }
+ // instantiate CalibrationManager
+ static CalibrationManager* Cal = CalibrationManager::getInstance();
+}
///////////////////////////////////////////////////////////////////////////
void TMinosPhysics::ReadAnalysisConfig() {
@@ -188,100 +642,158 @@ void TMinosPhysics::ReadAnalysisConfig() {
}
}
}
-}
-
-
-///////////////////////////////////////////////////////////////////////////
-void TMinosPhysics::Clear() {
- DetectorNumber.clear();
- Energy.clear();
- Time.clear();
}
+double TMinosPhysics::distance2(double x,double y,double z, double *p) {
+ // distance line point is D= | (xp-x0) cross ux |
+ // where ux is direction of line and x0 is a point in the line (like t = 0)
+ ROOT::Math::XYZVector xp(x,y,z); //point of the track
+ ROOT::Math:: XYZVector x0(p[0], p[2], 0. );
+ ROOT::Math::XYZVector x1(p[0] + p[1], p[2] + p[3], 1. ); //line
+ ROOT::Math::XYZVector u = (x1-x0).Unit();
+ double d2 = ((xp-x0).Cross(u)) .Mag2();
+ return d2;
+}
-
-///////////////////////////////////////////////////////////////////////////
-void TMinosPhysics::ReadConfiguration(NPL::InputParser parser) {
- vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Minos");
- if(NPOptionManager::getInstance()->GetVerboseLevel())
- cout << "//// " << blocks.size() << " detectors found " << endl;
-
- vector<string> cart = {"POS","Shape"};
- vector<string> sphe = {"R","Theta","Phi","Shape"};
-
- for(unsigned int i = 0 ; i < blocks.size() ; i++){
- if(blocks[i]->HasTokenList(cart)){
- if(NPOptionManager::getInstance()->GetVerboseLevel())
- cout << endl << "//// Minos " << 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(NPOptionManager::getInstance()->GetVerboseLevel())
- cout << endl << "//// Minos " << 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;
- exit(1);
+void TMinosPhysics::SumDistance(int &, double *, double & sum, double * par, int) {
+ TMinosPhysics* phy = current_phy;
+ int nused=0;
+ double qtot=0;
+ sum = 0;
+
+ for(int i=0; i < phy->data_result.GetEntriesFast(); i++)
+ {
+ phy->minosdata_result = (TMinosResult*)phy->data_result.At(i);
+ if(phy->minosdata_result->n_Cluster==NclusterFit)
+ {
+ float x=phy->minosdata_result->x_mm;
+ float y=phy->minosdata_result->y_mm;
+ float z=phy->minosdata_result->z_mm;
+ float q=phy->minosdata_result->Chargemax;
+ //if(nused<2)cout<<minosdata_result->n_Cluster<<" "<<x<<" "<<y<<" "<<z<<" "<<q<<endl;
+ /* double d = Tracking_functions->distance2(x, y, z, par); */
+ double d = TMinosPhysics::distance2(x, y, z, par);
+ sum += d*q;
+ nused++;
+ qtot+=q;
+ }
}
- }
+ //sum/=nused;
+ sum/=qtot;
+ return;
}
///////////////////////////////////////////////////////////////////////////
-void TMinosPhysics::InitSpectra() {
- m_Spectra = new TMinosSpectra(m_NumberOfDetectors);
-}
-
-
-
-///////////////////////////////////////////////////////////////////////////
-void TMinosPhysics::FillSpectra() {
- m_Spectra -> FillRawSpectra(m_EventData);
- m_Spectra -> FillPreTreatedSpectra(m_PreTreatedData);
- m_Spectra -> FillPhysicsSpectra(m_EventPhysics);
-}
-
+void TMinosPhysics::Clear() {
+
+ Xpad.clear();
+ Ypad.clear();
+ Qpad.clear();
+ XpadNew.clear();
+ YpadNew.clear();
+ ZpadNew.clear();
+ QpadNew.clear();
+
+ clusterringboolTemp.clear();
+ clusterringbool.clear();
+ clusternbr.clear();
+ clusterpads.clear();
+ hfit->Reset();
+
+ xin.clear();
+ yin.clear();
+ zin.clear();
+ qin.clear();
+ xout.clear();
+ yout.clear();
+ zout.clear();
+ xoutprime.clear();
+ youtprime.clear();
+ zoutprime.clear();
+ qout.clear();
+
+ trackclusternbr.clear();
+ tracknbr.clear();
+ TOTxoutprime.clear();
+ TOTyoutprime.clear();
+ TOTzoutprime.clear();
+ TOTqout.clear();
+
+ lenght.clear();
+ chargeTot.clear();
+ parFit1.clear();
+ parFit2.clear();
+ parFit3.clear();
+ parFit4.clear();
+
+ Xvertex.clear();
+ Yvertex.clear();
+ Zvertex.clear();
+ sumTheta.clear();
+
+ grxz.clear();
+ gryz.clear();
+
+ Theta1.clear();
+ Theta2.clear();
+
+ hfit->Reset();
+
+
+ filled=0;
+ indexfill = 0;
+ ChargeBin = 0.;
+ maxCharge = 0.;
+ Iteration=0;
+ filter_result=0;
+ fit2DStatus=0;
+ trackNbr=0;
+ trackNbr_FINAL=0;
+ x_mm = 0.; y_mm = 0.; z_mm = 0.; q_pad = 0.; t_pad = 0.;
+ array_final=0;
+ ringsum=0;
+ zmax=0.;
-///////////////////////////////////////////////////////////////////////////
-void TMinosPhysics::CheckSpectra() {
- m_Spectra->CheckSpectra();
}
-
-
///////////////////////////////////////////////////////////////////////////
-void TMinosPhysics::ClearSpectra() {
- // To be done
-}
-
+void TMinosPhysics::ReadConfiguration(NPL::InputParser parser) {
+ /* vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Minos"); */
+ /* if(NPOptionManager::getInstance()->GetVerboseLevel()) */
+ /* cout << "//// " << blocks.size() << " detectors found " << endl; */
+ /* vector<string> cart = {"POS","Shape"}; */
+ /* vector<string> sphe = {"R","THETA","POS","Phi","TargetLength"}; */
-///////////////////////////////////////////////////////////////////////////
-map< string , TH1*> TMinosPhysics::GetSpectra() {
- if(m_Spectra)
- return m_Spectra->GetMapHisto();
- else{
- map< string , TH1*> empty;
- return empty;
- }
-}
-
-///////////////////////////////////////////////////////////////////////////
-void TMinosPhysics::WriteSpectra() {
- m_Spectra->WriteSpectra();
+ /* for(unsigned int i = 0 ; i < blocks.size() ; i++){ */
+ /* if(blocks[i]->HasTokenList(cart)){ */
+ /* if(NPOptionManager::getInstance()->GetVerboseLevel()) */
+ /* cout << endl << "//// Minos " << 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(NPOptionManager::getInstance()->GetVerboseLevel()) */
+ /* cout << endl << "//// Minos " << i+1 << endl; */
+ /* double R = blocks[i]->GetDouble("R","mm"); */
+ /* double Theta = blocks[i]->GetDouble("THETA","deg"); */
+ /* TVector3 POS = blocks[i]->GetTVector3("POS","cm"); */
+ /* double Phi = blocks[i]->GetDouble("PHI","deg"); */
+ /* double TargetLenght = blocks[i]->GetDouble("TargetLength","mm"); */
+ /* AddDetector(R,Theta,POS,Phi,TargetLenght); */
+ /* } */
+ /* else{ */
+ /* cout << "ERROR: check your input file formatting " << endl; */
+ /* exit(1); */
+ /* } */
+ /* } */
}
-
///////////////////////////////////////////////////////////////////////////
void TMinosPhysics::AddParameterToCalibrationManager() {
CalibrationManager* Cal = CalibrationManager::getInstance();
@@ -291,33 +803,26 @@ void TMinosPhysics::AddParameterToCalibrationManager() {
}
}
-
-
///////////////////////////////////////////////////////////////////////////
void TMinosPhysics::InitializeRootInputRaw() {
TChain* inputChain = RootInput::getInstance()->GetChain();
inputChain->SetBranchStatus("Minos", true );
inputChain->SetBranchAddress("Minos", &m_EventData );
+ if( strlen(inputChain->GetTree()->GetName())==13) SimulationBool= true; // Test if TreeName==SimulatedTree
}
-
-
///////////////////////////////////////////////////////////////////////////
void TMinosPhysics::InitializeRootInputPhysics() {
TChain* inputChain = RootInput::getInstance()->GetChain();
inputChain->SetBranchAddress("Minos", &m_EventPhysics);
}
-
-
///////////////////////////////////////////////////////////////////////////
void TMinosPhysics::InitializeRootOutput() {
TTree* outputTree = RootOutput::getInstance()->GetTree();
outputTree->Branch("Minos", "TMinosPhysics", &m_EventPhysics);
}
-
-
////////////////////////////////////////////////////////////////////////////////
// Construct Method to be pass to the DetectorFactory //
////////////////////////////////////////////////////////////////////////////////
@@ -325,8 +830,6 @@ NPL::VDetector* TMinosPhysics::Construct() {
return (NPL::VDetector*) new TMinosPhysics();
}
-
-
////////////////////////////////////////////////////////////////////////////////
// Registering the construct method to the factory //
////////////////////////////////////////////////////////////////////////////////
diff --git a/NPLib/Detectors/Minos/TMinosPhysics.h b/NPLib/Detectors/Minos/TMinosPhysics.h
index f35fa2777..23d9d57e5 100644
--- a/NPLib/Detectors/Minos/TMinosPhysics.h
+++ b/NPLib/Detectors/Minos/TMinosPhysics.h
@@ -1,20 +1,20 @@
#ifndef TMinosPHYSICS_H
#define TMinosPHYSICS_H
/*****************************************************************************
- * Copyright (C) 2009-2018 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2018 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: Elidiano Tronchin contact address: tronchin@lpccaen.in2p3.fr *
+ * Original Author: Cyril Lenain contact address: lenain@lpccaen.in2p3.fr *
* *
- * Creation Date : October 2018 *
+ * Creation Date : 2019 *
* Last update : *
*---------------------------------------------------------------------------*
* Decription: *
- * This class hold Minos Treated data *
+ * This class hold Minos Treated data *
* *
*---------------------------------------------------------------------------*
* Comment: *
@@ -32,16 +32,23 @@ using namespace std;
#include "TObject.h"
#include "TH1.h"
#include "TVector3.h"
+#include "TF1.h"
+#include "TH2F.h"
+#include "TMinuit.h"
+
+
// NPTool headers
#include "TMinosData.h"
-#include "TMinosSpectra.h"
#include "NPCalibrationManager.h"
#include "NPVDetector.h"
#include "NPInputParser.h"
-// forward declaration
-class TMinosSpectra;
+#include "Tracking.h"
+#include "TMinosResult.h"
+#include "TClonesArray.h"
+using namespace NPL;
+// forward declaration
class TMinosPhysics : public TObject, public NPL::VDetector {
//////////////////////////////////////////////////////////////
@@ -49,91 +56,86 @@ class TMinosPhysics : public TObject, public NPL::VDetector {
public:
TMinosPhysics();
~TMinosPhysics() {};
-
-
//////////////////////////////////////////////////////////////
// Inherited from TObject and overriden to avoid warnings
public:
void Clear();
- void Clear(const Option_t*) {};
-
- //////////////////////////////////////////////////////////////
+ void Clear(const Option_t*) {};
+
+ /* static void SumDistance(int &, double *, double & sum, double * par, int); */
+
+ //////////////////////////////////////////////////////////////
// data obtained after BuildPhysicalEvent() and stored in
// output ROOT file
public:
- vector<int> DetectorNumber;
- vector<double> Energy;
- vector<double> Time;
+
+ vector<double> X_Pad;
+ vector<double> Y_Pad;
+ vector<double> Z_Pad;
+ vector<double> Q_Pad;
+ vector<double> T_Pad;
/// A usefull method to bundle all operation to add a detector
- void AddDetector(TVector3 POS, string shape);
- void AddDetector(double R, double Theta, double Phi, string shape);
+ void AddDetector(TVector3 POS, string shape);
+ void AddDetector(double R, double Theta,TVector3 POS, double Phi, double TargetLenght);
//////////////////////////////////////////////////////////////
// methods inherited from the VDetector ABC class
public:
// read stream from ConfigFile to pick-up detector parameters
void ReadConfiguration(NPL::InputParser);
-
// add parameters to the CalibrationManger
void AddParameterToCalibrationManager();
-
// method called event by event, aiming at extracting the
// physical information from detector
void BuildPhysicalEvent();
-
// same as BuildPhysicalEvent() method but with a simpler
// treatment
void BuildSimplePhysicalEvent();
-
// same as above but for online analysis
void BuildOnlinePhysicalEvent() {BuildPhysicalEvent();};
-
// activate raw data object and branches from input TChain
// in this method mother branches (Detector) AND daughter leaves
// (fDetector_parameter) have to be activated
void InitializeRootInputRaw();
-
// activate physics data object and branches from input TChain
// in this method mother branches (Detector) AND daughter leaves
// (fDetector_parameter) have to be activated
void InitializeRootInputPhysics();
-
// create branches of output ROOT file
void InitializeRootOutput();
-
// clear the raw and physical data objects event by event
void ClearEventPhysics() {Clear();}
void ClearEventData() {m_EventData->Clear();}
-
- // methods related to the TMinosSpectra class
- // instantiate the TMinosSpectra class and
// declare list of histograms
- void InitSpectra();
-
- // fill the spectra
- void FillSpectra();
-
- // used for Online mainly, sanity check for histograms and
- // change their color if issues are found, for example
- void CheckSpectra();
-
- // used for Online only, clear all the spectra
- void ClearSpectra();
-
- // write spectra to ROOT output file
- void WriteSpectra();
-
//////////////////////////////////////////////////////////////
// specific methods to Minos array
public:
+
+ NPL::Tracking* Tracking_functions; //!
+
+ /* int NclusterFit; //! */
+ // fit function for the Q(t) signals at each pad
+ static double conv_fit(double *x, double *p);
+ static double distance2(double x,double y,double z, double *p);
+ static void SumDistance(int &, double *, double & sum, double * par, int);
+
+
// remove bad channels, calibrate the data and apply thresholds
void PreTreat();
// clear the pre-treated object
- void ClearPreTreatedData() {m_PreTreatedData->Clear();}
+ void ClearPreTreatedData() {
+
+ m_PreTreatedData->Clear();
+ Q_Pad.clear();
+ X_Pad.clear();
+ Y_Pad.clear();
+ Z_Pad.clear();
+ T_Pad.clear();
+ }
// read the user configuration file. If no file is found, load standard one
void ReadAnalysisConfig();
@@ -141,40 +143,182 @@ class TMinosPhysics : public TObject, public NPL::VDetector {
// give and external TMinosData object to TMinosPhysics.
// needed for online analysis for example
void SetRawDataPointer(TMinosData* rawDataPointer) {m_EventData = rawDataPointer;}
-
+
// objects are not written in the TTree
private:
TMinosData* m_EventData; //!
TMinosData* m_PreTreatedData; //!
TMinosPhysics* m_EventPhysics; //!
+ TMinuit* min;//!
+
+TClonesArray data_result;//!
+TMinosResult *minosdata_result;//!
+
+ TClonesArray fitdata;//!
// getters for raw and pre-treated data object
public:
TMinosData* GetRawData() const {return m_EventData;}
TMinosData* GetPreTreatedData() const {return m_PreTreatedData;}
+ vector<double> GetPad_X() {return X_Pad;}
+ vector<double> GetPad_Y() {return Y_Pad;}
+ vector<double> GetPad_Z() {return Z_Pad;}
+ vector<double> GetParFit1() {return parFit1;}
+ vector<double> GetParFit2() {return parFit2;}
+ vector<double> GetParFit3() {return parFit3;}
+ vector<double> GetParFit4() {return parFit4;}
+ double GetVertexZ() {return Zvertex[0];}
+
+ /* TClonesArray fitdata; */
+
+ /* fitdata.SetClass("TMinosClust"); */
+
// parameters used in the analysis
private:
+
+
+ bool SimulationBool=false;//!
+
// thresholds
double m_E_RAW_Threshold; //!
double m_E_Threshold; //!
+ // used in PreTreatedData()
+ vector<int> clusterringboolTemp;//!
+ vector<int> clusterringbool;//!
+ vector<int> clusternbr;//!
+ vector<int> clusterpads;//!
+
+ vector<double> Xpad;//!
+ vector<double> Ypad;//!
+ vector<double> Qpad;//!
+ vector<double> XpadNew;
+ vector<double> YpadNew;
+ vector<double> QpadNew;//!
+ vector<double> ZpadNew;//!
+
+ vector<double> XpadTemp;//!
+ vector<double> YpadTemp;//!
+ vector<double> QpadTemp;//!
+ vector<double> ZpadTemp;//!
+
+ double x_mm;//!
+ double y_mm;//!
+ double z_mm;//!
+ double q_pad;//!
+ double t_pad;//!
+ double Chi2;//!
+ double maxCharge;//!
+ double ChargeBin;//!
+
+ double hfit_max; //!
+ double hfit_max_T; //!
+ double T_min; //!
+ double T_max; //!
+
+
+ int Iteration;///
+ int filter_result;//!
+ int filled;//!
+ int fit2DStatus;//!
+ int indexfill;//!
+ int trackNbr;//!
+
+ TF1 *fit_function;//!
+ /* TH1F *hfit = new TH1F("hfit","hfit",512,0,512);//! */
+ TH1F *hfit ;//!
+
+ int npoint_temp=0;//!
+ int cluster_temp;//!
+ int ringsum;//!
+ int ringtouch[19];//!
+
+ TGraph *gryztmp;//!
+ TGraph* grxztmp;//!
+
+ double zmax;//!
+ vector<double> xin;//!
+ vector<double> yin;//!
+ vector<double> zin;//!
+ vector<double> qin;//!
+ vector<double> xout;//!
+ vector<double> yout;//!
+ vector<double> zout;//!
+ vector<double> qout;//!
+ vector<double> xoutprime;//!
+ vector<double> youtprime;//!
+ vector<double> zoutprime;//!
+ vector<double> qoutprime;//!
+
+ vector<double> TOTxoutprime;
+ vector<double> TOTyoutprime;
+ vector<double> TOTzoutprime;
+ vector<double> TOTqout;
+
+ vector<double> Xvertex;
+ vector<double> Yvertex;
+ vector<double> Zvertex;
+
+ vector<double> sumTheta;
+ vector<double> Theta1;
+ vector<double> Theta2;
+
+ vector<int> trackclusternbr;//!
+ vector<int> tracknbr;//!
+
+ double Tot_tracks=0; //!
+ int trackNbr_FINAL;//!
+ vector<TGraph> gryz;//!
+ vector<TGraph> grxz;//!
+
+ int npoint;//!
+ int array_final;//!
+ TVector3 point;//!
+
+ // Variables only filled when trackNbr_FINAL>=1
+ int allevt_2pfiltered;//!
+ Double_t pStart[4];//!
+ double parFit_temp[4];//!
+ double err_temp[4];//!
+ Double_t chi[2];//!
+ Int_t fitStatus[2];//!
+ Double_t arglist[10];//!
+ Int_t iflag;//!
+ int nvpar;//!
+ int nparx;//!
+ double amin;//!
+ double edm;//!
+ double errdef;//!
+
+ vector<double> lenght;
+ vector<double> chargeTot;
+ vector<double> parFit1;
+ vector<double> parFit2;
+ vector<double> parFit3;
+ vector<double> parFit4;
+ vector<double> errFit1_local;
+ vector<double> errFit2_local;
+ vector<double> errFit3_local;
+ vector<double> errFit4_local;
+
+
+ double xv;//!
+ double yv;//!
+ double zv;//!
+
+ double Dist_min;//!
+ double Theta_tr1;//!
+ double Theta_tr2;//!
// number of detectors
private:
int m_NumberOfDetectors; //!
- // spectra class
- private:
- TMinosSpectra* m_Spectra; // !
-
- // spectra getter
- public:
- map<string, TH1*> GetSpectra();
-
// Static constructor to be passed to the Detector Factory
public:
static NPL::VDetector* Construct();
ClassDef(TMinosPhysics,1) // MinosPhysics structure
};
+
#endif
diff --git a/NPLib/Detectors/Minos/TMinosResult.cxx b/NPLib/Detectors/Minos/TMinosResult.cxx
new file mode 100644
index 000000000..6ec3a453f
--- /dev/null
+++ b/NPLib/Detectors/Minos/TMinosResult.cxx
@@ -0,0 +1,37 @@
+//----------! Written by C. Santamaria / CEA Saclay !----------
+#include "TMinosResult.h"
+
+ClassImp(TMinosResult)
+
+TMinosResult::TMinosResult()
+: TObject(), x_mm(0), y_mm(0), z_mm(0), Chargemax(0), n_Cluster(0), n_Pads(0), z_max(0){
+// Constructor
+}
+TMinosResult::TMinosResult(Double_t xmm, Double_t ymm, Double_t zmm , Double_t chargemax, Double_t ncluster, Double_t npads, Double_t zmax)
+ : TObject()
+{
+ // Constructor
+ x_mm = xmm;
+ y_mm= ymm;
+ z_mm = zmm;
+ Chargemax = chargemax;
+ n_Cluster=ncluster;
+ n_Pads = npads;
+ z_max = zmax;
+}
+TMinosResult::~TMinosResult()
+{
+// Destructor
+}
+
+void TMinosResult::Set(Double_t xmm, Double_t ymm, Double_t zmm, Double_t chargemax, Double_t ncluster, Double_t npads, Double_t zmax)
+{
+ x_mm = xmm;
+ y_mm= ymm;
+ z_mm = zmm;
+ Chargemax = chargemax;
+ n_Cluster=ncluster;
+ n_Pads = npads;
+ z_max = zmax;
+}
+
diff --git a/NPLib/Detectors/Minos/TMinosResult.h b/NPLib/Detectors/Minos/TMinosResult.h
new file mode 100644
index 000000000..38e04f1e9
--- /dev/null
+++ b/NPLib/Detectors/Minos/TMinosResult.h
@@ -0,0 +1,25 @@
+//----------! Written by C. Santamaria / CEA Saclay !----------
+#ifndef TMINOSRESULT_H
+#define TMINOSRESULT_H
+
+//#include "TClonesArray.h"
+#include "TObject.h"
+
+class TMinosResult : public TObject {
+ public:
+ TMinosResult();
+ TMinosResult(Double_t xmm, Double_t ymm, Double_t zmm , Double_t chargemax, Double_t ncluster, Double_t npads, Double_t zmax);
+ void Set(Double_t xmm, Double_t ymm, Double_t zmm , Double_t chargemax, Double_t ncluster, Double_t npads, Double_t zmax);
+ virtual ~TMinosResult();
+ public:
+ Double_t x_mm;
+ Double_t y_mm;
+ Double_t z_mm;
+ Double_t Chargemax;
+ Double_t n_Cluster;
+ Double_t n_Pads;
+ Double_t z_max;
+
+ ClassDef(TMinosResult,1);
+};
+#endif // end of #ifndef TMINOSRESULT_H
diff --git a/NPLib/TrackReconstruction/CMakeLists.txt b/NPLib/TrackReconstruction/CMakeLists.txt
index 4ae6ea5f7..7b71a04dc 100644
--- a/NPLib/TrackReconstruction/CMakeLists.txt
+++ b/NPLib/TrackReconstruction/CMakeLists.txt
@@ -2,7 +2,9 @@ add_custom_command(OUTPUT NPRansacDict.cxx COMMAND ../scripts/build_dict.sh NPRa
add_custom_command(OUTPUT NPClusterDict.cxx COMMAND ../scripts/build_dict.sh NPCluster.h NPClusterDict.cxx NPCluster.rootmap libNPTrackReconstruction.so NPTrackReconstructionLinkDef.h DEPENDS NPCluster.h)
-add_library(NPTrackReconstruction SHARED NPRansac.cxx NPCluster.cxx NPTrack.cxx NPRansacDict.cxx NPClusterDict.cxx)
+add_custom_command(OUTPUT TrackingDict.cxx COMMAND ../scripts/build_dict.sh Tracking.h TrackingDict.cxx Tracking.rootmap libNPTrackReconstruction.so NPTrackReconstructionLinkDef.h DEPENDS Tracking.h)
+
+add_library(NPTrackReconstruction SHARED NPRansac.cxx NPCluster.cxx NPTrack.cxx Tracking.cxx NPRansacDict.cxx NPClusterDict.cxx TrackingDict.cxx)
target_link_libraries(NPTrackReconstruction ${ROOT_LIBRARIES} NPCore)
-install(FILES NPRansac.h NPCluster.h NPTrack.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
+install(FILES NPRansac.h NPCluster.h NPTrack.h Tracking.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
diff --git a/NPLib/TrackReconstruction/NPTrackReconstructionLinkDef.h b/NPLib/TrackReconstruction/NPTrackReconstructionLinkDef.h
index 485465f52..8fc0faa4b 100644
--- a/NPLib/TrackReconstruction/NPTrackReconstructionLinkDef.h
+++ b/NPLib/TrackReconstruction/NPTrackReconstructionLinkDef.h
@@ -1,4 +1,5 @@
#ifdef __CINT__
#pragma link C++ defined_in "./NPRansac.h";
#pragma link C++ defined_in "./NPCluster.h";
+#pragma link C++ defined_in "./Tracking.h";
#endif
diff --git a/NPLib/TrackReconstruction/Tracking.cxx b/NPLib/TrackReconstruction/Tracking.cxx
new file mode 100644
index 000000000..8a93d564e
--- /dev/null
+++ b/NPLib/TrackReconstruction/Tracking.cxx
@@ -0,0 +1,488 @@
+//----------! Developed by C. Santamaria / CEA Saclay !----------
+//----------! Version date :: 2014/11/12 !----------
+// Tracking Functions
+
+#include <fstream>
+#include <string>
+#include <iostream>
+#include <vector>
+#include "TTree.h"
+#include "TFile.h"
+#include "TClonesArray.h"
+#include <vector>
+#include "TF1.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TMath.h"
+#include "TGraph.h"
+// #include "/opt/root/v5.34.14/include/Math/Vector3D.h"
+#include "/opt/cern/root/root_v5.34.36/math/genvector/inc/Math/Vector3D.h"
+/* #include "TMinosResult.h" */
+/* #include "Vector3D.h" */
+#include "Tracking.h"
+
+using namespace std;
+
+/* TClonesArray data_result; */
+/* TMinosResult *minosdata_result; */
+
+ClassImp(NPL::Tracking)
+
+NPL::Tracking::Tracking()
+{
+// Constructor
+}
+
+NPL::Tracking::~Tracking()
+{
+// Destructor
+}
+
+// Fit function for E(T)
+/* double NPL::Tracking::conv_fit(double *x, double *p) { */
+/* double val; */
+/* if(!(x[0]<p[1] || x[0]>512.)) val = p[0] * exp(-3.*(x[0]-p[1])/p[2]) * sin((x[0]-p[1])/p[2]) * pow((x[0]-p[1])/p[2], 3) + 250; */
+/* //else val = p[3]; */
+/* else val = 250; */
+/* return(val); */
+/* } */
+
+int NPL::Tracking::Hough_modified(vector<double> *x,vector<double> *y,vector<double> *q, vector<double> *x_out,vector<double> *y_out, vector<double> *q_out, vector<int> *ringbool) {
+ double bint1=2.;
+ double bint2=2.;
+ int maxt = 360.;
+ int mint = 0.;
+ int nt1=(maxt-mint)/bint1;
+ int nt2=(maxt-mint)/bint1;
+
+ double PI = TMath::Pi();
+ double Rint = 45.2;
+ double Rext = 45.2 + 18*2.1;
+ int filter_result = 0;
+
+ TH2F *hp_xy = new TH2F("hp_xy","hp_xy",nt1,mint,maxt,nt2,mint,maxt);
+ TH2F *hpDiag_xy = new TH2F("hpDiag_xy","hpDiag_xy",nt1,mint,maxt,nt2,mint,maxt);
+
+ double max_xy;
+
+ vector<double> xTemp, yTemp, qTemp;
+
+ double theta1, theta2, xt1, yt1, xt2, yt2;
+ double line0=0., line1=0.;
+ double delta=0., AA=0., BB=0., CC=0.;
+ double maxtheta1=0., maxtheta2=0., xmax1=0., ymax1=0., xmax2=0., ymax2=0.;
+ double par0=0., par1=0.;
+ double r_mm=0.;
+ int ringsum=0;
+ bool maxfound = false;
+
+ for(unsigned int i=0;i<x->size();i++){
+ xTemp.push_back(x->at(i));
+ yTemp.push_back(y->at(i));
+ qTemp.push_back(q->at(i));
+
+ //Loop of indices
+ for(int j=0;j<nt1;j++){
+
+ theta1 = (j+0.5)*bint1 + mint;
+ xt1 = Rint * TMath::Cos(theta1*PI/180.);
+ yt1 = Rint * TMath::Sin(theta1*PI/180.);
+ line1 = (yt1 - y->at(i))/(xt1 - x->at(i));
+ line0 = yt1 - xt1 * line1;
+ AA = 1 + line1*line1;
+ BB = 2*line0*line1;
+ CC = line0*line0 - Rext*Rext;
+
+ delta = BB*BB - 4*AA*CC;
+
+ if(delta>=0){
+ xt2 = (-BB - sqrt(delta))/(2*AA);
+ yt2 = line0 + line1*xt2;
+
+ if(xt2<=0) theta2= 180 - asin(yt2/Rext)*180/PI;
+
+ else if(xt2>0){
+ if(yt2>0) theta2= asin(yt2/Rext)*180/PI;
+ else if(yt2<=0) theta2= 360 + asin(yt2/Rext)*180/PI;
+ }
+
+ if( (xt1*x->at(i) + yt1*y->at(i))>=0 && (xt2*x->at(i) + yt2*y->at(i))>=0 && (xt1*xt2+yt1*yt2)>=0){
+ hp_xy->Fill(theta1,theta2);
+ if(abs(theta1-theta2)<=10) hpDiag_xy->Fill(theta1,theta2);
+ }
+ else{
+ if(delta!=0){
+ xt2 = (-BB + sqrt(delta))/(2*AA);
+ yt2 = line0 + line1*xt2;
+
+ if(xt2<=0) theta2= 180 - asin(yt2/Rext)*180/PI;
+
+ else if(xt2>0){
+ if(yt2>0) theta2= asin(yt2/Rext)*180/PI;
+ else if(yt2<=0) theta2= 360 + asin(yt2/Rext)*180/PI;
+ }
+
+ if( (xt1*x->at(i) + yt1*y->at(i))>=0 && (xt2*x->at(i) + yt2*y->at(i))>=0 && (xt1*xt2+yt1*yt2)>=0){
+ hp_xy->Fill(theta1,theta2);
+ if(abs(theta1-theta2)<=10) hpDiag_xy->Fill(theta1,theta2);
+ }
+ }
+ }
+ } // end if delta>=0
+ } // end loop on indices
+ } // end loop on pads
+
+ x->clear();
+ y->clear();
+ q->clear();
+
+ if(hpDiag_xy->GetMaximum()>=10) max_xy = hpDiag_xy->GetMaximum();
+// cout << "Max taken in diag... withh value=" << max_xy << endl;
+ else max_xy = hp_xy->GetMaximum();
+
+ for(int ii=0; ii<nt1; ii++){
+ if(maxfound ==true) break;
+ for(int jj=0; jj<nt2; jj++){
+ if(hp_xy->GetBinContent(ii+1, jj+1) == max_xy){
+ maxtheta1 = (ii+0.5)*bint1 + mint;
+ maxtheta2 = (jj+0.5)*bint2 + mint;
+ maxfound = true;
+ //cout << "xy: theta max are " << maxtheta1 << " , " << maxtheta2 << endl;
+ }
+ if(maxfound ==true) break;
+ }
+ }
+
+ xmax1 = Rint * TMath::Cos(maxtheta1*PI/180.);
+ ymax1 = Rint * TMath::Sin(maxtheta1*PI/180.);
+ xmax2 = Rext * TMath::Cos(maxtheta2*PI/180.);
+ ymax2 = Rext * TMath::Sin(maxtheta2*PI/180.);
+
+ // xy PEAK
+ par1 = (ymax2-ymax1)/(xmax2-xmax1);
+ par0 = (ymax1 - xmax1*par1);
+ //Selection of x,y points IN the maxmean+/-1 found in Obertelli transform of xy plane
+ for(unsigned int i=0;i<xTemp.size();i++){
+ if( (abs(par1*xTemp[i]-yTemp[i]+par0)/sqrt(1+par1*par1))<= 6 && ((xmax1*xTemp[i] + ymax1*yTemp[i]) >= 0) && ((xmax2*xTemp[i] + ymax2*yTemp[i]) >= 0) && ((xmax1*xmax2 + ymax1*ymax2) >= 0)){
+// couti<< "Taken points= " << xTemp[i] << " , " << yTemp[i] << " , " << zTemp[i] << endl;
+// hcnew_xy->Fill(xTemp[i],yTemp[i],qTemp[i]);
+ x_out->push_back(xTemp[i]);
+ y_out->push_back(yTemp[i]);
+ q_out->push_back(qTemp[i]);
+ filter_result++;
+ /* r_mm = sqrt(xTemp[i]*xTemp[i]+yTemp[i]*yTemp[i]); */
+ /* if(r_mm<(45.2+5*2.1)) ringsum++; */ // commented by Cyril
+
+ }
+ else{
+ x->push_back(xTemp[i]);
+ y->push_back(yTemp[i]);
+ q->push_back(qTemp[i]);
+
+ }
+ }
+
+
+ /* cout << "ringsum = " << ringsum << endl; */
+ /* cout << "Filter_result = " << filter_result << endl; */
+ for(int ip=0; ip<filter_result; ip++){
+ /* if(ringsum>2) ringbool->push_back(1); */
+ /* else ringbool->push_back(0); */
+ ringbool->push_back(1);
+ }
+
+ delete hp_xy;
+ delete hpDiag_xy;
+
+ return filter_result;
+
+}
+
+double NPL::Tracking::FitFunction(double *x, double *p) {
+ double val=p[0]+p[1]*x[0];
+ return(val);
+}
+void NPL::Tracking::FindStart(double pStart[4], double chi[2], int fitStatus[2],TGraph *grxz, TGraph *gryz) {
+ double par1D[2];
+ TF1 *myfit1 = new TF1("myfit1","[0]+[1]*x", -100,500);
+ myfit1->SetParameter(0,0);
+ myfit1->SetParameter(1,10);
+ fitStatus[0] =0;
+ grxz->Fit(myfit1,"RQM");
+ chi[0]=myfit1->GetChisquare();
+ par1D[0]=myfit1->GetParameter(0);
+ par1D[1]=myfit1->GetParameter(1);
+ pStart[0]=par1D[0];
+ pStart[1]=par1D[1];
+ fitStatus[1] =0;
+ gryz->Fit(myfit1,"RQM");
+ chi[1]=myfit1->GetChisquare();
+ par1D[0]=myfit1->GetParameter(0);
+ par1D[1]=myfit1->GetParameter(1);
+ pStart[2]=par1D[0];
+ pStart[3]=par1D[1];
+ //AC 07/12/14
+ delete myfit1;
+}
+
+// Calculation of the distance line-point
+double NPL::Tracking::distance2(double x,double y,double z, double *p) {
+ // distance line point is D= | (xp-x0) cross ux |
+ // where ux is direction of line and x0 is a point in the line (like t = 0)
+ ROOT::Math::XYZVector xp(x,y,z); //point of the track
+ ROOT::Math:: XYZVector x0(p[0], p[2], 0. );
+ ROOT::Math::XYZVector x1(p[0] + p[1], p[2] + p[3], 1. ); //line
+ ROOT::Math::XYZVector u = (x1-x0).Unit();
+ double d2 = ((xp-x0).Cross(u)) .Mag2();
+ return d2;
+}
+
+//void Hough_3D(TCanvas *c1, vector<double> *x,vector<double> *y,vector<double> *z,vector<double> *q, vector<double> *x_out,vector<double> *y_out,vector<double> *z_out,vector<double> *q_out) {
+void NPL::Tracking::Hough_3D(vector<double> *x,vector<double> *y,vector<double> *z,vector<double> *q,vector<double> *x_out,vector<double> *y_out,vector<double> *z_out,vector<double> *q_out) {
+
+ int nt_xy=180;
+ int nt_xz=180;
+ int nt_yz=180;
+ int nr_xy=45;
+ int nr_xz=300;
+ int nr_yz=300;
+ double bint_xy=2.;
+ double bint_xz=2.;
+ double bint_yz=2.;
+ double binr_xy=3.;
+ double binr_xz=3.;
+ double binr_yz=3.;
+ int nt,nr;
+ double PI = TMath::Pi();
+
+ double rho_xy,rho_xz,rho_yz;
+ double theta_xy,theta_xz,theta_yz;
+
+ TH2F *hp_xy = new TH2F("hp_xy","hp_xy",nt_xy,0,180,nr_xy,-1*nr_xy,nr_xy);
+ TH2F *hp_xz = new TH2F("hp_xz","hp_xz",nt_xz,0,180,nr_xz,-1*nr_xz,nr_xz);
+ TH2F *hp_yz = new TH2F("hp_yz","hp_yz",nt_yz,0,180,nr_yz,-1*nr_yz,nr_yz);
+
+ // int npeaks_xy, npeaks_xz, npeaks_yz;
+ vector<double> thetapeaks_xy, rpeaks_xy, thetapeaks_xz, rpeaks_xz, thetapeaks_yz, rpeaks_yz;
+ double max_xy, max_xz, max_yz;
+ double rmean_xy=0, thetamean_xy=0, rmean_xz=0, thetamean_xz=0, rmean_yz=0, thetamean_yz=0;
+
+ double r0_xy=0., r0_xz=0., r0_yz=0., rmin_xy=0., rmin_xz=0., rmin_yz=0., rmax_xy=0., rmax_xz=0., rmax_yz=0.;
+ double tmin=0., tmax=0.;
+ double rinf=0., rsup=0.;
+
+ nt=nt_xy;
+ nr=nr_xy;
+ if(nt<nt_xz)nt=nt_xz;
+ if(nr<nr_xz)nr=nr_xz;
+ if(nt<nt_yz)nt=nt_yz;
+ if(nr<nr_yz)nr=nr_yz;
+
+ for(unsigned int i=0;i<x->size();i++)
+ {
+ //Fill coordinate space histograms for plots
+ //Loop of indices and fill Histograms
+ for(int j=0;j<nt;j++)
+ {
+
+ //xy
+ theta_xy = j*180./nt_xy;
+
+ /* cout <<y->at(i) << endl; */
+ rho_xy = x->at(i)*TMath::Cos(theta_xy*PI/180.)+y->at(i)*TMath::Sin(theta_xy*PI/180.);
+ if(abs(theta_xy)<180. && abs(rho_xy)<nr_xy)
+ {
+ hp_xy->Fill(theta_xy,rho_xy);
+ }
+
+ //xz
+ theta_xz = j*180./nt_xz;
+ rho_xz = z->at(i)*TMath::Cos(theta_xz*PI/180.)+x->at(i)*TMath::Sin(theta_xz*PI/180.);
+ if(abs(theta_xz)<180. && abs(rho_xz)<nr_xz)
+ {
+ hp_xz->Fill(theta_xz,rho_xz);
+ }
+
+ //yz
+ theta_yz = j*180./nt_yz;
+ rho_yz = z->at(i)*TMath::Cos(theta_yz*PI/180.)+y->at(i)*TMath::Sin(theta_yz*PI/180.);
+ if(abs(theta_yz)<180. && abs(rho_yz)<nr_yz)
+ {
+ hp_yz->Fill(theta_yz,rho_yz);
+ }
+ }
+ }
+
+ max_xy = hp_xy->GetMaximum();
+ max_xz = hp_xz->GetMaximum();
+ max_yz = hp_yz->GetMaximum();
+
+ for(int ii=0; ii<nt; ii++)
+ {
+ for(int jj=0; jj<nr; jj++)
+ {
+ if(hp_xy->GetBinContent(ii+1, jj+1) == max_xy && jj<nr_xy)
+ {
+ thetapeaks_xy.push_back((ii+0.5)*nt_xy/nt);
+ rpeaks_xy.push_back((jj+0.5)*2 - nr_xy);
+ rmean_xy += rpeaks_xy.back();
+ thetamean_xy += thetapeaks_xy.back();
+ }
+ if(hp_xz->GetBinContent(ii+1, jj+1) == max_xz)
+ {
+ thetapeaks_xz.push_back((ii+0.5)*nt_xz/nt);
+ rpeaks_xz.push_back((jj+0.5)*2 - nr_xz);
+ rmean_xz += rpeaks_xz.back();
+ thetamean_xz += thetapeaks_xz.back();
+ }
+ if(hp_yz->GetBinContent(ii+1, jj+1) == max_yz)
+ {
+ thetapeaks_yz.push_back((ii+0.5)*nt_yz/nt);
+ rpeaks_yz.push_back((jj+0.5)*2 - nr_yz);
+ rmean_yz += rpeaks_yz.back();
+ thetamean_yz += thetapeaks_yz.back();
+ }
+ }
+ }
+
+ // xy PEAK
+ rmean_xy = rmean_xy / rpeaks_xy.size();
+ thetamean_xy = thetamean_xy / thetapeaks_xy.size();
+
+ // xz PEAK
+ rmean_xz = rmean_xz / rpeaks_xz.size();
+ thetamean_xz = thetamean_xz / thetapeaks_xz.size();
+
+ // yz PEAK
+ rmean_yz = rmean_yz / rpeaks_yz.size();
+ thetamean_yz = thetamean_yz / thetapeaks_yz.size();
+
+ rmean_xy = rpeaks_xy[0];
+ thetamean_xy = thetapeaks_xy[0];
+ rmean_xz = rpeaks_xz[0];
+ thetamean_xz = thetapeaks_xz[0];
+ rmean_yz = rpeaks_yz[0];
+ thetamean_yz = thetapeaks_yz[0];
+
+ //Selection of x,y,z points COMMON to the 3 maxmean+/-1 found in Hough spaces for xy, xz and yz spaces
+ for(unsigned int i=0;i<x->size();i++)
+ {
+ r0_xy = x->at(i)*TMath::Cos(thetamean_xy*PI/180.)+y->at(i)*TMath::Sin(thetamean_xy*PI/180.);
+ tmin = thetamean_xy-bint_xy;
+ tmax = thetamean_xy+bint_xy;
+ if((tmin)<0) tmin = tmin + 180.;
+ if((tmax)>180) tmax = tmax - 180.;
+ rmin_xy = x->at(i)*TMath::Cos(tmin*PI/180.)+y->at(i)*TMath::Sin(tmin*PI/180.);
+ rmax_xy = x->at(i)*TMath::Cos(tmax*PI/180.)+y->at(i)*TMath::Sin(tmax*PI/180.);
+
+ rinf = min( rmean_xy - binr_xy, rmean_xy + binr_xy);
+ rsup = max( rmean_xy - binr_xy, rmean_xy + binr_xy);
+ if((r0_xy>=rinf || rmin_xy>=rinf || rmax_xy>=rinf) && (r0_xy<=rsup || rmin_xy<=rsup || rmax_xy<=rsup))
+ {
+ r0_xz = z->at(i)*TMath::Cos(thetamean_xz*PI/180.)+x->at(i)*TMath::Sin(thetamean_xz*PI/180.);
+ tmin = thetamean_xz-bint_xz;
+ tmax = thetamean_xz+bint_xz;
+ if((tmin)<0) tmin = tmin + 180.;
+ if((tmax)>180) tmax = tmax - 180.;
+ rmin_xz = z->at(i)*TMath::Cos(tmin*PI/180.)+x->at(i)*TMath::Sin(tmin*PI/180.);
+ rmax_xz = z->at(i)*TMath::Cos(tmax*PI/180.)+x->at(i)*TMath::Sin(tmax*PI/180.);
+
+ rinf = min( rmean_xz - binr_xz, rmean_xz + binr_xz);
+ rsup = max( rmean_xz - binr_xz, rmean_xz + binr_xz);
+
+ if((r0_xz>=rinf || rmin_xz>=rinf || rmax_xz>=rinf) && (r0_xz<=rsup || rmin_xz<=rsup || rmax_xz<=rsup))
+ {
+ r0_yz = z->at(i)*TMath::Cos(thetamean_yz*PI/180.)+y->at(i)*TMath::Sin(thetamean_yz*PI/180.);
+ tmin = thetamean_yz-bint_yz;
+ tmax = thetamean_yz+bint_yz;
+ if((tmin)<0) tmin = tmin + 180.;
+ if((tmax)>180) tmax = tmax - 180.;
+ rmin_yz = z->at(i)*TMath::Cos(tmin*PI/180.)+y->at(i)*TMath::Sin(tmin*PI/180.);
+ rmax_yz = z->at(i)*TMath::Cos(tmax*PI/180.)+y->at(i)*TMath::Sin(tmax*PI/180.);
+
+ rinf = min( rmean_yz - binr_yz, rmean_yz + binr_yz);
+ rsup = max( rmean_yz - binr_yz, rmean_yz + binr_yz);
+
+ if((r0_yz>=rinf || rmin_yz>=rinf || rmax_yz>=rinf) && (r0_yz<=rsup || rmin_yz<=rsup || rmax_yz<=rsup))
+ {
+ x_out->push_back(x->at(i));
+ y_out->push_back(y->at(i));
+ z_out->push_back(z->at(i));
+ q_out->push_back(q->at(i));
+ }
+ }
+ }
+
+ }
+ delete hp_xy;
+ delete hp_xz;
+ delete hp_yz;
+
+}
+
+// Calculation of the minimal distance between 2 lines in 3D space & calculation of mid-point=>vertex of interaction
+void NPL::Tracking::vertex(double *p, double *pp, double &xv,double &yv,double &zv,double &min_dist, double &Theta_tr1, double &Theta_tr2) {
+ double a1 = p[0];
+ double a2 = p[2];
+ double b1 = p[1];
+ double b2 = p[3];
+ double ap1 = pp[0];
+ double ap2 = pp[2];
+ double bp1 = pp[1];
+ double bp2 = pp[3];
+
+ double alpha, beta, A, B, C;
+
+ alpha = (bp1*(a1-ap1)+bp2*(a2-ap2))/(bp1*bp1 + bp2*bp2 + 1);
+ beta = (bp1*b1+bp2*b2+1)/(bp1*bp1 + bp2*bp2 + 1);
+
+ A = beta*(bp1*bp1 + bp2*bp2 + 1) - (bp1*b1 + bp2*b2 + 1);
+ B = (b1*b1 + b2*b2 + 1) - beta*(bp1*b1+bp2*b2+1);
+ C = beta*(bp1*(ap1-a1) + bp2*(ap2-a2)) - (b1*(ap1-a1) + b2*(ap2-a2));
+
+ double sol1, solf1;
+ double x,y,z,xp,yp,zp;
+
+
+ sol1 = -(A*alpha + C)/(A*beta + B);
+ solf1 = alpha + beta* sol1;
+
+ x = a1 + b1*sol1;
+ y = a2 + b2*sol1;
+ z = sol1;
+ xp = ap1 + bp1*solf1;
+ yp = ap2 + bp2*solf1;
+ zp = solf1;
+
+ xv = (x+xp)/2.;
+ yv = (y+yp)/2.;
+ zv = (z+zp)/2.;
+
+ Theta_tr1 = TMath::ATan(TMath::Sqrt((x-a1)*(x-a1)+(y-a2)*(y-a2))/TMath::Abs(sol1));
+ Theta_tr2 = TMath::ATan(TMath::Sqrt((xp-ap1)*(xp-ap1)+(yp-ap2)*(yp-ap2))/TMath::Abs(solf1));
+ Theta_tr1 = Theta_tr1*180./TMath::Pi();
+ Theta_tr2 = Theta_tr2*180./TMath::Pi();
+
+ //cout << "Vertex 1st :" << x << "," << y << "," << z << endl;
+ //cout << "Vertex 2nd :" << xp << "," << yp << "," << zp << endl;
+ //cout << "Vertex middle :" << xv << "," << yv << "," << zv << endl;
+
+ //cout << "min dist " << sqrt(pow((x-xp),2) + pow((y-yp),2) + pow((z-zp),2)) << endl;
+
+}
+
+void NPL::Tracking::ParFor_Vertex(double *a, double *b, double *parFit) {
+ // input 2 points in 3D, output parFit format for the use in vertex function
+ double APX, APY, APZ, X0, Y0, APX0, APY0;
+
+ APX=a[0]-b[0]; APY=a[1]-b[1]; APZ=a[2]-b[2];
+ APX0=APX/APZ; APY0=APY/APZ;
+ X0=a[0]-a[2]*APX0; Y0=a[1]-a[2]*APY0; //Z0=0
+
+ parFit[0]=X0;
+ parFit[1]=APX0;
+ parFit[2]=Y0;
+ parFit[3]=APY0;
+}
+
+
diff --git a/NPLib/TrackReconstruction/Tracking.h b/NPLib/TrackReconstruction/Tracking.h
new file mode 100644
index 000000000..ee158a5bb
--- /dev/null
+++ b/NPLib/TrackReconstruction/Tracking.h
@@ -0,0 +1,33 @@
+//----------! Developed by C. Santamaria / CEA Saclay !----------
+//----------! Version date :: 2014/11/12 !----------
+// Tracking Functions
+#ifndef TRACKING_H
+#define TRACKING_H
+
+#include <fstream>
+#include <string>
+#include <iostream>
+#include <vector>
+#include "TGraph.h"
+
+using namespace std;
+
+namespace NPL{
+class Tracking {
+ public:
+ Tracking();
+ virtual ~Tracking();
+ /* double conv_fit(double *x, double *p); */
+ int Hough_modified(vector<double> *x,vector<double> *y,vector<double> *q, vector<double> *x_out,vector<double> *y_out, vector<double> *q_out, vector<int> *ringbool);
+ double FitFunction(double *x, double *p);
+ void FindStart(double pStart[4], double chi[2], int fitStatus[2],TGraph *grxz, TGraph *gryz);
+ double distance2(double x,double y,double z, double *p);
+ void Hough_3D(vector<double> *x,vector<double> *y,vector<double> *z,vector<double> *q,vector<double> *x_out,vector<double> *y_out,vector<double> *z_out,vector<double> *q_out);
+ void vertex(double *p, double *pp, double &xv,double &yv,double &zv, double &min_dist, double &Theta_tr1, double &Theta_tr2);
+ void ParFor_Vertex(double *a, double *b, double *parFit);
+
+ ClassDef(Tracking,1);
+
+};
+}
+#endif
--
GitLab