diff --git a/NPLib/Core/NPDetectorManager.cxx b/NPLib/Core/NPDetectorManager.cxx index 99dd6dd01a8426ee628ead945e1383936485d9c1..1c8274a717eeafc557b8f2cd7ea756cd99d811c9 100644 --- a/NPLib/Core/NPDetectorManager.cxx +++ b/NPLib/Core/NPDetectorManager.cxx @@ -50,6 +50,7 @@ using namespace NPUNITS; // Default Constructor NPL::DetectorManager::DetectorManager(){ + ROOT::EnableThreadSafety(); m_BuildPhysicalPtr = &NPL::VDetector::BuildPhysicalEvent; m_ClearEventPhysicsPtr = &NPL::VDetector::ClearEventPhysics; m_ClearEventDataPtr = &NPL::VDetector::ClearEventData ; diff --git a/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.cxx b/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.cxx index 716841301a1a525cb60141ec17b0115efe84da49..cff7c21ef0d281766f6d39963f63ff16cedea8d3 100644 --- a/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.cxx +++ b/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.cxx @@ -42,7 +42,6 @@ ClassImp(TSamuraiFDC2Physics) /////////////////////////////////////////////////////////////////////////// TSamuraiFDC2Physics::TSamuraiFDC2Physics(){ m_EventData = new TSamuraiFDC2Data ; - m_PreTreatedData = new TSamuraiFDC2Data ; m_EventPhysics = this ; //m_Spectra = NULL; ToTThreshold_L = 180; @@ -62,12 +61,13 @@ void TSamuraiFDC2Physics::BuildPhysicalEvent(){ PreTreat(); // RemoveNoise(); - // Map[detector & plane angle, vector of spatial information] - static map<std::pair<unsigned int,double>, vector<double> > X ; - static map<std::pair<unsigned int,double>, vector<double> > Z ; - static map<std::pair<unsigned int,double>, vector<double> > R ; + // Map[plane angle, vector of spatial information] + static map<double, vector<double> > X ; + static map<double, vector<double> > Z ; + static map<double, vector<double> > R ; static int det,layer,wire; 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){ @@ -75,20 +75,21 @@ void TSamuraiFDC2Physics::BuildPhysicalEvent(){ layer = Layer[i]; wire = Wire[i]; SamuraiDCIndex idx(det,layer,wire); - std::pair<unsigned int, double> p(det,Wire_Angle[idx]); - X[p].push_back(Wire_X[idx]); - Z[p].push_back(Wire_Z[idx]); - R[p].push_back(DriftLength[i]); + X[Wire_Angle[idx]].push_back(Wire_X[idx]); + Z[Wire_Angle[idx]].push_back(Wire_Z[idx]); + R[Wire_Angle[idx]].push_back(DriftLength[i]); } } // Reconstruct the vector for each of the plane of each of the detector static double X0,X100,a,b; // store the BuildTrack2D results - static map<std::pair<unsigned int,double>, TVector3 > VX0 ; - static map<std::pair<unsigned int,double>, TVector3 > VX100 ; - static map<std::pair<unsigned int,double>, double > D ;// the minimum distance - static unsigned int uid=0; + static map<double, TVector3 > VX0 ; + static map<double, TVector3 > VX100 ; + static map<double, double > D ;// the minimum distance + static unsigned int uid; uid=0; + VX0.clear();VX100.clear(),D.clear(); + for(auto it = X.begin();it!=X.end();++it){ #if __cplusplus > 199711L && NPMULTITHREADING m_reconstruction.AddPlan(uid++,X[it->first],Z[it->first],R[it->first]); @@ -98,7 +99,6 @@ void TSamuraiFDC2Physics::BuildPhysicalEvent(){ } #if __cplusplus > 199711L && NPMULTITHREADING - // do all plan at once in parallele, return when all plan are done m_reconstruction.BuildTrack2D(); uid=0; @@ -123,51 +123,51 @@ void TSamuraiFDC2Physics::BuildPhysicalEvent(){ Mult+=X[it->first].size(); // Position at z=0 TVector3 P(X0,0,0); - P.RotateZ(it->first.second); + P.RotateZ(it->first); VX0[it->first]=P; // Position at z=100 TVector3 P100= TVector3(X100,0,0); - P100.RotateZ(it->first.second); + P100.RotateZ(it->first); VX100[it->first]=P100; } // Reconstruct the central position (z=0) for each detector - static map<unsigned int,vector<TVector3> > C ; - static map<unsigned int,vector<double> > W ; // weight based on D + static vector<TVector3> C ; + static vector<double> W ; // weight based on D C.clear(),W.clear(); TVector3 P; for(auto it1 = VX0.begin();it1!=VX0.end();++it1){ for(auto it2 = it1;it2!=VX0.end();++it2){ - if(it1!=it2 && it1->first.first==it2->first.first){// different plane, same detector - m_reconstruction.ResolvePlane(it1->second,it1->first.second,it2->second,it2->first.second,P); + if(it1!=it2){// different plane, same detector + m_reconstruction.ResolvePlane(it1->second,it1->first,it2->second,it2->first,P); - if(P.X()!=-10000 && D[it1->first]<PowerThreshold&& D[it2->first]<PowerThreshold){ - C[it1->first.first].push_back(P); + if(P.X()!=-10000 && D[it1->first]<PowerThreshold && D[it2->first]<PowerThreshold){ + C.push_back(P); // Mean pos are weighted based on the the sum of distance from track // to hit obtained during the minimisation - W[it1->first.first].push_back(1./sqrt(D[it1->first]*D[it2->first])); + W.push_back(1./sqrt(D[it1->first]*D[it2->first])); } } } } // Reconstruct the position at z=100 for each detector - static map<unsigned int,vector<TVector3> > C100 ; + static vector<TVector3> C100 ; C100.clear(); for(auto it1 = VX100.begin();it1!=VX100.end();++it1){ for(auto it2 = it1;it2!=VX100.end();++it2){ - if(it1!=it2 && it1->first.first==it2->first.first){// different plane, same detector - m_reconstruction.ResolvePlane(it1->second,it1->first.second,it2->second,it2->first.second,P); + if(it1!=it2){// different plane + m_reconstruction.ResolvePlane(it1->second,it1->first,it2->second,it2->first,P); if(P.X()!=-10000&& D[it1->first]<PowerThreshold && D[it2->first]<PowerThreshold) - C100[it1->first.first].push_back(P); + C100.push_back(P); } } } // Build the Reference position by averaging all possible pair - size = C[2].size(); + size = C.size(); static double PosX100,PosY100,norm; if(size){ PosX=0; @@ -176,11 +176,11 @@ void TSamuraiFDC2Physics::BuildPhysicalEvent(){ PosY100=0; norm=0; for(unsigned int i = 0 ; i < size ; i++){ - PosX+= C[2][i].X()*W[2][i]; - PosY+= C[2][i].Y()*W[2][i]; - PosX100+= C100[2][i].X()*W[2][i]; - PosY100+= C100[2][i].Y()*W[2][i]; - norm+=W[2][i]; + PosX+= C[i].X()*W[i]; + PosY+= C[i].Y()*W[i]; + PosX100+= C100[i].X()*W[i]; + PosY100+= C100[i].Y()*W[i]; + norm+=W[i]; // cout << C[2][i].X() << " (" << C[2][i].Y() << ") "; } // cout << endl; @@ -195,8 +195,8 @@ void TSamuraiFDC2Physics::BuildPhysicalEvent(){ devX=0; devY=0; for(unsigned int i = 0 ; i < size ; i++){ - devX+=W[2][i]*(C[2][i].X()-PosX)*(C[2][i].X()-PosX); - devY+=W[2][i]*(C[2][i].Y()-PosY)*(C[2][i].Y()-PosY); + devX+=W[i]*(C[i].X()-PosX)*(C[i].X()-PosX); + devY+=W[i]*(C[i].Y()-PosY)*(C[i].Y()-PosY); } devX=sqrt(devX/((size-1)*norm)); @@ -212,19 +212,14 @@ void TSamuraiFDC2Physics::BuildPhysicalEvent(){ PhiY=(PosY100-PosY)/100.; Dir=TVector3(PosX100-PosX,PosY100-PosY,100).Unit(); } - +//cout <<size << " " << PosX << endl; return; } /////////////////////////////////////////////////////////////////////////// void TSamuraiFDC2Physics::PreTreat(){ - ClearPreTreatedData(); static CalibrationManager* Cal = CalibrationManager::getInstance(); static string channel; - // one map per detector - map<unsigned int, vector<double> > X ; - map<unsigned int, vector<double> > Z ; - map<unsigned int, vector<double> > R ; - + //m_EventData->Print(); unsigned int size = m_EventData->Mult(); for(unsigned int i = 0 ; i < size ; i++){ // EDGE=1 is the leading edge, IE, the real time. @@ -264,7 +259,6 @@ void TSamuraiFDC2Physics::PreTreat(){ DriftLength.push_back(10-Cal->ApplySigmoid(channel,etime)); } } - } return; } @@ -400,7 +394,6 @@ void TSamuraiFDC2Physics::InitSpectra(){ /////////////////////////////////////////////////////////////////////////// void TSamuraiFDC2Physics::FillSpectra(){ // m_Spectra -> FillRawSpectra(m_EventData); - // m_Spectra -> FillPreTreatedSpectra(m_PreTreatedData); // m_Spectra -> FillPhysicsSpectra(m_EventPhysics); } /////////////////////////////////////////////////////////////////////////// diff --git a/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.h b/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.h index 80576f85d0f7c176cac4e62079c6bd394bfc5722..5de457b9f4d9a0393d17be4ee14b84023dda0b1f 100644 --- a/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.h +++ b/NPLib/Detectors/Samurai/TSamuraiFDC2Physics.h @@ -158,15 +158,11 @@ class TSamuraiFDC2Physics : public TObject, public NPL::VDetector{ public: // Specific to SamuraiFDC2 Array - // Clear The PreTeated object - void ClearPreTreatedData() {m_PreTreatedData->Clear();} - // Remove bad channel, calibrate the data and apply threshold void PreTreat();//! // Retrieve raw and pre-treated data TSamuraiFDC2Data* GetRawData() const {return m_EventData;} - TSamuraiFDC2Data* GetPreTreatedData() const {return m_PreTreatedData;} double GetPosX(){return PosX;} double GetPosY(){return PosY;} @@ -177,7 +173,6 @@ class TSamuraiFDC2Physics : public TObject, public NPL::VDetector{ private: // Root Input and Output tree classes TSamuraiFDC2Data* m_EventData;//! - TSamuraiFDC2Data* m_PreTreatedData;//! TSamuraiFDC2Physics* m_EventPhysics;//! diff --git a/NPLib/TrackReconstruction/NPDCReconstructionMT.cxx b/NPLib/TrackReconstruction/NPDCReconstructionMT.cxx index 09e8e8c7757406a54864c6130949d6f986d873e5..2d466e04b4630f202d570a9e699948c625c4d22e 100644 --- a/NPLib/TrackReconstruction/NPDCReconstructionMT.cxx +++ b/NPLib/TrackReconstruction/NPDCReconstructionMT.cxx @@ -26,12 +26,14 @@ #include "TError.h" #include "TGraph.h" #include "TVector3.h" +#include "TROOT.h" using namespace std; using namespace NPL; //////////////////////////////////////////////////////////////////////////////// DCReconstructionMT::DCReconstructionMT(unsigned int number_thread){ + ROOT::EnableThreadSafety(); m_nbr_thread= number_thread; // force loading of the minimizer plugin ahead ROOT::Math::Minimizer* mini=ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad"); @@ -146,8 +148,9 @@ double DCReconstructionMT::SumD(const double* parameter ){ double ab= a*b; double a2=a*a; unsigned int id = parameter[2]; + unsigned int size = sizeX[id]; double c,d,r,x,z,p; - for(unsigned int i = 0 ; i < sizeX[id] ; i++){ + for(unsigned int i = 0 ; i < size ; i++){ c = (*fitX[id])[i]; d = (*fitZ[id])[i]; r = (*fitR[id])[i]; @@ -159,7 +162,7 @@ double DCReconstructionMT::SumD(const double* parameter ){ } // return normalized power - return P/sizeX[id]; + return P/size; }