diff --git a/NPLib/Detectors/Actar/TActarPhysics.cxx b/NPLib/Detectors/Actar/TActarPhysics.cxx
index 463ac59dfd7275ccb3bfaafef2a3ecff464e0a1f..ddf0d8a1e7767834a121946162b41b520f5332bf 100644
--- a/NPLib/Detectors/Actar/TActarPhysics.cxx
+++ b/NPLib/Detectors/Actar/TActarPhysics.cxx
@@ -70,6 +70,7 @@ fPadSizeY(2),
fDriftVelocity(40),
fPressure(100),
fGas("iC4H10"),
+m_NumberOfPadSilicon(20),
m_NumberOfDetectors(0) {
}
@@ -99,33 +100,36 @@ void TActarPhysics::BuildSimplePhysicalEvent() {
///////////////////////////////////////////////////////////////////////////
void TActarPhysics::BuildPhysicalEvent() {
-
+
PreTreat();
-
+
//unsigned int mysize = m_PreTreatedData->GetPadMult();
-
+
/*for (unsigned int e = 0; e < mysize; e++) {
PadX.push_back(m_PreTreatedData->GetPadX(e));
PadY.push_back(m_PreTreatedData->GetPadY(e));
PadZ.push_back(m_PreTreatedData->GetPadZ(e));
PadCharge.push_back(m_PreTreatedData->GetPadCharge(e));
}*/
-
if(fRecoRansac && PadX.size()>fHitThreshold){
m_Ransac->Init(PadX, PadY, PadZ, PadCharge);
m_Track = m_Ransac->SimpleRansac();
}
-
- else if(fRecoCluster && PadX.size()>fHitThreshold){
+
+ else if(fRecoCluster){
+ if(PadX.size()>fHitThreshold){
m_Cluster->Init(PadX, PadY, PadZ, PadCharge);
m_Cluster->FindTracks();
-
+ }
+
+ if(BeamPadX.size()>fHitThreshold){
m_Cluster->Init(BeamPadX, BeamPadY, BeamPadZ, BeamPadCharge);
m_Cluster->FindTracks();
-
- m_Track = m_Cluster->GetTracks();
-
- m_Cluster->Clear();
+ }
+
+ m_Track = m_Cluster->GetTracks();
+
+ m_Cluster->Clear();
}
TrackMult = m_Track.size();
@@ -135,25 +139,25 @@ void TActarPhysics::BuildPhysicalEvent() {
void TActarPhysics::PreTreat() {
// This method typically applies thresholds and calibrations
// Might test for disabled channels for more complex detector
-
+
// clear pre-treated object
ClearPreTreatedData();
-
+
CleanPads();
-
+
// instantiate CalibrationManager
static CalibrationManager* Cal = CalibrationManager::getInstance();
-
-
+
+
unsigned int mysize = m_EventReduced->CoboAsad.size();
-
+
for (unsigned int it = 0; it < mysize ; ++it) {
int co=m_EventReduced->CoboAsad[it].globalchannelid>>11;
int as=(m_EventReduced->CoboAsad[it].globalchannelid - (co<<11))>>9;
int ag=(m_EventReduced->CoboAsad[it].globalchannelid - (co<<11)-(as<<9))>>7;
int ch=m_EventReduced->CoboAsad[it].globalchannelid - (co<<11)-(as<<9)-(ag<<7);
int where=co*NumberOfASAD*NumberOfAGET*NumberOfChannel + as*NumberOfAGET*NumberOfChannel + ag*NumberOfChannel + ch;
-
+
if(co!=31){
unsigned int vector_size = m_EventReduced->CoboAsad[it].peakheight.size();
for(unsigned int hh=0; hh<vector_size; hh++){
@@ -173,7 +177,7 @@ void TActarPhysics::PreTreat() {
PadY.push_back(TABLE[5][where]);
PadZ.push_back(m_EventReduced->CoboAsad[it].peaktime[hh]);
}
-
+
}
else if(fRecoRansac){
PadCharge.push_back(m_EventReduced->CoboAsad[it].peakheight[hh]);
@@ -189,14 +193,25 @@ void TActarPhysics::PreTreat() {
else if(co==31){
unsigned int vector_size = m_EventReduced->CoboAsad[it].peakheight.size();
for(unsigned int hit=0;hit<vector_size;hit++){
-
- Si_Number.push_back(m_EventReduced->CoboAsad[it].peaktime[hit]);
- Si_E.push_back(m_EventReduced->CoboAsad[it].peakheight[hit]);
-
- int vxi_parameter = m_EventReduced->CoboAsad[it].peaktime[hit];
- if(Si_map[vxi_parameter]<21 && Si_map[vxi_parameter]>0){
- m_PreTreatedData->SetSiliconEnergy(m_EventReduced->CoboAsad[it].peakheight[hit]);
- m_PreTreatedData->SetSiliconDetectorNumber(Si_map[vxi_parameter]);
+ if(fInputTreeName=="ACTAR_TTree"){
+ int vxi_parameter = m_EventReduced->CoboAsad[it].peaktime[hit];
+
+ if(Si_map[vxi_parameter]<21 && Si_map[vxi_parameter]>0){
+ double RawEnergy = m_EventReduced->CoboAsad[it].peakheight[hit];
+
+ static string name;
+ name = "ActarSi/D" ;
+ name+= NPL::itoa( Si_map[vxi_parameter] - 1) ;
+ name+= "_E" ;
+ double CalEnergy = CalibrationManager::getInstance()->ApplyCalibration( name, RawEnergy );
+
+ Si_Number.push_back(Si_map[vxi_parameter]);
+ Si_E.push_back(CalEnergy);
+ }
+ }
+ else{
+ Si_Number.push_back(m_EventReduced->CoboAsad[it].peaktime[hit]);
+ Si_E.push_back(m_EventReduced->CoboAsad[it].peakheight[hit]);
}
}
}
@@ -210,7 +225,7 @@ bool TActarPhysics::IsBeamZone(int X, int Y)
if( (X>=fXBeamMin && X<=fXBeamMax) && (Y>=fYBeamMin && Y<=fYBeamMax) ){
isBeam=true;
}
-
+
return isBeam;
}
@@ -225,7 +240,7 @@ bool TActarPhysics::GoodHit(int iX, int iY)
}
}
}
-
+
return bHit;
}
@@ -241,8 +256,8 @@ void TActarPhysics::CleanPads()
int ag=(m_EventReduced->CoboAsad[it].globalchannelid - (co<<11)-(as<<9))>>7;
int ch=m_EventReduced->CoboAsad[it].globalchannelid - (co<<11)-(as<<9)-(ag<<7);
int where=co*NumberOfASAD*NumberOfAGET*NumberOfChannel + as*NumberOfAGET*NumberOfChannel + ag*NumberOfChannel + ch;
-
-
+
+
if(co!=31){
unsigned int vector_size = m_EventReduced->CoboAsad[it].peakheight.size();
for(unsigned int hh=0; hh < vector_size; hh++){
@@ -279,7 +294,7 @@ void TActarPhysics::ReadAnalysisConfig() {
}
}
VXIConfigFile.close();
-
+
// Lookup table //
bool ReadingLookupTable = false;
string LT_FileName = "./configs/LT.dat";
@@ -297,23 +312,23 @@ void TActarPhysics::ReadAnalysisConfig() {
ReadingLookupTable = true;
}
LTConfigFile.close();
-
-
+
+
bool ReadingStatus = false;
-
+
// ACTAR config file //
string FileName = "./configs/ConfigActar.dat";
-
+
// open analysis config file
ifstream AnalysisConfigFile;
AnalysisConfigFile.open(FileName.c_str());
-
+
if (!AnalysisConfigFile.is_open()) {
cout << " No ConfigActar.dat found: Default parameter loaded for Analayis " << FileName << endl;
return;
}
cout << "/// Loading user parameter for Analysis from ConfigActar.dat ///" << endl;
-
+
// Save it in a TAsciiFile
TAsciiFile* asciiConfig = RootOutput::getInstance()->GetAsciiFileAnalysisConfig();
asciiConfig->AppendLine("%%% ConfigActar.dat %%%");
@@ -324,54 +339,54 @@ void TActarPhysics::ReadAnalysisConfig() {
while (!AnalysisConfigFile.eof()) {
// Pick-up next line
getline(AnalysisConfigFile, LineBuffer);
-
+
// search for "header"
string name = "ConfigActar";
if (LineBuffer.compare(0, name.length(), name) == 0){
ReadingStatus = true;
cout << "**** ConfigActar found" << endl;
}
-
+
// loop on tokens and data
while (ReadingStatus ) {
whatToDo="";
AnalysisConfigFile >> whatToDo;
-
+
// Search for comment symbol (%)
if (whatToDo.compare(0, 1, "%") == 0) {
AnalysisConfigFile.ignore(numeric_limits<streamsize>::max(), '\n' );
}
-
+
else if (whatToDo.compare(0,11,"RecoRansac=") == 0) {
AnalysisConfigFile >> DataBuffer;
fRecoRansac = atoi(DataBuffer.c_str());
cout << "/// Reco using Ransac= " << " " << fRecoRansac << " ///" << endl;
}
-
+
else if (whatToDo.compare(0,12,"RecoCluster=") == 0) {
AnalysisConfigFile >> DataBuffer;
fRecoCluster = atoi(DataBuffer.c_str());
cout << "/// Reco using Cluster= " << " " << fRecoCluster << " ///" << endl;
}
-
+
else if (whatToDo.compare(0,9,"RecoVisu=") == 0) {
AnalysisConfigFile >> DataBuffer;
fRecoVisu = atoi(DataBuffer.c_str());
cout << "/// Visualisation= " << " " << fRecoVisu << " ///" << endl;
}
-
+
else if (whatToDo.compare(0,14,"HIT_THRESHOLD=") == 0) {
AnalysisConfigFile >> DataBuffer;
fHitThreshold = atoi(DataBuffer.c_str());
cout << "/// Hit Threshold= " << " " << fHitThreshold << " ///" << endl;
}
-
+
else if (whatToDo.compare(0,12,"Q_THRESHOLD=") == 0) {
AnalysisConfigFile >> DataBuffer;
fQ_Threshold = atof(DataBuffer.c_str());
cout << "/// Q Threshold= " << " " << fQ_Threshold << " ///" << endl;
}
-
+
else if (whatToDo.compare(0,12,"T_THRESHOLD=") == 0) {
AnalysisConfigFile >> DataBuffer;
fT_Threshold = atof(DataBuffer.c_str());
@@ -383,79 +398,86 @@ void TActarPhysics::ReadAnalysisConfig() {
//check_padsX=true;
cout << "/// Number Of Pads X= " << fNumberOfPadsX << " ///" << endl;
}
-
+
else if(whatToDo.compare(0,14,"NumberOfPadsY=")==0){
AnalysisConfigFile >> DataBuffer;
fNumberOfPadsY = atoi(DataBuffer.c_str());
//check_padsY=true;
cout << "/// Number Of Pads Y= " << fNumberOfPadsY << " ///" << endl;
}
-
+
else if(whatToDo.compare(0,9,"PadSizeX=")==0){
AnalysisConfigFile >> DataBuffer;
fPadSizeX = atof(DataBuffer.c_str());
//check_sizeX=true;
cout << "/// Pad Size X= " << fPadSizeX << " ///" << endl;
}
-
+
else if(whatToDo.compare(0,9,"PadSizeY=")==0){
AnalysisConfigFile >> DataBuffer;
fPadSizeY = atof(DataBuffer.c_str());
//check_sizeY=true;
cout << "/// Pad Size Y= " << fPadSizeY << " ///" << endl;
}
-
+
else if(whatToDo.compare(0,9,"XBeamMin=")==0){
AnalysisConfigFile >> DataBuffer;
fXBeamMin = atof(DataBuffer.c_str());
cout << "/// X Beam Min= " << fXBeamMin << " ///" << endl;
}
-
+
else if(whatToDo.compare(0,9,"XBeamMax=")==0){
AnalysisConfigFile >> DataBuffer;
fXBeamMax = atof(DataBuffer.c_str());
cout << "/// X Beam Max= " << fXBeamMax << " ///" << endl;
}
-
+
else if(whatToDo.compare(0,9,"YBeamMin=")==0){
AnalysisConfigFile >> DataBuffer;
fYBeamMin = atof(DataBuffer.c_str());
cout << "/// Y Beam Min= " << fYBeamMin << " ///" << endl;
}
-
+
else if(whatToDo.compare(0,9,"YBeamMax=")==0){
AnalysisConfigFile >> DataBuffer;
fYBeamMax = atof(DataBuffer.c_str());
cout << "/// Y Beam Max= " << fYBeamMax << " ///" << endl;
}
-
+
+ else if(whatToDo.compare(0,19,"NumberOfPadSilicon=")==0){
+ AnalysisConfigFile >> DataBuffer;
+ m_NumberOfPadSilicon = atof(DataBuffer.c_str());
+ //check_pressure=true;
+ cout << "/// Number of Silicons= " << m_NumberOfPadSilicon << " ///" << endl;
+ }
+
else if(whatToDo.compare(0,9,"Pressure=")==0){
AnalysisConfigFile >> DataBuffer;
fPressure = atof(DataBuffer.c_str());
//check_pressure=true;
cout << "/// Pressure= " << fPressure << " ///" << endl;
}
-
+
else if(whatToDo.compare(0,14,"DriftVelocity=")==0){
AnalysisConfigFile >> DataBuffer;
fDriftVelocity = atof(DataBuffer.c_str());
//check_driftvelocity=true;
cout << "/// Drift Velocity= " << fDriftVelocity << " ///" << endl;
}
-
+
else if(whatToDo.compare(0,4,"Gas=")==0){
AnalysisConfigFile >> DataBuffer;
fGas = DataBuffer.c_str();
//check_gas=true;
cout << "/// Gas Type= " << fGas << " ///" << endl;
}
-
+
else {
ReadingStatus = false;
}
}
}
-
+
if(fRecoRansac){
m_Ransac = new NPL::Ransac(fNumberOfPadsX,fNumberOfPadsY,fRecoVisu);
}
@@ -484,7 +506,7 @@ void TActarPhysics::Clear() {
BeamPadY.clear();
BeamPadZ.clear();
BeamPadCharge.clear();
-
+
PadX.clear();
PadY.clear();
PadZ.clear();
@@ -492,7 +514,7 @@ void TActarPhysics::Clear() {
m_Track.clear();
Si_Number.clear();
Si_E.clear();
-
+
for(int i=0; i<fNumberOfPadsX; i++){
for(int j=0; j<fNumberOfPadsY; j++){
Hit[i][j] = 0;
@@ -507,15 +529,15 @@ void TActarPhysics::ReadConfiguration(NPL::InputParser parser) {
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Actar");
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 << "//// Actar " << i+1 << endl;
-
+
TVector3 Pos = blocks[i]->GetTVector3("POS","mm");
string Shape = blocks[i]->GetString("Shape");
AddDetector(Pos,Shape);
@@ -604,6 +626,10 @@ void TActarPhysics::AddParameterToCalibrationManager() {
Cal->AddParameter("Actar", "D"+ NPL::itoa(i+1)+"_CHARGE","Actar_D"+ NPL::itoa(i+1)+"_CHARGE");
Cal->AddParameter("Actar", "D"+ NPL::itoa(i+1)+"_TIME","Actar_D"+ NPL::itoa(i+1)+"_TIME");
}
+
+ for(int i=0; i<m_NumberOfPadSilicon; i++){
+ Cal->AddParameter("ActarSi","D"+NPL::itoa(i)+"_E","ActarSi_D"+NPL::itoa(i)+"_E");
+ }
}
@@ -621,6 +647,8 @@ void TActarPhysics::InitializeRootInputRaw() {
inputChain->SetBranchStatus("data", true );
inputChain->SetBranchStatus("ACTAR_TTree", true );
inputChain->SetBranchAddress("data", &m_EventReduced );
+
+ fInputTreeName = inputChain->GetTree()->GetName();
}
@@ -661,6 +689,6 @@ extern "C"{
NPL::DetectorFactory::getInstance()->AddDetector("Actar",TActarPhysics::Construct);
}
};
-
+
proxy_Actar p_Actar;
}
diff --git a/NPLib/Detectors/Actar/TActarPhysics.h b/NPLib/Detectors/Actar/TActarPhysics.h
index a80f629d5f91973d61d8d54a7e6b6e236bc6fbc3..505fbd367072f302b70fb3829fc2b077f910041a 100644
--- a/NPLib/Detectors/Actar/TActarPhysics.h
+++ b/NPLib/Detectors/Actar/TActarPhysics.h
@@ -58,15 +58,15 @@ class TActarPhysics : public TObject, public NPL::VDetector {
public:
TActarPhysics();
~TActarPhysics() {};
-
-
+
+
//////////////////////////////////////////////////////////////
// Inherited from TObject and overriden to avoid warnings
public:
void Clear();
void Clear(const Option_t*) {};
-
-
+
+
//////////////////////////////////////////////////////////////
// data obtained after BuildPhysicalEvent() and stored in
// output ROOT file
@@ -82,94 +82,94 @@ public:
vector<double> Si_E;
vector<int> Si_Number;
int TrackMult;
-
+
/// 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);
-
+
public:
int GetTrackMult() {return m_Track.size();}
vector<NPL::Track> GetTracks() {return m_Track;}
-
-
+
+
//////////////////////////////////////////////////////////////
// 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 TActarSpectra class
// instantiate the TActarSpectra 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 Actar array
public:
// remove bad channels, calibrate the data and apply thresholds
void PreTreat();
-
+
// clear the pre-treated object
void ClearPreTreatedData() {m_PreTreatedData->Clear();}
-
+
// read the user configuration file. If no file is found, load standard one
void ReadAnalysisConfig();
-
+
void CleanPads();
-
+
bool GoodHit(int iX, int iY);
-
+
bool IsBeamZone(int X, int Y);
-
+
// give and external TActarData object to TActarPhysics.
// needed for online analysis for example
void SetRawDataPointer(TActarData* rawDataPointer) {m_EventData = rawDataPointer;}
-
+
// objects are not written in the TTree
private:
TActarData* m_EventData; //!
@@ -179,12 +179,12 @@ private:
NPL::Ransac* m_Ransac; //!
NPL::Cluster* m_Cluster; //!
vector<NPL::Track> m_Track; //!
-
+
// getters for raw and pre-treated data object
public:
TActarData* GetRawData() const {return m_EventData;}
TActarData* GetPreTreatedData() const {return m_PreTreatedData;}
-
+
double GetDriftVelocity() {return fDriftVelocity;}
double GetPadSizeX() {return fPadSizeX;}
double GetPadSizeY() {return fPadSizeY;}
@@ -192,7 +192,7 @@ public:
int GetNumberOfPadsY() {return fNumberOfPadsY;}
double GetPRessure() {return fPressure;}
string GetGasName() {return fGas;}
-
+
// parameters used in the analysis
private:
// thresholds
@@ -214,19 +214,21 @@ private:
bool fRecoCluster; //!
bool fRecoVisu; //!
map<int, int> Si_map; //!
-
+ string fInputTreeName; //!
+
int TABLE[6][NumberOfCobo*NumberOfASAD*NumberOfAGET*NumberOfChannel]; //!
int Hit[128][128]; //!
-
+
// number of detectors
private:
int m_NumberOfDetectors; //!
-
+ int m_NumberOfPadSilicon; //!
+
// spectra class
private:
TActarSpectra* m_Spectra; //!
-
+
//spme getters and setters
public:
void SetRansacParameter(string filename);
@@ -235,17 +237,17 @@ public:
bool GetRansacStatus() {return fRecoRansac;}
NPL::Cluster* GetClusterObject() {return m_Cluster;}
bool GetClusterStatus() {return fRecoCluster;}
-
+
// spectra getter
public:
map<string, TH1*> GetSpectra();
vector<TCanvas*> GetCanvas();
-
+
// Static constructor to be passed to the Detector Factory
public:
static NPL::VDetector* Construct();
-
+
ClassDef(TActarPhysics,1) // ActarPhysics structure
};
#endif
diff --git a/NPLib/TrackReconstruction/NPCluster.cxx b/NPLib/TrackReconstruction/NPCluster.cxx
index 3739e30332667dc4efa8fbddd7ba38301cf0de1a..77421d2d4210a99b1bfbdda2ceb4fb7aae532fef 100644
--- a/NPLib/TrackReconstruction/NPCluster.cxx
+++ b/NPLib/TrackReconstruction/NPCluster.cxx
@@ -34,7 +34,7 @@ Cluster::Cluster()
fClusterThreshold = 15;
fVisu= 0;
-
+
LocBigVoxels = new int[BigVoxelSize];
}
@@ -47,11 +47,11 @@ Cluster::Cluster(int NumberOfPadsX, int NumberOfPadsY, bool Visu)
fNumberOfTracksMax = 20;
fClusterDistance = 5;
fClusterThreshold = 15;
-
+
fVisu = 0;
-
+
LocBigVoxels = new int[BigVoxelSize];
-
+
if(fVisu==1){
m_ServerSocket = new TServerSocket(9092,true,100);
//pl = new TGraph2D();
@@ -63,17 +63,17 @@ Cluster::Cluster(int NumberOfPadsX, int NumberOfPadsY, bool Visu)
cout << endl;
cout << "/// ServerSocket valid ///" << endl;
m_ServerSocket->SetCompressionSettings(1);
-
+
// Add server socket to monitor so we are notified when a client needs to be
// accepted
m_Monitor = new TMonitor();
m_Monitor->Add(m_ServerSocket,TMonitor::kRead|TMonitor::kWrite);
// Create a list to contain all client connections
m_Sockets = new TList;
-
+
// Create a list to contain all client connections
m_Histo = new TList;
-
+
}
}
@@ -87,10 +87,10 @@ Cluster::~Cluster()
void Cluster::ReadParameterValue(string filename)
{
bool ReadingStatus = false;
-
+
ifstream ClusterParamFile;
ClusterParamFile.open(filename.c_str());
-
+
if(!ClusterParamFile.is_open()){
cout << "/// No Paramter File found for Cluster parameters ///" << endl;
cout << "/// Using the default parameter:" << endl;
@@ -113,7 +113,7 @@ void Cluster::ReadParameterValue(string filename)
if (whatToDo.compare(0, 1, "%") == 0) {
ClusterParamFile.ignore(numeric_limits<streamsize>::max(), '\n' );
}
-
+
else if (whatToDo.compare(0,18,"NumberOfTracksMax=") == 0) {
ClusterParamFile >> DataBuffer;
fNumberOfTracksMax = atoi(DataBuffer.c_str());
@@ -141,15 +141,17 @@ void Cluster::ReadParameterValue(string filename)
void Cluster::FindTracks()
{
for(unsigned int i=0; i<fOriginalCloudSize; i++){
+ if(vX[i]>=0 && vX[i]<fNumberOfPadsX && vY[i]>=0 && vY[i]<fNumberOfPadsY){
FillBigVoxels(vX[i],vY[i],(int)vZ[i],(int)vQ[i]);
+ }
}
-
+
RegroupZone();
-
+
FitZone();
-
+
RegroupTrack();
-
+
FitTrack();
}
@@ -163,76 +165,86 @@ void Cluster::FillBigVoxels(int x, int y, int t, int q)
{
BigVoxels[int(x/4 + y/4*fNumberOfPadsX/4 + t/8*fNumberOfPadsX/4*fNumberOfPadsY/4)].push_back(q);
IDBigVoxels[int(x/4 + y/4*fNumberOfPadsX/4 + t/8*fNumberOfPadsX/4*fNumberOfPadsY/4)].push_back(x + y*fNumberOfPadsX + t*fNumberOfPadsX*fNumberOfPadsY);
-
+
for(int tB=TMath::Max(t/8-1.,0.);tB<=TMath::Min(t/8+1.,510/8.);tB++){
for(int yB=TMath::Max(y/4-1.,0.);yB<=TMath::Min(y/4+1.,(fNumberOfPadsY-1)/4.);yB++){
for(int xB=TMath::Max(x/4-1.,0.);xB<=TMath::Min(x/4+1.,(fNumberOfPadsX-1)/4.);xB++){
if(BigVoxels[xB+yB*fNumberOfPadsX/4+tB*fNumberOfPadsX/4*fNumberOfPadsY/4].size()>0 && (xB!=x/4||yB!=y/4||tB!=t/8) && LocBigVoxels[x/4+y/4*fNumberOfPadsX/4+t/8*fNumberOfPadsX/4*fNumberOfPadsY/4]==-1 && LocBigVoxels[xB+yB*fNumberOfPadsX/4+tB*fNumberOfPadsX/4*fNumberOfPadsY/4]==-1){
//(yB+y!=fNumberOfPadsY/4-1) corresponds to the differentiation between each side of the pad plane
// Case where neither the current box nor the neighbour studied are associated to any track zone yet
+ //cout << 1 << endl;
it++;
-
-
+
NPL::Track BigTrack;
BigTrack.SetPointX(x/4);
BigTrack.SetPointY(y/4);
BigTrack.SetPointZ(t/8);
-
+
BigTrack.SetPointX(xB);
BigTrack.SetPointY(yB);
BigTrack.SetPointZ(tB);
-
+
LocBigVoxels[xB+yB*fNumberOfPadsX/4+tB*fNumberOfPadsX/4*fNumberOfPadsY/4]=it;
LocBigVoxels[x/4+y/4*fNumberOfPadsX/4+t/8*fNumberOfPadsX/4*fNumberOfPadsY/4]=it;
TrackZone.push_back(BigTrack);
BigTrack.Clear();
+
+
}
else if(BigVoxels[xB+yB*fNumberOfPadsX/4+tB*fNumberOfPadsX/4*fNumberOfPadsY/4].size()>0 && (xB!=x/4||yB!=y/4||tB!=t/8) && LocBigVoxels[x/4+y/4*fNumberOfPadsX/4+t/8*fNumberOfPadsX/4*fNumberOfPadsY/4]==-1 && LocBigVoxels[xB+yB*fNumberOfPadsX/4+tB*fNumberOfPadsX/4*fNumberOfPadsY/4]>-1){
// Case where the neighbour is already associated with a track zone but not the current box
+ //cout << 2 << endl;
LocBigVoxels[x/4+y/4*fNumberOfPadsX/4+t/8*fNumberOfPadsX/4*fNumberOfPadsY/4]=LocBigVoxels[xB+yB*fNumberOfPadsX/4+tB*fNumberOfPadsX/4*fNumberOfPadsY/4];
TrackZone.at(LocBigVoxels[xB+yB*fNumberOfPadsX/4+tB*fNumberOfPadsX/4*fNumberOfPadsY/4]).SetPointX(x/4);
TrackZone.at(LocBigVoxels[xB+yB*fNumberOfPadsX/4+tB*fNumberOfPadsX/4*fNumberOfPadsY/4]).SetPointY(y/4);
TrackZone.at(LocBigVoxels[xB+yB*fNumberOfPadsX/4+tB*fNumberOfPadsX/4*fNumberOfPadsY/4]).SetPointZ(t/8);
+
+
}
else if(BigVoxels[xB+yB*fNumberOfPadsX/4+tB*fNumberOfPadsX/4*fNumberOfPadsY/4].size()>0 && (xB!=x/4||yB!=y/4||tB!=t/8) && LocBigVoxels[x/4+y/4*fNumberOfPadsX/4+t/8*fNumberOfPadsX/4*fNumberOfPadsY/4]>-1 && LocBigVoxels[xB+yB*fNumberOfPadsX/4+tB*fNumberOfPadsX/4*fNumberOfPadsY/4]==-1){
// Case where the current box is already associated with a track zone but not the neighbour
+ //cout << 3 << endl;
LocBigVoxels[xB+yB*fNumberOfPadsX/4+tB*fNumberOfPadsX/4*fNumberOfPadsY/4]=LocBigVoxels[x/4+y/4*fNumberOfPadsX/4+t/8*fNumberOfPadsX/4*fNumberOfPadsY/4];
TrackZone.at(LocBigVoxels[xB+yB*fNumberOfPadsX/4+tB*fNumberOfPadsX/4*fNumberOfPadsY/4]).SetPointX(xB);
TrackZone.at(LocBigVoxels[xB+yB*fNumberOfPadsX/4+tB*fNumberOfPadsX/4*fNumberOfPadsY/4]).SetPointY(yB);
TrackZone.at(LocBigVoxels[xB+yB*fNumberOfPadsX/4+tB*fNumberOfPadsX/4*fNumberOfPadsY/4]).SetPointZ(tB);
+
+
}
else if(LocBigVoxels[x/4+y/4*fNumberOfPadsX/4+t/8*fNumberOfPadsX/4*fNumberOfPadsY/4]!=LocBigVoxels[xB+yB*fNumberOfPadsX/4+tB*fNumberOfPadsX/4*fNumberOfPadsY/4] && (xB!=x/4||yB!=y/4||tB!=t/8) && LocBigVoxels[x/4+y/4*fNumberOfPadsX/4+t/8*fNumberOfPadsX/4*fNumberOfPadsY/4]>-1 && LocBigVoxels[xB+yB*fNumberOfPadsX/4+tB*fNumberOfPadsX/4*fNumberOfPadsY/4]>-1){
// Case where both boxes are associated to a track zone, we keep the smallest possible id.
-
+ //cout << 4 << endl;
unsigned int maxid=TMath::Max(LocBigVoxels[x/4+y/4*fNumberOfPadsX/4+t/8*fNumberOfPadsX/4*fNumberOfPadsY/4],LocBigVoxels[xB+yB*fNumberOfPadsX/4+tB*fNumberOfPadsX/4*fNumberOfPadsY/4]);
-
+
unsigned int minid=TMath::Min(LocBigVoxels[x/4+y/4*fNumberOfPadsX/4+t/8*fNumberOfPadsX/4*fNumberOfPadsY/4],LocBigVoxels[xB+yB*fNumberOfPadsX/4+tB*fNumberOfPadsX/4*fNumberOfPadsY/4]);
-
+
unsigned int mysize=TrackZone.at(maxid).GetXPoints().size();
-
+
int element = TrackZone.at(maxid).GetXPoints().at(0) + fNumberOfPadsX/4*TrackZone.at(maxid).GetYPoints().at(0) + fNumberOfPadsX/4*fNumberOfPadsY/4*TrackZone.at(maxid).GetZPoints().at(0);
int Loc=LocBigVoxels[element];
-
+
for(unsigned int i=0;i<mysize;i++){
int element_max = TrackZone.at(maxid).GetXPoints().at(mysize-i-1) + fNumberOfPadsX/4*TrackZone.at(maxid).GetYPoints().at(mysize-i-1) + fNumberOfPadsX/4*fNumberOfPadsY/4*TrackZone.at(maxid).GetZPoints().at(mysize-i-1);
-
+
int element_min = TrackZone.at(minid).GetXPoints().at(0) + fNumberOfPadsX/4*TrackZone.at(minid).GetYPoints().at(0) + fNumberOfPadsX/4*fNumberOfPadsY/4*TrackZone.at(minid).GetZPoints().at(0);
-
+
LocBigVoxels[element_max]=LocBigVoxels[element_min];
-
- TrackZone.at(minid).GetXPoints().push_back(TrackZone.at(maxid).GetXPoints().at(mysize-i-1));
- TrackZone.at(minid).GetYPoints().push_back(TrackZone.at(maxid).GetYPoints().at(mysize-i-1));
- TrackZone.at(minid).GetZPoints().push_back(TrackZone.at(maxid).GetZPoints().at(mysize-i-1));
-
+
+ TrackZone.at(minid).SetPointX(TrackZone.at(maxid).GetXPoints().at(mysize-i-1));
+ TrackZone.at(minid).SetPointY(TrackZone.at(maxid).GetYPoints().at(mysize-i-1));
+ TrackZone.at(minid).SetPointZ(TrackZone.at(maxid).GetZPoints().at(mysize-i-1));
+
TrackZone.at(maxid).GetXPoints().pop_back();
TrackZone.at(maxid).GetYPoints().pop_back();
TrackZone.at(maxid).GetZPoints().pop_back();
+
}
-
+
for(unsigned int i=Loc+1;i<TrackZone.size();i++){
for(unsigned int j=0;j<TrackZone.at(i).GetXPoints().size();j++){
+
element = TrackZone.at(i).GetXPoints().at(j) + fNumberOfPadsX/4*TrackZone.at(i).GetYPoints().at(j) + fNumberOfPadsX/4*fNumberOfPadsY/4*TrackZone.at(i).GetZPoints().at(j);
LocBigVoxels[element]--;
}
@@ -267,7 +279,7 @@ void Cluster::RegroupZone()
x2 = TrackZone.at(j).GetPointX(k2);
y2 = TrackZone.at(j).GetPointY(k2);
z2 = TrackZone.at(j).GetPointZ(k2);
-
+
//if(fabs(x1-x2)+fabs(y1-y2)+fabs(z1-z2)==1 && y1+y2!=fNumberOfPadsY/4-1) it++;
if(fabs(x1-x2)+fabs(y1-y2)+fabs(z1-z2)==1) it++;
}
@@ -313,14 +325,14 @@ void Cluster::FitZone()
NPL::Track myTrack;
for(unsigned int j=0;j<TrackZone.at(i).GetXPoints().size();j++){
int element = TrackZone.at(i).GetPointX(j) + fNumberOfPadsX/4*TrackZone.at(i).GetPointY(j) + fNumberOfPadsX/4*fNumberOfPadsY/4*TrackZone.at(i).GetPointZ(j);
-
+
unsigned int mysize = IDBigVoxels[element].size();
-
+
for(unsigned int k=0;k<mysize;k++){
int x, y, z, q;
q = BigVoxels[element].at(k);
ExtractCoordinate(IDBigVoxels[element].at(k), x, y, z);
-
+
myTrack.SetPointX(x);
myTrack.SetPointY(y);
myTrack.SetPointZ(z);
@@ -329,7 +341,7 @@ void Cluster::FitZone()
}
if(myTrack.GetXPoints().size()>0){
Fit3D(myTrack);
-
+
TrackZone[i].SetXm(myTrack.GetXm());
TrackZone[i].SetYm(myTrack.GetYm());
TrackZone[i].SetZm(myTrack.GetZm());
@@ -348,14 +360,14 @@ void Cluster::FitZone()
void Cluster::RegroupTrack()
{
ReorderTracks(TrackZone);
-
+
unsigned int TrackZoneSize = TrackZone.size();
int ite=0;
-
+
for(unsigned int id=0; id<TrackZoneSize; id++){
if(TrackZone[id].GetXPoints().size()>0){
NPL::Track newTrack;
-
+
int TrackZoneNumberOfPoints=0;
int xT, xTZ;
int yT, yTZ;
@@ -367,7 +379,7 @@ void Cluster::RegroupTrack()
int element2 = TrackZone.at(i).GetPointX(j) + fNumberOfPadsX/4*TrackZone.at(i).GetPointY(j) + fNumberOfPadsX/4*fNumberOfPadsY/4*TrackZone.at(i).GetPointZ(j);
ExtractCoordinate(IDBigVoxels[element2].at(0),xTZ,yTZ,zTZ);
-
+
//We verify that we are on the same side.
if((yTZ-fNumberOfPadsY/2.)*(yT-fNumberOfPadsY/2.)>=0){
if(i==id) TrackZoneNumberOfPoints+= IDBigVoxels[element2].size();
@@ -376,12 +388,12 @@ void Cluster::RegroupTrack()
ExtractCoordinate(IDBigVoxels[element2].at(k),xTZ,yTZ,zTZ);
if(yTZ-fNumberOfPadsY/2.-1) yTZ=fNumberOfPadsY/2.;
if( (i==id && BigVoxels[element2].at(k)>4050) || (((i!=id && (xTZ>=100 || fabs(yTZ-fNumberOfPadsY/2)>=6)) || i==id) && (yTZ-fNumberOfPadsY/2.)*(yT-fNumberOfPadsY/2.+0.1)>=0 && IDBigVoxels[element2].at(k)>0 && Distance3D(TrackZone.at(id),IDBigVoxels[element2].at(k))<fClusterDistance) ){
-
+
ite++;
-
+
int x_track, y_track, z_track;
int q_track = BigVoxels[element2].at(k);
-
+
ExtractCoordinate(IDBigVoxels[element2].at(k),x_track,y_track,z_track);
newTrack.SetPointX(x_track);
newTrack.SetPointY(y_track);
@@ -397,7 +409,7 @@ void Cluster::RegroupTrack()
if((float)ite>=TrackZoneNumberOfPoints*0.9){
vTrack.push_back(newTrack);
}
-
+
newTrack.Clear();
}
}
@@ -413,7 +425,7 @@ void Cluster::FitTrack()
}
else to_erase.push_back(i);
}
-
+
for(int i=to_erase.size()-1; i>=0; i--) vTrack.erase(vTrack.begin()+to_erase[i]);
}
@@ -424,17 +436,17 @@ double Cluster::Distance3D(NPL::Track aTrack, int ID)
double xf, yf, zf;
double distance;
ExtractCoordinate(ID,x,y,z);
-
+
TVector3 V1;
TVector3 V2;
-
+
V1.SetXYZ(aTrack.GetXm()-aTrack.GetXh(), aTrack.GetYm()-aTrack.GetYh(), aTrack.GetZm()-aTrack.GetZh());
V2.SetXYZ(aTrack.GetXm()-x, aTrack.GetYm()-y, aTrack.GetZm()-z);
TVector3 V2CV1 = V2.Cross(V1);
-
+
distance = V2CV1.Mag()/V1.Mag();
-
+
return distance;
}
@@ -478,12 +490,12 @@ void Cluster::RemoveTrack(vector<NPL::Track> &Ephem, int id)
void Cluster::Init(vector<int> v1, vector<int> v2, vector<double> v3, vector<double> v4)
{
Reset();
-
+
vX = v1;
vY = v2;
vZ = v3;
vQ = v4;
-
+
fOriginalCloudSize = vX.size();
double TotalCharge=0;
for(unsigned int i=0; i< vQ.size(); i++){
@@ -502,7 +514,7 @@ void Cluster::Reset()
vZ.clear();
vQ.clear();
TrackZone.clear();
-
+
for(int i=0; i<BigVoxelSize; i++){
BigVoxels[i].clear();
IDBigVoxels[i].clear();
@@ -521,16 +533,16 @@ void Cluster::Clear()
double Cluster::Fit3D(NPL::Track &aTrack)
{
// adapted from: http://fr.scribd.com/doc/31477970/Regressions-et-trajectoires-3D
-
+
vector<int> X;
vector<int> Y;
vector<int> Z;
-
+
X = aTrack.GetXPoints();
Y = aTrack.GetYPoints();
Z = aTrack.GetZPoints();
-
-
+
+
int R, C;
double Q;
double Xm,Ym,Zm;
@@ -542,10 +554,10 @@ double Cluster::Fit3D(NPL::Track &aTrack)
double c0,c1,c2;
double p,q,r,dm2;
double rho,phi;
-
+
Q=Xm=Ym=Zm=0.;
Sxx=Syy=Szz=Sxy=Sxz=Syz=0.;
-
+
for (unsigned int i=0; i<X.size(); i++)
{
Q++;
@@ -559,7 +571,7 @@ double Cluster::Fit3D(NPL::Track &aTrack)
Sxz+=X[i]*Z[i];
Syz+=Y[i]*Z[i];
}
-
+
Xm/=Q;
Ym/=Q;
Zm/=Q;
@@ -575,26 +587,26 @@ double Cluster::Fit3D(NPL::Track &aTrack)
Sxy-=(Xm*Ym);
Sxz-=(Xm*Zm);
Syz-=(Ym*Zm);
-
+
theta=0.5*atan((2.*Sxy)/(Sxx-Syy));
-
+
K11=(Syy+Szz)*pow(cos(theta),2)+(Sxx+Szz)*pow(sin(theta),2)-2.*Sxy*cos(theta)*sin(theta);
K22=(Syy+Szz)*pow(sin(theta),2)+(Sxx+Szz)*pow(cos(theta),2)+2.*Sxy*cos(theta)*sin(theta);
K12=-Sxy*(pow(cos(theta),2)-pow(sin(theta),2))+(Sxx-Syy)*cos(theta)*sin(theta);
K10=Sxz*cos(theta)+Syz*sin(theta);
K01=-Sxz*sin(theta)+Syz*cos(theta);
K00=Sxx+Syy;
-
+
c2=-K00-K11-K22;
c1=K00*K11+K00*K22+K11*K22-K01*K01-K10*K10;
c0=K01*K01*K11+K10*K10*K22-K00*K11*K22;
-
-
+
+
p=c1-pow(c2,2)/3.;
q=2.*pow(c2,3)/27.-c1*c2/3.+c0;
r=pow(q/2.,2)+pow(p,3)/27.;
-
-
+
+
if(r>0) dm2=-c2/3.+pow(-q/2.+sqrt(r),1./3.)+pow(-q/2.-sqrt(r),1./3.);
if(r<0)
{
@@ -602,14 +614,14 @@ double Cluster::Fit3D(NPL::Track &aTrack)
phi=acos(-q/(2.*rho));
dm2=min(-c2/3.+2.*pow(rho,1./3.)*cos(phi/3.),min(-c2/3.+2.*pow(rho,1./3.)*cos((phi+2.*TMath::Pi())/3.),-c2/3.+2.*pow(rho,1./3.)*cos((phi+4.*TMath::Pi())/3.)));
}
-
+
a=-K10*cos(theta)/(K11-dm2)+K01*sin(theta)/(K22-dm2);
b=-K10*sin(theta)/(K11-dm2)-K01*cos(theta)/(K22-dm2);
-
+
Xh=((1.+b*b)*Xm-a*b*Ym+a*Zm)/(1.+a*a+b*b);
Yh=((1.+a*a)*Ym-a*b*Xm+b*Zm)/(1.+a*a+b*b);
Zh=((a*a+b*b)*Zm+a*Xm+b*Ym)/(1.+a*a+b*b);
-
+
aTrack.SetXm(Xm);
aTrack.SetYm(Ym);
aTrack.SetZm(Zm);
@@ -617,36 +629,7 @@ double Cluster::Fit3D(NPL::Track &aTrack)
aTrack.SetYh(Yh);
aTrack.SetZh(Zh);
aTrack.SetSlopesAndOffsets();
-
-
- return(fabs(dm2/Q));
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+ return(fabs(dm2/Q));
+}