diff --git a/Examples/Example4/Analysis.cxx b/Examples/Example4/Analysis.cxx index a89dbc82f8a69343a42001e7f72ea0a09f08ccb4..a9a9a1ab0840c11b788bc4100169b431e09696e6 100644 --- a/Examples/Example4/Analysis.cxx +++ b/Examples/Example4/Analysis.cxx @@ -78,7 +78,6 @@ void Analysis::TreatEvent(){ TVector3 vX = TVector3(1,0,0); TVector3 aTrack, vB; - if(TrackMult>1){ vTrack = Actar->GetTracks(); double scalarproduct=0; @@ -106,7 +105,8 @@ void Analysis::TreatEvent(){ double ZBeamPoint = vTrack[BeamTrack].GetZh(); TVector3 vBeamPos = TVector3(XBeamPoint,YBeamPoint,ZBeamPoint); - vB = TVector3(XBeam*PadSizeX, YBeam*PadSizeY,ZBeam*DriftVelocity); + //vB = TVector3(XBeam*PadSizeX, YBeam*PadSizeY,ZBeam*DriftVelocity); + vB = TVector3(XBeam*PadSizeX, YBeam*PadSizeY,ZBeam); BeamAngle = (vX.Angle(vB))*180/TMath::Pi(); for(unsigned int i=0; i<TrackMult; i++){ @@ -117,9 +117,11 @@ void Analysis::TreatEvent(){ double vertex_x = vTrack[i].GetVertexPostion(vBeam,vBeamPos).X()*PadSizeX; double vertex_y = vTrack[i].GetVertexPostion(vBeam,vBeamPos).Y()*PadSizeY; - double vertex_z = vTrack[i].GetVertexPostion(vBeam,vBeamPos).Z()*DriftVelocity; + //double vertex_z = vTrack[i].GetVertexPostion(vBeam,vBeamPos).Z()*DriftVelocity; + double vertex_z = vTrack[i].GetVertexPostion(vBeam,vBeamPos).Z(); - aTrack = TVector3(Xdir*PadSizeX, Ydir*PadSizeY, Zdir*DriftVelocity); + //aTrack = TVector3(Xdir*PadSizeX, Ydir*PadSizeY, Zdir*DriftVelocity); + aTrack = TVector3(Xdir*PadSizeX, Ydir*PadSizeY, Zdir); double angle = vX.Angle(aTrack)*180/TMath::Pi(); //double angle = vB.Angle(aTrack)*180/TMath::Pi(); if(angle>90) angle = 180-angle; @@ -130,8 +132,10 @@ void Analysis::TreatEvent(){ double y2 = vTrack[i].GetYh()*PadSizeY-0.5*NumberOfPadsY*PadSizeY; //double z1 = -(vTrack[i].GetZm()-256)*DriftVelocity; //double z2 = -(vTrack[i].GetZh()-256)*DriftVelocity; - double z1 = vTrack[i].GetZm()*DriftVelocity; - double z2 = vTrack[i].GetZh()*DriftVelocity; + //double z1 = vTrack[i].GetZm()*DriftVelocity; + double z1 = vTrack[i].GetZm(); + //double z2 = vTrack[i].GetZh()*DriftVelocity; + double z2 = vTrack[i].GetZh(); if(vertex_x>0 && vertex_x<256){ diff --git a/Examples/Example4/configs/ConfigActar.dat b/Examples/Example4/configs/ConfigActar.dat index e3f69115de47822e40dc83ffe7ddc01bae6461ea..b1cbe69c7d64792ea68662e16039e2a124f77efc 100644 --- a/Examples/Example4/configs/ConfigActar.dat +++ b/Examples/Example4/configs/ConfigActar.dat @@ -1,7 +1,7 @@ ConfigActar -RecoRansac= 0 -RecoCluster= 1 -RecoVisu= 0 +RecoRansac= 1 +RecoCluster= 0 +RecoVisu= 1 HIT_THRESHOLD= 2 Q_THRESHOLD= 0 T_THRESHOLD= 0 @@ -17,4 +17,4 @@ Gas= iC4H10 %Gas= D2 %Pressure in mbar Pressure= 100 -DriftVelocity= 2.66 +DriftVelocity= 0.8e-2 diff --git a/Examples/Example4/macro/CheckReco.C b/Examples/Example4/macro/CheckReco.C index 59a59b28bf5e6df82e9764854baca1285089b81a..1ac16248ef717c556517ff3ee92441f9f6c3f520 100644 --- a/Examples/Example4/macro/CheckReco.C +++ b/Examples/Example4/macro/CheckReco.C @@ -17,7 +17,7 @@ void LoadCut() void LoadChain() { chain = new TChain("PhysicsTree"); - chain->Add("../../Outputs/Analysis/Example4.root"); + chain->Add("../../../Outputs/Analysis/Example4.root"); } //////////////////////////////////////////////// diff --git a/Examples/Example4/macro/ShowResult.C b/Examples/Example4/macro/ShowResult.C index 79904dca34220881593a386cae17111e504015b4..ec564894c7642904094b1662d03c2534b0dffb94 100644 --- a/Examples/Example4/macro/ShowResult.C +++ b/Examples/Example4/macro/ShowResult.C @@ -18,7 +18,7 @@ void LoadCut() void LoadChain() { chain = new TChain("PhysicsTree"); - chain->Add("../../Outputs/Analysis/Example4.root"); + chain->Add("../../../Outputs/Analysis/Example4.root"); } //////////////////////////////////////////////// diff --git a/NPLib/Detectors/Actar/TActarPhysics.cxx b/NPLib/Detectors/Actar/TActarPhysics.cxx index ddf0d8a1e7767834a121946162b41b520f5332bf..1d83d2156d6516e45da21d1929464c06467467e8 100644 --- a/NPLib/Detectors/Actar/TActarPhysics.cxx +++ b/NPLib/Detectors/Actar/TActarPhysics.cxx @@ -45,55 +45,55 @@ using namespace std; ClassImp(TActarPhysics) -/////////////////////////////////////////////////////////////////////////// + /////////////////////////////////////////////////////////////////////////// TActarPhysics::TActarPhysics() -: m_EventData(new TActarData), -m_EventReduced(new MEventReduced), -m_PreTreatedData(new TActarData), -m_EventPhysics(this), - -m_Spectra(0), -fRecoRansac(1), -fRecoCluster(0), -fRecoVisu(0), -fHitThreshold(20), -fQ_Threshold(0), -fT_Threshold(0), -fXBeamMin(0), -fXBeamMax(128), -fYBeamMin(60), -fYBeamMax(67), -fNumberOfPadsX(128), -fNumberOfPadsY(128), -fPadSizeX(2), -fPadSizeY(2), -fDriftVelocity(40), -fPressure(100), -fGas("iC4H10"), -m_NumberOfPadSilicon(20), -m_NumberOfDetectors(0) { -} + : m_EventData(new TActarData), + m_EventReduced(new MEventReduced), + m_PreTreatedData(new TActarData), + m_EventPhysics(this), + + m_Spectra(0), + fRecoRansac(1), + fRecoCluster(0), + fRecoVisu(0), + fHitThreshold(20), + fQ_Threshold(0), + fT_Threshold(0), + fXBeamMin(0), + fXBeamMax(128), + fYBeamMin(60), + fYBeamMax(67), + fNumberOfPadsX(128), + fNumberOfPadsY(128), + fPadSizeX(2), + fPadSizeY(2), + fDriftVelocity(40), + fPressure(100), + fGas("iC4H10"), + m_NumberOfPadSilicon(20), + m_NumberOfDetectors(0) { + } /////////////////////////////////////////////////////////////////////////// /// A usefull method to bundle all operation to add a detector void TActarPhysics::AddDetector(TVector3 , string ){ - // In That simple case nothing is done - // Typically for more complex detector one would calculate the relevant - // positions (stripped silicon) or angles (gamma array) - m_NumberOfDetectors++; + // In That simple case nothing is done + // Typically for more complex detector one would calculate the relevant + // positions (stripped silicon) or angles (gamma array) + m_NumberOfDetectors++; } /////////////////////////////////////////////////////////////////////////// void TActarPhysics::AddDetector(double R, double Theta, double Phi, string shape){ - // 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); + // 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 TActarPhysics::BuildSimplePhysicalEvent() { - BuildPhysicalEvent(); + BuildPhysicalEvent(); } @@ -101,570 +101,570 @@ void TActarPhysics::BuildSimplePhysicalEvent() { /////////////////////////////////////////////////////////////////////////// void TActarPhysics::BuildPhysicalEvent() { - PreTreat(); + PreTreat(); - //unsigned int mysize = m_PreTreatedData->GetPadMult(); + //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)); + /*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(); + if(fRecoRansac && PadX.size()>fHitThreshold){ + //m_Ransac->Init(PadX, PadY, PadZ, PadCharge); + m_Ransac->Init(PadX, PadY, PadZ, PadCharge); + m_Track = m_Ransac->SimpleRansac(); + } + + else if(fRecoCluster){ + if(PadX.size()>fHitThreshold){ + m_Cluster->Init(PadX, PadY, PadZ, PadCharge); + m_Cluster->FindTracks(); } - 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(); - } + if(BeamPadX.size()>fHitThreshold){ + m_Cluster->Init(BeamPadX, BeamPadY, BeamPadZ, BeamPadCharge); + m_Cluster->FindTracks(); + } - m_Track = m_Cluster->GetTracks(); + m_Track = m_Cluster->GetTracks(); - m_Cluster->Clear(); - } + m_Cluster->Clear(); + } - TrackMult = m_Track.size(); + TrackMult = m_Track.size(); } /////////////////////////////////////////////////////////////////////////// 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++){ - if(m_EventReduced->CoboAsad[it].peakheight[hh]>fQ_Threshold && m_EventReduced->CoboAsad[it].peaktime[hh]>fT_Threshold){ - if(GoodHit(TABLE[4][where],TABLE[5][where])){ - //if(Hit[TABLE[4][where]][TABLE[5][where]]<2){ - if(fRecoCluster){ - if(IsBeamZone(TABLE[4][where],TABLE[5][where])){ - BeamPadCharge.push_back(m_EventReduced->CoboAsad[it].peakheight[hh]); - BeamPadX.push_back(TABLE[4][where]); - BeamPadY.push_back(TABLE[5][where]); - BeamPadZ.push_back(m_EventReduced->CoboAsad[it].peaktime[hh]); - } - else{ - PadCharge.push_back(m_EventReduced->CoboAsad[it].peakheight[hh]); - PadX.push_back(TABLE[4][where]); - 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]); - PadX.push_back(TABLE[4][where]); - PadY.push_back(TABLE[5][where]); - PadZ.push_back(m_EventReduced->CoboAsad[it].peaktime[hh]); - } - //} - } - } + // 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(); + //cout << vector_size << endl; + for(unsigned int hh=0; hh<vector_size; hh++){ + if(m_EventReduced->CoboAsad[it].peakheight[hh]>fQ_Threshold && m_EventReduced->CoboAsad[it].peaktime[hh]>fT_Threshold){ + if(GoodHit(TABLE[4][where],TABLE[5][where])){ + //if(Hit[TABLE[4][where]][TABLE[5][where]]<2){ + if(fRecoCluster){ + if(IsBeamZone(TABLE[4][where],TABLE[5][where])){ + BeamPadCharge.push_back(m_EventReduced->CoboAsad[it].peakheight[hh]); + BeamPadX.push_back(TABLE[4][where]); + BeamPadY.push_back(TABLE[5][where]); + BeamPadZ.push_back(int(m_EventReduced->CoboAsad[it].peaktime[hh]*fDriftVelocity)); + } + else{ + PadCharge.push_back(m_EventReduced->CoboAsad[it].peakheight[hh]); + PadX.push_back(TABLE[4][where]); + PadY.push_back(TABLE[5][where]); + PadZ.push_back(int(m_EventReduced->CoboAsad[it].peaktime[hh]*fDriftVelocity)); + } + } - } - else if(co==31){ - unsigned int vector_size = m_EventReduced->CoboAsad[it].peakheight.size(); - for(unsigned int hit=0;hit<vector_size;hit++){ - 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]); - } + else if(fRecoRansac){ + PadCharge.push_back(m_EventReduced->CoboAsad[it].peakheight[hh]); + PadX.push_back(TABLE[4][where]); + PadY.push_back(TABLE[5][where]); + PadZ.push_back(int(m_EventReduced->CoboAsad[it].peaktime[hh]*fDriftVelocity)); } + //} + } + } + } + } + else if(co==31){ + unsigned int vector_size = m_EventReduced->CoboAsad[it].peakheight.size(); + for(unsigned int hit=0;hit<vector_size;hit++){ + 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]); + } + } } + } } /////////////////////////////////////////////////////////////////////////// bool TActarPhysics::IsBeamZone(int X, int Y) { - bool isBeam=false; - if( (X>=fXBeamMin && X<=fXBeamMax) && (Y>=fYBeamMin && Y<=fYBeamMax) ){ - isBeam=true; - } + bool isBeam=false; + if( (X>=fXBeamMin && X<=fXBeamMax) && (Y>=fYBeamMin && Y<=fYBeamMax) ){ + isBeam=true; + } - return isBeam; + return isBeam; } /////////////////////////////////////////////////////////////////////////// bool TActarPhysics::GoodHit(int iX, int iY) { - bool bHit = true; - if(Hit[iX][iY]<2){ - if(Hit[iX+1][iY]>0 || Hit[iX+1][iY-1]>0 || Hit[iX+1][iY+1]>0){ - if(Hit[iX+2][iY]>0 || Hit[iX+2][iY-1]>0 || Hit[iX+2][iY+1]>0){ - bHit = true; - } - } + bool bHit = true; + if(Hit[iX][iY]<2){ + if(Hit[iX+1][iY]>0 || Hit[iX+1][iY-1]>0 || Hit[iX+1][iY+1]>0){ + if(Hit[iX+2][iY]>0 || Hit[iX+2][iY-1]>0 || Hit[iX+2][iY+1]>0){ + bHit = true; + } } + } - return bHit; + return bHit; } /////////////////////////////////////////////////////////////////////////// void TActarPhysics::CleanPads() { - 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++){ - if(m_EventReduced->CoboAsad[it].peakheight[hh]>fQ_Threshold && m_EventReduced->CoboAsad[it].peaktime[hh]>fT_Threshold){ - Hit[TABLE[4][where]][TABLE[5][where]] += 1; - /* - m_PreTreatedData->SetCharge(m_EventReduced->CoboAsad[it].peakheight[hh]); - m_PreTreatedData->SetPadX(TABLE[4][where]); - m_PreTreatedData->SetPadY(TABLE[5][where]); - m_PreTreatedData->SetPadZ(m_EventReduced->CoboAsad[it].peaktime[hh]);*/ - } - } + 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++){ + if(m_EventReduced->CoboAsad[it].peakheight[hh]>fQ_Threshold && m_EventReduced->CoboAsad[it].peaktime[hh]>fT_Threshold){ + Hit[TABLE[4][where]][TABLE[5][where]] += 1; + /* + m_PreTreatedData->SetCharge(m_EventReduced->CoboAsad[it].peakheight[hh]); + m_PreTreatedData->SetPadX(TABLE[4][where]); + m_PreTreatedData->SetPadY(TABLE[5][where]); + m_PreTreatedData->SetPadZ(m_EventReduced->CoboAsad[it].peaktime[hh]);*/ } + } } + } } /////////////////////////////////////////////////////////////////////////// void TActarPhysics::ReadAnalysisConfig() { - // VXI ACTION FILE // - string VXI_FileName = "./configs/ACTION_Si_config.dat"; - ifstream VXIConfigFile; - VXIConfigFile.open(VXI_FileName.c_str()); - if(!VXIConfigFile.is_open()){ - cout << "No VXI ACTION FILE ./configs/ACTION_Si_config.dat found!" << endl; - return; - } - else{ - cout << "/// Using VXI ACTION FILE: " << VXI_FileName << " ///" << endl; - string token; - int vxi_param, si_nbr; - for(int i=0; i<20; i++){ - VXIConfigFile >> token >> vxi_param >> si_nbr; - Si_map[vxi_param] = si_nbr+1; - } + // VXI ACTION FILE // + string VXI_FileName = "./configs/ACTION_Si_config.dat"; + ifstream VXIConfigFile; + VXIConfigFile.open(VXI_FileName.c_str()); + if(!VXIConfigFile.is_open()){ + cout << "No VXI ACTION FILE ./configs/ACTION_Si_config.dat found!" << endl; + return; + } + else{ + cout << "/// Using VXI ACTION FILE: " << VXI_FileName << " ///" << endl; + string token; + int vxi_param, si_nbr; + for(int i=0; i<20; i++){ + VXIConfigFile >> token >> vxi_param >> si_nbr; + Si_map[vxi_param] = si_nbr+1; } - VXIConfigFile.close(); - - // Lookup table // - bool ReadingLookupTable = false; - string LT_FileName = "./configs/LT.dat"; - ifstream LTConfigFile; - LTConfigFile.open(LT_FileName.c_str()); - if(!LTConfigFile.is_open()){ - cout << "No Lookup Table in ./configs/LT.dat found!" << endl; - return; + } + VXIConfigFile.close(); + + // Lookup table // + bool ReadingLookupTable = false; + string LT_FileName = "./configs/LT.dat"; + ifstream LTConfigFile; + LTConfigFile.open(LT_FileName.c_str()); + if(!LTConfigFile.is_open()){ + cout << "No Lookup Table in ./configs/LT.dat found!" << endl; + return; + } + else{ + cout << "/// Using LookupTable from: " << LT_FileName << " ///" << endl; + for(int i=0;i<NumberOfCobo*NumberOfASAD*NumberOfAGET*NumberOfChannel;i++){ + LTConfigFile >> TABLE[0][i] >> TABLE[1][i] >> TABLE[2][i] >> TABLE[3][i] >> TABLE[4][i] >> TABLE[5][i]; } - else{ - cout << "/// Using LookupTable from: " << LT_FileName << " ///" << endl; - for(int i=0;i<NumberOfCobo*NumberOfASAD*NumberOfAGET*NumberOfChannel;i++){ - LTConfigFile >> TABLE[0][i] >> TABLE[1][i] >> TABLE[2][i] >> TABLE[3][i] >> TABLE[4][i] >> TABLE[5][i]; - } - 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; + 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 %%%"); + asciiConfig->Append(FileName.c_str()); + asciiConfig->AppendLine(""); + // read analysis config file + string LineBuffer,DataBuffer,whatToDo; + 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; } - cout << "/// Loading user parameter for Analysis from ConfigActar.dat ///" << endl; - - // Save it in a TAsciiFile - TAsciiFile* asciiConfig = RootOutput::getInstance()->GetAsciiFileAnalysisConfig(); - asciiConfig->AppendLine("%%% ConfigActar.dat %%%"); - asciiConfig->Append(FileName.c_str()); - asciiConfig->AppendLine(""); - // read analysis config file - string LineBuffer,DataBuffer,whatToDo; - 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; + // 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; - } + // Search for comment symbol (%) + if (whatToDo.compare(0, 1, "%") == 0) { + AnalysisConfigFile.ignore(numeric_limits<streamsize>::max(), '\n' ); + } - 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,11,"RecoRansac=") == 0) { + AnalysisConfigFile >> DataBuffer; + fRecoRansac = atoi(DataBuffer.c_str()); + cout << "/// Reco using Ransac= " << " " << fRecoRansac << " ///" << 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,12,"RecoCluster=") == 0) { + AnalysisConfigFile >> DataBuffer; + fRecoCluster = atoi(DataBuffer.c_str()); + cout << "/// Reco using Cluster= " << " " << fRecoCluster << " ///" << 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,9,"RecoVisu=") == 0) { + AnalysisConfigFile >> DataBuffer; + fRecoVisu = atoi(DataBuffer.c_str()); + cout << "/// Visualisation= " << " " << fRecoVisu << " ///" << 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,14,"HIT_THRESHOLD=") == 0) { + AnalysisConfigFile >> DataBuffer; + fHitThreshold = atoi(DataBuffer.c_str()); + cout << "/// Hit Threshold= " << " " << fHitThreshold << " ///" << endl; + } - else if (whatToDo.compare(0,12,"T_THRESHOLD=") == 0) { - AnalysisConfigFile >> DataBuffer; - fT_Threshold = atof(DataBuffer.c_str()); - cout << "/// T Threshold= " << " " << fT_Threshold << " ///" << endl; - } - else if(whatToDo.compare(0,14,"NumberOfPadsX=")==0){ - AnalysisConfigFile >> DataBuffer; - fNumberOfPadsX = atoi(DataBuffer.c_str()); - //check_padsX=true; - cout << "/// Number Of Pads X= " << fNumberOfPadsX << " ///" << 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,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,12,"T_THRESHOLD=") == 0) { + AnalysisConfigFile >> DataBuffer; + fT_Threshold = atof(DataBuffer.c_str()); + cout << "/// T Threshold= " << " " << fT_Threshold << " ///" << endl; + } + else if(whatToDo.compare(0,14,"NumberOfPadsX=")==0){ + AnalysisConfigFile >> DataBuffer; + fNumberOfPadsX = atoi(DataBuffer.c_str()); + //check_padsX=true; + cout << "/// Number Of Pads X= " << fNumberOfPadsX << " ///" << 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,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,"PadSizeY=")==0){ - AnalysisConfigFile >> DataBuffer; - fPadSizeY = atof(DataBuffer.c_str()); - //check_sizeY=true; - cout << "/// Pad Size Y= " << fPadSizeY << " ///" << 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,"XBeamMin=")==0){ - AnalysisConfigFile >> DataBuffer; - fXBeamMin = atof(DataBuffer.c_str()); - cout << "/// X Beam Min= " << fXBeamMin << " ///" << 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,"XBeamMax=")==0){ - AnalysisConfigFile >> DataBuffer; - fXBeamMax = atof(DataBuffer.c_str()); - cout << "/// X Beam Max= " << fXBeamMax << " ///" << 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,"YBeamMin=")==0){ - AnalysisConfigFile >> DataBuffer; - fYBeamMin = atof(DataBuffer.c_str()); - cout << "/// Y Beam Min= " << fYBeamMin << " ///" << 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,"YBeamMax=")==0){ - AnalysisConfigFile >> DataBuffer; - fYBeamMax = atof(DataBuffer.c_str()); - cout << "/// Y Beam Max= " << fYBeamMax << " ///" << 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,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,"YBeamMax=")==0){ + AnalysisConfigFile >> DataBuffer; + fYBeamMax = atof(DataBuffer.c_str()); + cout << "/// Y Beam Max= " << fYBeamMax << " ///" << 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,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,14,"DriftVelocity=")==0){ - AnalysisConfigFile >> DataBuffer; - fDriftVelocity = atof(DataBuffer.c_str()); - //check_driftvelocity=true; - cout << "/// Drift Velocity= " << fDriftVelocity << " ///" << 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,4,"Gas=")==0){ - AnalysisConfigFile >> DataBuffer; - fGas = DataBuffer.c_str(); - //check_gas=true; - cout << "/// Gas Type= " << fGas << " ///" << 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 { - ReadingStatus = false; - } - } - } + else if(whatToDo.compare(0,4,"Gas=")==0){ + AnalysisConfigFile >> DataBuffer; + fGas = DataBuffer.c_str(); + //check_gas=true; + cout << "/// Gas Type= " << fGas << " ///" << endl; + } - if(fRecoRansac){ - m_Ransac = new NPL::Ransac(fNumberOfPadsX,fNumberOfPadsY,fRecoVisu); - } - else if(fRecoCluster){ - m_Cluster = new NPL::Cluster(fNumberOfPadsX,fNumberOfPadsY,fRecoVisu); + else { + ReadingStatus = false; + } } + } + + if(fRecoRansac){ + m_Ransac = new NPL::Ransac(fNumberOfPadsX,fNumberOfPadsY,fRecoVisu); + } + else if(fRecoCluster){ + m_Cluster = new NPL::Cluster(fNumberOfPadsX,fNumberOfPadsY,fRecoVisu); + } } /////////////////////////////////////////////////////////////////////////// void TActarPhysics::SetRansacParameter(string filename){ - if(fRecoRansac){ - m_Ransac->ReadParameterValue(filename); - } + if(fRecoRansac){ + m_Ransac->ReadParameterValue(filename); + } } /////////////////////////////////////////////////////////////////////////// void TActarPhysics::SetClusterParameter(string filename){ - if(fRecoCluster){ - m_Cluster->ReadParameterValue(filename); - } + if(fRecoCluster){ + m_Cluster->ReadParameterValue(filename); + } } /////////////////////////////////////////////////////////////////////////// void TActarPhysics::Clear() { - BeamPadX.clear(); - BeamPadY.clear(); - BeamPadZ.clear(); - BeamPadCharge.clear(); - - PadX.clear(); - PadY.clear(); - PadZ.clear(); - PadCharge.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; - } + BeamPadX.clear(); + BeamPadY.clear(); + BeamPadZ.clear(); + BeamPadCharge.clear(); + + PadX.clear(); + PadY.clear(); + PadZ.clear(); + PadCharge.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; } + } } /////////////////////////////////////////////////////////////////////////// void TActarPhysics::ReadConfiguration(NPL::InputParser parser) { - vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Actar"); - if(NPOptionManager::getInstance()->GetVerboseLevel()) - cout << "//// " << blocks.size() << " detectors found " << endl; + 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"}; + 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; + 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); - } - else if(blocks[i]->HasTokenList(sphe)){ - if(NPOptionManager::getInstance()->GetVerboseLevel()) - cout << endl << "//// Actar " << 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); - } + 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 << "//// Actar " << 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 TActarPhysics::InitSpectra() { - m_Spectra = new TActarSpectra(m_NumberOfDetectors); + m_Spectra = new TActarSpectra(m_NumberOfDetectors); } /////////////////////////////////////////////////////////////////////////// void TActarPhysics::FillSpectra() { - m_Spectra -> FillRawSpectra(m_EventData); - m_Spectra -> FillPreTreatedSpectra(m_PreTreatedData); - m_Spectra -> FillPhysicsSpectra(m_EventPhysics); + m_Spectra -> FillRawSpectra(m_EventData); + m_Spectra -> FillPreTreatedSpectra(m_PreTreatedData); + m_Spectra -> FillPhysicsSpectra(m_EventPhysics); } /////////////////////////////////////////////////////////////////////////// void TActarPhysics::CheckSpectra() { - m_Spectra->CheckSpectra(); + m_Spectra->CheckSpectra(); } /////////////////////////////////////////////////////////////////////////// void TActarPhysics::ClearSpectra() { - // To be done + // To be done } /////////////////////////////////////////////////////////////////////////// map< string , TH1*> TActarPhysics::GetSpectra() { - if(m_Spectra) - return m_Spectra->GetMapHisto(); - else{ - map< string , TH1*> empty; - return empty; - } + if(m_Spectra) + return m_Spectra->GetMapHisto(); + else{ + map< string , TH1*> empty; + return empty; + } } //////////////////////////////////////////////////////////////////////////////// vector<TCanvas*> TActarPhysics::GetCanvas() { - //if(m_Spectra) - //return m_Spectra->GetCanvas(); - //else{ - vector<TCanvas*> empty; - return empty; - //} + //if(m_Spectra) + //return m_Spectra->GetCanvas(); + //else{ + vector<TCanvas*> empty; + return empty; + //} } /////////////////////////////////////////////////////////////////////////// void TActarPhysics::WriteSpectra() { - m_Spectra->WriteSpectra(); + m_Spectra->WriteSpectra(); } /////////////////////////////////////////////////////////////////////////// void TActarPhysics::AddParameterToCalibrationManager() { - CalibrationManager* Cal = CalibrationManager::getInstance(); - for (int i = 0; i < m_NumberOfDetectors; ++i) { - 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"); - } + CalibrationManager* Cal = CalibrationManager::getInstance(); + for (int i = 0; i < m_NumberOfDetectors; ++i) { + 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"); + } } /////////////////////////////////////////////////////////////////////////// /*void TActarPhysics::InitializeRootInputRaw() { - TChain* inputChain = RootInput::getInstance()->GetChain(); - inputChain->SetBranchStatus("Actar", true ); - inputChain->SetBranchAddress("Actar", &m_EventData ); -}*/ + TChain* inputChain = RootInput::getInstance()->GetChain(); + inputChain->SetBranchStatus("Actar", true ); + inputChain->SetBranchAddress("Actar", &m_EventData ); + }*/ /////////////////////////////////////////////////////////////////////////// void TActarPhysics::InitializeRootInputRaw() { - TChain* inputChain = RootInput::getInstance()->GetChain(); - inputChain->SetBranchStatus("data", true ); - inputChain->SetBranchStatus("ACTAR_TTree", true ); - inputChain->SetBranchAddress("data", &m_EventReduced ); + TChain* inputChain = RootInput::getInstance()->GetChain(); + inputChain->SetBranchStatus("data", true ); + inputChain->SetBranchStatus("ACTAR_TTree", true ); + inputChain->SetBranchAddress("data", &m_EventReduced ); - fInputTreeName = inputChain->GetTree()->GetName(); + fInputTreeName = inputChain->GetTree()->GetName(); } /////////////////////////////////////////////////////////////////////////// void TActarPhysics::InitializeRootInputPhysics() { - TChain* inputChain = RootInput::getInstance()->GetChain(); - inputChain->SetBranchAddress("Actar", &m_EventPhysics); + TChain* inputChain = RootInput::getInstance()->GetChain(); + inputChain->SetBranchAddress("Actar", &m_EventPhysics); } /////////////////////////////////////////////////////////////////////////// void TActarPhysics::InitializeRootOutput() { - TTree* outputTree = RootOutput::getInstance()->GetTree(); - outputTree->Branch("Actar", "TActarPhysics", &m_EventPhysics); + TTree* outputTree = RootOutput::getInstance()->GetTree(); + outputTree->Branch("Actar", "TActarPhysics", &m_EventPhysics); } @@ -673,7 +673,7 @@ void TActarPhysics::InitializeRootOutput() { // Construct Method to be pass to the DetectorFactory // //////////////////////////////////////////////////////////////////////////////// NPL::VDetector* TActarPhysics::Construct() { - return (NPL::VDetector*) new TActarPhysics(); + return (NPL::VDetector*) new TActarPhysics(); } @@ -682,13 +682,13 @@ NPL::VDetector* TActarPhysics::Construct() { // Registering the construct method to the factory // //////////////////////////////////////////////////////////////////////////////// extern "C"{ - class proxy_Actar{ + class proxy_Actar{ public: - proxy_Actar(){ - NPL::DetectorFactory::getInstance()->AddToken("Actar","Actar"); - NPL::DetectorFactory::getInstance()->AddDetector("Actar",TActarPhysics::Construct); - } - }; + proxy_Actar(){ + NPL::DetectorFactory::getInstance()->AddToken("Actar","Actar"); + NPL::DetectorFactory::getInstance()->AddDetector("Actar",TActarPhysics::Construct); + } + }; - proxy_Actar p_Actar; + proxy_Actar p_Actar; } diff --git a/NPLib/Detectors/FissionChamber/CMakeLists.txt b/NPLib/Detectors/FissionChamber/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..b40b59a7c005105be9b3dbfd4bb36a9300608f1c --- /dev/null +++ b/NPLib/Detectors/FissionChamber/CMakeLists.txt @@ -0,0 +1,6 @@ +add_custom_command(OUTPUT TFissionChamberPhysicsDict.cxx COMMAND ../../scripts/build_dict.sh TFissionChamberPhysics.h TFissionChamberPhysicsDict.cxx TFissionChamberPhysics.rootmap libNPFissionChamber.dylib DEPENDS TFissionChamberPhysics.h) +add_custom_command(OUTPUT TFissionChamberDataDict.cxx COMMAND ../../scripts/build_dict.sh TFissionChamberData.h TFissionChamberDataDict.cxx TFissionChamberData.rootmap libNPFissionChamber.dylib DEPENDS TFissionChamberData.h) +add_library(NPFissionChamber SHARED TFissionChamberSpectra.cxx TFissionChamberData.cxx TFissionChamberPhysics.cxx TFissionChamberDataDict.cxx TFissionChamberPhysicsDict.cxx ) +target_link_libraries(NPFissionChamber ${ROOT_LIBRARIES} NPCore) +install(FILES TFissionChamberData.h TFissionChamberPhysics.h TFissionChamberSpectra.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY}) + diff --git a/NPLib/Detectors/FissionChamber/TFissionChamberData.cxx b/NPLib/Detectors/FissionChamber/TFissionChamberData.cxx new file mode 100644 index 0000000000000000000000000000000000000000..13453c44c95eb2aa9ed10e02fd58e7f49ff50303 --- /dev/null +++ b/NPLib/Detectors/FissionChamber/TFissionChamberData.cxx @@ -0,0 +1,79 @@ +/***************************************************************************** + * Copyright (C) 2009-2020 this file is part of the NPTool Project * + * * + * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * + * For the list of contributors see $NPTOOL/Licence/Contributors * + *****************************************************************************/ + +/***************************************************************************** + * Original Author: Pierre Morfouace contact address: pierre.morfouace2@cea.fr * + * * + * Creation Date : September 2020 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class hold FissionChamber Raw data * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ +#include "TFissionChamberData.h" + +#include <iostream> +#include <fstream> +#include <sstream> +#include <string> +using namespace std; + +ClassImp(TFissionChamberData) + + +////////////////////////////////////////////////////////////////////// +TFissionChamberData::TFissionChamberData() { +} + + + +////////////////////////////////////////////////////////////////////// +TFissionChamberData::~TFissionChamberData() { +} + + + +////////////////////////////////////////////////////////////////////// +void TFissionChamberData::Clear() { + // Energy + fFissionChamber_E_DetectorNbr.clear(); + fFissionChamber_Energy.clear(); + // Time + fFissionChamber_T_DetectorNbr.clear(); + fFissionChamber_Time.clear(); +} + + + +////////////////////////////////////////////////////////////////////// +void TFissionChamberData::Dump() const { + // This method is very useful for debuging and worth the dev. + cout << "XXXXXXXXXXXXXXXXXXXXXXXX New Event [TFissionChamberData::Dump()] XXXXXXXXXXXXXXXXX" << endl; + + // Energy + size_t mysize = fFissionChamber_E_DetectorNbr.size(); + cout << "FissionChamber_E_Mult: " << mysize << endl; + + for (size_t i = 0 ; i < mysize ; i++){ + cout << "DetNbr: " << fFissionChamber_E_DetectorNbr[i] + << " Energy: " << fFissionChamber_Energy[i]; + } + + // Time + mysize = fFissionChamber_T_DetectorNbr.size(); + cout << "FissionChamber_T_Mult: " << mysize << endl; + + for (size_t i = 0 ; i < mysize ; i++){ + cout << "DetNbr: " << fFissionChamber_T_DetectorNbr[i] + << " Time: " << fFissionChamber_Time[i]; + } +} diff --git a/NPLib/Detectors/FissionChamber/TFissionChamberData.h b/NPLib/Detectors/FissionChamber/TFissionChamberData.h new file mode 100644 index 0000000000000000000000000000000000000000..637a869eab02bbe196901540734a5558e2e14f6d --- /dev/null +++ b/NPLib/Detectors/FissionChamber/TFissionChamberData.h @@ -0,0 +1,104 @@ +#ifndef __FissionChamberDATA__ +#define __FissionChamberDATA__ +/***************************************************************************** + * Copyright (C) 2009-2020 this file is part of the NPTool Project * + * * + * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * + * For the list of contributors see $NPTOOL/Licence/Contributors * + *****************************************************************************/ + +/***************************************************************************** + * Original Author: Pierre Morfouace contact address: pierre.morfouace2@cea.fr * + * * + * Creation Date : September 2020 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class hold FissionChamber Raw data * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +// STL +#include <vector> +using namespace std; + +// ROOT +#include "TObject.h" + +class TFissionChamberData : public TObject { + ////////////////////////////////////////////////////////////// + // data members are hold into vectors in order + // to allow multiplicity treatment + private: + // Energy + vector<UShort_t> fFissionChamber_E_DetectorNbr; + vector<Double_t> fFissionChamber_Energy; + + // Time + vector<UShort_t> fFissionChamber_T_DetectorNbr; + vector<Double_t> fFissionChamber_Time; + + + ////////////////////////////////////////////////////////////// + // Constructor and destructor + public: + TFissionChamberData(); + ~TFissionChamberData(); + + + ////////////////////////////////////////////////////////////// + // Inherited from TObject and overriden to avoid warnings + public: + void Clear(); + void Clear(const Option_t*) {}; + void Dump() const; + + + ////////////////////////////////////////////////////////////// + // Getters and Setters + // Prefer inline declaration to avoid unnecessary called of + // frequently used methods + // add //! to avoid ROOT creating dictionnary for the methods + public: + ////////////////////// SETTERS //////////////////////// + // Energy + inline void SetEnergy(const UShort_t& DetNbr,const Double_t& Energy){ + fFissionChamber_E_DetectorNbr.push_back(DetNbr); + fFissionChamber_Energy.push_back(Energy); + };//! + + // Time + inline void SetTime(const UShort_t& DetNbr,const Double_t& Time) { + fFissionChamber_T_DetectorNbr.push_back(DetNbr); + fFissionChamber_Time.push_back(Time); + };//! + + + ////////////////////// GETTERS //////////////////////// + // Energy + inline UShort_t GetMultEnergy() const + {return fFissionChamber_E_DetectorNbr.size();} + inline UShort_t GetE_DetectorNbr(const unsigned int &i) const + {return fFissionChamber_E_DetectorNbr[i];}//! + inline Double_t Get_Energy(const unsigned int &i) const + {return fFissionChamber_Energy[i];}//! + + // Time + inline UShort_t GetMultTime() const + {return fFissionChamber_T_DetectorNbr.size();} + inline UShort_t GetT_DetectorNbr(const unsigned int &i) const + {return fFissionChamber_T_DetectorNbr[i];}//! + inline Double_t Get_Time(const unsigned int &i) const + {return fFissionChamber_Time[i];}//! + + + ////////////////////////////////////////////////////////////// + // Required for ROOT dictionnary + ClassDef(TFissionChamberData,1) // FissionChamberData structure +}; + +#endif diff --git a/NPLib/Detectors/FissionChamber/TFissionChamberPhysics.cxx b/NPLib/Detectors/FissionChamber/TFissionChamberPhysics.cxx new file mode 100644 index 0000000000000000000000000000000000000000..88959478444aef21f3ed6494e64e2db8a6b94fd3 --- /dev/null +++ b/NPLib/Detectors/FissionChamber/TFissionChamberPhysics.cxx @@ -0,0 +1,344 @@ +/***************************************************************************** + * Copyright (C) 2009-2020 this file is part of the NPTool Project * + * * + * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * + * For the list of contributors see $NPTOOL/Licence/Contributors * + *****************************************************************************/ + +/***************************************************************************** + * Original Author: Pierre Morfouace contact address: pierre.morfouace2@cea.fr * + * * + * Creation Date : September 2020 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class hold FissionChamber Treated data * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +#include "TFissionChamberPhysics.h" + +// STL +#include <sstream> +#include <iostream> +#include <cmath> +#include <stdlib.h> +#include <limits> +using namespace std; + +// NPL +#include "RootInput.h" +#include "RootOutput.h" +#include "NPDetectorFactory.h" +#include "NPOptionManager.h" + +// ROOT +#include "TChain.h" + +ClassImp(TFissionChamberPhysics) + + +/////////////////////////////////////////////////////////////////////////// +TFissionChamberPhysics::TFissionChamberPhysics() + : m_EventData(new TFissionChamberData), + m_PreTreatedData(new TFissionChamberData), + m_EventPhysics(this), + m_Spectra(0), + m_E_RAW_Threshold(0), // adc channels + m_E_Threshold(0), // MeV + m_NumberOfDetectors(0) { +} + +/////////////////////////////////////////////////////////////////////////// +/// A usefull method to bundle all operation to add a detector +void TFissionChamberPhysics::AddDetector(TVector3 , string ){ + // In That simple case nothing is done + // Typically for more complex detector one would calculate the relevant + // positions (stripped silicon) or angles (gamma array) + m_NumberOfDetectors++; +} + +/////////////////////////////////////////////////////////////////////////// +void TFissionChamberPhysics::AddDetector(double R, double Theta, double Phi, string shape){ + // 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 TFissionChamberPhysics::BuildSimplePhysicalEvent() { + BuildPhysicalEvent(); +} + + + +/////////////////////////////////////////////////////////////////////////// +void TFissionChamberPhysics::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)); + } + } + } +} + +/////////////////////////////////////////////////////////////////////////// +void TFissionChamberPhysics::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("FissionChamber/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("FissionChamber/TIME"+NPL::itoa(m_EventData->GetT_DetectorNbr(i)),m_EventData->Get_Time(i)); + m_PreTreatedData->SetTime(m_EventData->GetT_DetectorNbr(i), Time); + } +} + + + +/////////////////////////////////////////////////////////////////////////// +void TFissionChamberPhysics::ReadAnalysisConfig() { + bool ReadingStatus = false; + + // path to file + string FileName = "./configs/ConfigFissionChamber.dat"; + + // open analysis config file + ifstream AnalysisConfigFile; + AnalysisConfigFile.open(FileName.c_str()); + + if (!AnalysisConfigFile.is_open()) { + cout << " No ConfigFissionChamber.dat found: Default parameter loaded for Analayis " << FileName << endl; + return; + } + cout << " Loading user parameter for Analysis from ConfigFissionChamber.dat " << endl; + + // Save it in a TAsciiFile + TAsciiFile* asciiConfig = RootOutput::getInstance()->GetAsciiFileAnalysisConfig(); + asciiConfig->AppendLine("%%% ConfigFissionChamber.dat %%%"); + asciiConfig->Append(FileName.c_str()); + asciiConfig->AppendLine(""); + // read analysis config file + string LineBuffer,DataBuffer,whatToDo; + while (!AnalysisConfigFile.eof()) { + // Pick-up next line + getline(AnalysisConfigFile, LineBuffer); + + // search for "header" + string name = "ConfigFissionChamber"; + if (LineBuffer.compare(0, name.length(), name) == 0) + ReadingStatus = true; + + // 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=="E_RAW_THRESHOLD") { + AnalysisConfigFile >> DataBuffer; + m_E_RAW_Threshold = atof(DataBuffer.c_str()); + cout << whatToDo << " " << m_E_RAW_Threshold << endl; + } + + else if (whatToDo=="E_THRESHOLD") { + AnalysisConfigFile >> DataBuffer; + m_E_Threshold = atof(DataBuffer.c_str()); + cout << whatToDo << " " << m_E_Threshold << endl; + } + + else { + ReadingStatus = false; + } + } + } +} + + + +/////////////////////////////////////////////////////////////////////////// +void TFissionChamberPhysics::Clear() { + DetectorNumber.clear(); + Energy.clear(); + Time.clear(); +} + + + +/////////////////////////////////////////////////////////////////////////// +void TFissionChamberPhysics::ReadConfiguration(NPL::InputParser parser) { + vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("FissionChamber"); + 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 << "//// FissionChamber " << 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 << "//// FissionChamber " << 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 TFissionChamberPhysics::InitSpectra() { + m_Spectra = new TFissionChamberSpectra(m_NumberOfDetectors); +} + + + +/////////////////////////////////////////////////////////////////////////// +void TFissionChamberPhysics::FillSpectra() { + m_Spectra -> FillRawSpectra(m_EventData); + m_Spectra -> FillPreTreatedSpectra(m_PreTreatedData); + m_Spectra -> FillPhysicsSpectra(m_EventPhysics); +} + + + +/////////////////////////////////////////////////////////////////////////// +void TFissionChamberPhysics::CheckSpectra() { + m_Spectra->CheckSpectra(); +} + + + +/////////////////////////////////////////////////////////////////////////// +void TFissionChamberPhysics::ClearSpectra() { + // To be done +} + + + +/////////////////////////////////////////////////////////////////////////// +map< string , TH1*> TFissionChamberPhysics::GetSpectra() { + if(m_Spectra) + return m_Spectra->GetMapHisto(); + else{ + map< string , TH1*> empty; + return empty; + } +} + +/////////////////////////////////////////////////////////////////////////// +void TFissionChamberPhysics::WriteSpectra() { + m_Spectra->WriteSpectra(); +} + + + +/////////////////////////////////////////////////////////////////////////// +void TFissionChamberPhysics::AddParameterToCalibrationManager() { + CalibrationManager* Cal = CalibrationManager::getInstance(); + for (int i = 0; i < m_NumberOfDetectors; ++i) { + Cal->AddParameter("FissionChamber", "D"+ NPL::itoa(i+1)+"_ENERGY","FissionChamber_D"+ NPL::itoa(i+1)+"_ENERGY"); + Cal->AddParameter("FissionChamber", "D"+ NPL::itoa(i+1)+"_TIME","FissionChamber_D"+ NPL::itoa(i+1)+"_TIME"); + } +} + + + +/////////////////////////////////////////////////////////////////////////// +void TFissionChamberPhysics::InitializeRootInputRaw() { + TChain* inputChain = RootInput::getInstance()->GetChain(); + inputChain->SetBranchStatus("FissionChamber", true ); + inputChain->SetBranchAddress("FissionChamber", &m_EventData ); +} + + + +/////////////////////////////////////////////////////////////////////////// +void TFissionChamberPhysics::InitializeRootInputPhysics() { + TChain* inputChain = RootInput::getInstance()->GetChain(); + inputChain->SetBranchAddress("FissionChamber", &m_EventPhysics); +} + + + +/////////////////////////////////////////////////////////////////////////// +void TFissionChamberPhysics::InitializeRootOutput() { + TTree* outputTree = RootOutput::getInstance()->GetTree(); + outputTree->Branch("FissionChamber", "TFissionChamberPhysics", &m_EventPhysics); +} + + + +//////////////////////////////////////////////////////////////////////////////// +// Construct Method to be pass to the DetectorFactory // +//////////////////////////////////////////////////////////////////////////////// +NPL::VDetector* TFissionChamberPhysics::Construct() { + return (NPL::VDetector*) new TFissionChamberPhysics(); +} + + + +//////////////////////////////////////////////////////////////////////////////// +// Registering the construct method to the factory // +//////////////////////////////////////////////////////////////////////////////// +extern "C"{ +class proxy_FissionChamber{ + public: + proxy_FissionChamber(){ + NPL::DetectorFactory::getInstance()->AddToken("FissionChamber","FissionChamber"); + NPL::DetectorFactory::getInstance()->AddDetector("FissionChamber",TFissionChamberPhysics::Construct); + } +}; + +proxy_FissionChamber p_FissionChamber; +} + diff --git a/NPLib/Detectors/FissionChamber/TFissionChamberPhysics.h b/NPLib/Detectors/FissionChamber/TFissionChamberPhysics.h new file mode 100644 index 0000000000000000000000000000000000000000..7c07cc1a9fbd8032af56a5cbb0ca5c6dba8f24fa --- /dev/null +++ b/NPLib/Detectors/FissionChamber/TFissionChamberPhysics.h @@ -0,0 +1,180 @@ +#ifndef TFissionChamberPHYSICS_H +#define TFissionChamberPHYSICS_H +/***************************************************************************** + * Copyright (C) 2009-2020 this file is part of the NPTool Project * + * * + * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * + * For the list of contributors see $NPTOOL/Licence/Contributors * + *****************************************************************************/ + +/***************************************************************************** + * Original Author: Pierre Morfouace contact address: pierre.morfouace2@cea.fr * + * * + * Creation Date : September 2020 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class hold FissionChamber Treated data * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +// C++ headers +#include <vector> +#include <map> +#include <string> +using namespace std; + +// ROOT headers +#include "TObject.h" +#include "TH1.h" +#include "TVector3.h" +// NPTool headers +#include "TFissionChamberData.h" +#include "TFissionChamberSpectra.h" +#include "NPCalibrationManager.h" +#include "NPVDetector.h" +#include "NPInputParser.h" +// forward declaration +class TFissionChamberSpectra; + + + +class TFissionChamberPhysics : public TObject, public NPL::VDetector { + ////////////////////////////////////////////////////////////// + // constructor and destructor + public: + TFissionChamberPhysics(); + ~TFissionChamberPhysics() {}; + + + ////////////////////////////////////////////////////////////// + // 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 + public: + vector<int> DetectorNumber; + vector<double> Energy; + vector<double> Time; + + /// 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); + + ////////////////////////////////////////////////////////////// + // 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 TFissionChamberSpectra class + // instantiate the TFissionChamberSpectra 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 FissionChamber 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(); + + // give and external TFissionChamberData object to TFissionChamberPhysics. + // needed for online analysis for example + void SetRawDataPointer(TFissionChamberData* rawDataPointer) {m_EventData = rawDataPointer;} + + // objects are not written in the TTree + private: + TFissionChamberData* m_EventData; //! + TFissionChamberData* m_PreTreatedData; //! + TFissionChamberPhysics* m_EventPhysics; //! + + // getters for raw and pre-treated data object + public: + TFissionChamberData* GetRawData() const {return m_EventData;} + TFissionChamberData* GetPreTreatedData() const {return m_PreTreatedData;} + + // parameters used in the analysis + private: + // thresholds + double m_E_RAW_Threshold; //! + double m_E_Threshold; //! + + // number of detectors + private: + int m_NumberOfDetectors; //! + + // spectra class + private: + TFissionChamberSpectra* m_Spectra; // ! + + // spectra getter + public: + map<string, TH1*> GetSpectra(); + + // Static constructor to be passed to the Detector Factory + public: + static NPL::VDetector* Construct(); + + ClassDef(TFissionChamberPhysics,1) // FissionChamberPhysics structure +}; +#endif diff --git a/NPLib/Detectors/FissionChamber/TFissionChamberSpectra.cxx b/NPLib/Detectors/FissionChamber/TFissionChamberSpectra.cxx new file mode 100644 index 0000000000000000000000000000000000000000..9f6f07f65f666ea78ff72ce757a72c02faa36b5b --- /dev/null +++ b/NPLib/Detectors/FissionChamber/TFissionChamberSpectra.cxx @@ -0,0 +1,174 @@ +/***************************************************************************** + * Copyright (C) 2009-2020 this file is part of the NPTool Project * + * * + * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * + * For the list of contributors see $NPTOOL/Licence/Contributors * + *****************************************************************************/ + +/***************************************************************************** + * Original Author: Pierre Morfouace contact address: pierre.morfouace2@cea.fr * + * * + * Creation Date : September 2020 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class hold FissionChamber Spectra * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +// class header +#include "TFissionChamberSpectra.h" + +// STL +#include <iostream> +#include <string> +using namespace std; + +// NPTool header +#include "NPOptionManager.h" + + + +//////////////////////////////////////////////////////////////////////////////// +TFissionChamberSpectra::TFissionChamberSpectra() + : fNumberOfDetectors(0) { + SetName("FissionChamber"); +} + + + +//////////////////////////////////////////////////////////////////////////////// +TFissionChamberSpectra::TFissionChamberSpectra(unsigned int NumberOfDetectors) { + if(NPOptionManager::getInstance()->GetVerboseLevel()>0) + cout << "************************************************" << endl + << "TFissionChamberSpectra : Initalizing control spectra for " + << NumberOfDetectors << " Detectors" << endl + << "************************************************" << endl ; + SetName("FissionChamber"); + fNumberOfDetectors = NumberOfDetectors; + + InitRawSpectra(); + InitPreTreatedSpectra(); + InitPhysicsSpectra(); +} + + + +//////////////////////////////////////////////////////////////////////////////// +TFissionChamberSpectra::~TFissionChamberSpectra() { +} + + + +//////////////////////////////////////////////////////////////////////////////// +void TFissionChamberSpectra::InitRawSpectra() { + static string name; + for (unsigned int i = 0; i < fNumberOfDetectors; i++) { // loop on number of detectors + // Energy + name = "FissionChamber"+NPL::itoa(i+1)+"_ENERGY_RAW"; + AddHisto1D(name, name, 4096, 0, 16384, "FissionChamber/RAW"); + // Time + name = "FissionChamber"+NPL::itoa(i+1)+"_TIME_RAW"; + AddHisto1D(name, name, 4096, 0, 16384, "FissionChamber/RAW"); + } // end loop on number of detectors +} + + + +//////////////////////////////////////////////////////////////////////////////// +void TFissionChamberSpectra::InitPreTreatedSpectra() { + static string name; + for (unsigned int i = 0; i < fNumberOfDetectors; i++) { // loop on number of detectors + // Energy + name = "FissionChamber"+NPL::itoa(i+1)+"_ENERGY_CAL"; + AddHisto1D(name, name, 500, 0, 25, "FissionChamber/CAL"); + // Time + name = "FissionChamber"+NPL::itoa(i+1)+"_TIME_CAL"; + AddHisto1D(name, name, 500, 0, 25, "FissionChamber/CAL"); + + + } // end loop on number of detectors +} + + + +//////////////////////////////////////////////////////////////////////////////// +void TFissionChamberSpectra::InitPhysicsSpectra() { + static string name; + // Kinematic Plot + name = "FissionChamber_ENERGY_TIME"; + AddHisto2D(name, name, 500, 0, 500, 500, 0, 50, "FissionChamber/PHY"); +} + + + +//////////////////////////////////////////////////////////////////////////////// +void TFissionChamberSpectra::FillRawSpectra(TFissionChamberData* RawData) { + static string name; + static string family; + + // Energy + unsigned int sizeE = RawData->GetMultEnergy(); + for (unsigned int i = 0; i < sizeE; i++) { + name = "FissionChamber"+NPL::itoa(RawData->GetE_DetectorNbr(i))+"_ENERGY_RAW"; + family = "FissionChamber/RAW"; + + FillSpectra(family,name,RawData->Get_Energy(i)); + } + + // Time + unsigned int sizeT = RawData->GetMultTime(); + for (unsigned int i = 0; i < sizeT; i++) { + name = "FissionChamber"+NPL::itoa(RawData->GetT_DetectorNbr(i))+"_TIME_RAW"; + family = "FissionChamber/RAW"; + + FillSpectra(family,name,RawData->Get_Time(i)); + } +} + + + +//////////////////////////////////////////////////////////////////////////////// +void TFissionChamberSpectra::FillPreTreatedSpectra(TFissionChamberData* PreTreatedData) { + static string name; + static string family; + + // Energy + unsigned int sizeE = PreTreatedData->GetMultEnergy(); + for (unsigned int i = 0; i < sizeE; i++) { + name = "FissionChamber"+NPL::itoa(PreTreatedData->GetE_DetectorNbr(i))+"_ENERGY_CAL"; + family = "FissionChamber/CAL"; + + FillSpectra(family,name,PreTreatedData->Get_Energy(i)); + } + + // Time + unsigned int sizeT = PreTreatedData->GetMultTime(); + for (unsigned int i = 0; i < sizeT; i++) { + name = "FissionChamber"+NPL::itoa(PreTreatedData->GetT_DetectorNbr(i))+"_TIME_CAL"; + family = "FissionChamber/CAL"; + + FillSpectra(family,name,PreTreatedData->Get_Time(i)); + } +} + + + +//////////////////////////////////////////////////////////////////////////////// +void TFissionChamberSpectra::FillPhysicsSpectra(TFissionChamberPhysics* Physics) { + static string name; + static string family; + family= "FissionChamber/PHY"; + + // Energy vs time + unsigned int sizeE = Physics->Energy.size(); + for(unsigned int i = 0 ; i < sizeE ; i++){ + name = "FissionChamber_ENERGY_TIME"; + FillSpectra(family,name,Physics->Energy[i],Physics->Time[i]); + } +} + diff --git a/NPLib/Detectors/FissionChamber/TFissionChamberSpectra.h b/NPLib/Detectors/FissionChamber/TFissionChamberSpectra.h new file mode 100644 index 0000000000000000000000000000000000000000..e60e5ed22deab295b64d62df3d5e4a58834b52e6 --- /dev/null +++ b/NPLib/Detectors/FissionChamber/TFissionChamberSpectra.h @@ -0,0 +1,62 @@ +#ifndef TFissionChamberSPECTRA_H +#define TFissionChamberSPECTRA_H +/***************************************************************************** + * Copyright (C) 2009-2020 this file is part of the NPTool Project * + * * + * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * + * For the list of contributors see $NPTOOL/Licence/Contributors * + *****************************************************************************/ + +/***************************************************************************** + * Original Author: Pierre Morfouace contact address: pierre.morfouace2@cea.fr * + * * + * Creation Date : September 2020 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class hold FissionChamber Spectra * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +// NPLib headers +#include "NPVSpectra.h" +#include "TFissionChamberData.h" +#include "TFissionChamberPhysics.h" + +// Forward Declaration +class TFissionChamberPhysics; + + +class TFissionChamberSpectra : public VSpectra { + ////////////////////////////////////////////////////////////// + // constructor and destructor + public: + TFissionChamberSpectra(); + TFissionChamberSpectra(unsigned int NumberOfDetectors); + ~TFissionChamberSpectra(); + + ////////////////////////////////////////////////////////////// + // Initialization methods + private: + void InitRawSpectra(); + void InitPreTreatedSpectra(); + void InitPhysicsSpectra(); + + ////////////////////////////////////////////////////////////// + // Filling methods + public: + void FillRawSpectra(TFissionChamberData*); + void FillPreTreatedSpectra(TFissionChamberData*); + void FillPhysicsSpectra(TFissionChamberPhysics*); + + ////////////////////////////////////////////////////////////// + // Detector parameters + private: + unsigned int fNumberOfDetectors; +}; + +#endif diff --git a/NPLib/TrackReconstruction/NPCluster.cxx b/NPLib/TrackReconstruction/NPCluster.cxx index 77421d2d4210a99b1bfbdda2ceb4fb7aae532fef..6face975737d39cb8fb9f9bd486803a31cc8f760 100644 --- a/NPLib/TrackReconstruction/NPCluster.cxx +++ b/NPLib/TrackReconstruction/NPCluster.cxx @@ -23,58 +23,58 @@ using namespace std; ClassImp(Cluster) -///////////////////////////////////////////////// + ///////////////////////////////////////////////// Cluster::Cluster() { - fNumberOfPadsX = 128; - fNumberOfPadsY = 128; - fNumberOfSample = 512; - fNumberOfTracksMax = 20; - fClusterDistance = 5; - fClusterThreshold = 15; + fNumberOfPadsX = 128; + fNumberOfPadsY = 128; + fNumberOfSample = 512; + fNumberOfTracksMax = 20; + fClusterDistance = 5; + fClusterThreshold = 15; - fVisu= 0; + fVisu= 0; - LocBigVoxels = new int[BigVoxelSize]; + LocBigVoxels = new int[BigVoxelSize]; } ///////////////////////////////////////////////// Cluster::Cluster(int NumberOfPadsX, int NumberOfPadsY, bool Visu) { - fNumberOfPadsX = NumberOfPadsX; - fNumberOfPadsY = NumberOfPadsY; - fNumberOfSample = 512; - 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(); - if(!m_ServerSocket->IsValid()){ - cout << endl; - cout << "/// ServerSocket invalid! ///" << endl; - exit(1); - } - cout << endl; - cout << "/// ServerSocket valid ///" << endl; - m_ServerSocket->SetCompressionSettings(1); + fNumberOfPadsX = NumberOfPadsX; + fNumberOfPadsY = NumberOfPadsY; + fNumberOfSample = 512; + 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(); + if(!m_ServerSocket->IsValid()){ + cout << endl; + cout << "/// ServerSocket invalid! ///" << endl; + exit(1); + } + 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; + // 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; + // Create a list to contain all client connections + m_Histo = new TList; - } + } } @@ -86,73 +86,74 @@ 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; - cout << "/// Number of tracks max= " << fNumberOfTracksMax << endl; - } - else{ - string LineBuffer, whatToDo, DataBuffer; - while(!ClusterParamFile.eof()){ - getline(ClusterParamFile,LineBuffer); - string name = "ConfigCluster"; - if(LineBuffer.compare(0, name.length(), name) == 0){ - cout << endl; - cout << "**** ConfigCluster found" << endl; - ReadingStatus = true; - } - while(ReadingStatus){ - whatToDo=""; - ClusterParamFile >> whatToDo; - // Search for comment symbol (%) - if (whatToDo.compare(0, 1, "%") == 0) { - ClusterParamFile.ignore(numeric_limits<streamsize>::max(), '\n' ); - } + 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; + cout << "/// Number of tracks max= " << fNumberOfTracksMax << endl; + } + else{ + string LineBuffer, whatToDo, DataBuffer; + while(!ClusterParamFile.eof()){ + getline(ClusterParamFile,LineBuffer); + string name = "ConfigCluster"; + if(LineBuffer.compare(0, name.length(), name) == 0){ + cout << endl; + cout << "**** ConfigCluster found" << endl; + ReadingStatus = true; + } + while(ReadingStatus){ + whatToDo=""; + ClusterParamFile >> whatToDo; + // Search for comment symbol (%) + 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()); - cout << "/// Number of tracks max= " << " " << fNumberOfTracksMax << " ///" << endl; - } - else if (whatToDo.compare(0,16,"ClusterDistance=") == 0) { - ClusterParamFile >> DataBuffer; - fClusterDistance = atoi(DataBuffer.c_str()); - cout << "/// Cluster Distance max= " << " " << fClusterDistance << " ///" << endl; - } - else if (whatToDo.compare(0,17,"ClusterThreshold=") == 0) { - ClusterParamFile >> DataBuffer; - fClusterThreshold = atoi(DataBuffer.c_str()); - cout << "/// Cluster Threshold= " << " " << fClusterThreshold << " ///" << endl; - } - else{ - ReadingStatus = false; - } - } + else if (whatToDo.compare(0,18,"NumberOfTracksMax=") == 0) { + ClusterParamFile >> DataBuffer; + fNumberOfTracksMax = atoi(DataBuffer.c_str()); + cout << "/// Number of tracks max= " << " " << fNumberOfTracksMax << " ///" << endl; + } + else if (whatToDo.compare(0,16,"ClusterDistance=") == 0) { + ClusterParamFile >> DataBuffer; + fClusterDistance = atoi(DataBuffer.c_str()); + cout << "/// Cluster Distance max= " << " " << fClusterDistance << " ///" << endl; } + else if (whatToDo.compare(0,17,"ClusterThreshold=") == 0) { + ClusterParamFile >> DataBuffer; + fClusterThreshold = atoi(DataBuffer.c_str()); + cout << "/// Cluster Threshold= " << " " << fClusterThreshold << " ///" << endl; + } + else{ + ReadingStatus = false; + } + } } + } } ////////////////////////////////////////////////// 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]); - } + 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(); + RegroupZone(); - FitZone(); + FitZone(); - RegroupTrack(); + RegroupTrack(); + + FitTrack(); - FitTrack(); } ////////////////////////////////////////////////// @@ -163,100 +164,102 @@ void Cluster::FindTracks() ////////////////////////////////////////////////// void Cluster::FillBigVoxels(int x, int y, int t, int q) { + if(t<512){ 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++; + 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); + NPL::Track BigTrack; + BigTrack.SetPointX(x/4); + BigTrack.SetPointY(y/4); + BigTrack.SetPointZ(t/8); - BigTrack.SetPointX(xB); - BigTrack.SetPointY(yB); - BigTrack.SetPointZ(tB); + 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(); + 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]; + } + 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); + 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]; + } + 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); + 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]); + } + 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 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(); + 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]; + 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); + 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); + 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]; + LocBigVoxels[element_max]=LocBigVoxels[element_min]; - 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(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(); + 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++){ + 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]--; - } - TrackZone.at(i-1) = TrackZone.at(i); - } - it--; - TrackZone.at(TrackZone.size()-1).Clear(); - TrackZone.pop_back(); - } + 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]--; + } + TrackZone.at(i-1) = TrackZone.at(i); } + it--; + TrackZone.at(TrackZone.size()-1).Clear(); + TrackZone.pop_back(); + } } + } } + } } ////////////////////////////////////////////////// @@ -264,51 +267,51 @@ void Cluster::FillBigVoxels(int x, int y, int t, int q) ////////////////////////////////////////////////// void Cluster::RegroupZone() { - int x1,y1,z1,x2,y2,z2; - if(TrackZone.size()>0){ - for(int i=0;i<TrackZone.size()-1;i++){ - for(int j=i+1;j<TrackZone.size();j++){ - it=0; - for(int k1=0;k1<TrackZone.at(i).GetXPoints().size();k1++){ - for(int k2=0;k2<TrackZone.at(j).GetXPoints().size();k2++){ - // Extract the coordinates of each big voxel of 2 TrackZone to check their closeness. - x1 = TrackZone.at(i).GetPointX(k1); - y1 = TrackZone.at(i).GetPointY(k1); - z1 = TrackZone.at(i).GetPointZ(k1); - - 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++; - } - if(it>0){ - //Threshold to decide when two TrackZone are actually only one or not. - for(int k2=0;k2<TrackZone.at(j).GetXPoints().size();k2++){ - //Add the neighbour TrackZone j to the considered TrackZone i. - TrackZone.at(i).SetPointX(TrackZone.at(j).GetPointX(k2)); - TrackZone.at(i).SetPointY(TrackZone.at(j).GetPointY(k2)); - TrackZone.at(i).SetPointZ(TrackZone.at(j).GetPointZ(k2)); - } - for(int k2=0;k2<TrackZone.at(j).GetXPoints().size();k2++){ - // Get rid of the remnant of this TrackZone j integrated to the afore mentioned TrackZone i. - TrackZone.at(j).GetXPoints().pop_back(); - TrackZone.at(j).GetYPoints().pop_back(); - TrackZone.at(j).GetZPoints().pop_back(); - } - RemoveTrack(TrackZone,j); - } - } + int x1,y1,z1,x2,y2,z2; + if(TrackZone.size()>0){ + for(int i=0;i<TrackZone.size()-1;i++){ + for(int j=i+1;j<TrackZone.size();j++){ + it=0; + for(int k1=0;k1<TrackZone.at(i).GetXPoints().size();k1++){ + for(int k2=0;k2<TrackZone.at(j).GetXPoints().size();k2++){ + // Extract the coordinates of each big voxel of 2 TrackZone to check their closeness. + x1 = TrackZone.at(i).GetPointX(k1); + y1 = TrackZone.at(i).GetPointY(k1); + z1 = TrackZone.at(i).GetPointZ(k1); + + 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++; + } + if(it>0){ + //Threshold to decide when two TrackZone are actually only one or not. + for(int k2=0;k2<TrackZone.at(j).GetXPoints().size();k2++){ + //Add the neighbour TrackZone j to the considered TrackZone i. + TrackZone.at(i).SetPointX(TrackZone.at(j).GetPointX(k2)); + TrackZone.at(i).SetPointY(TrackZone.at(j).GetPointY(k2)); + TrackZone.at(i).SetPointZ(TrackZone.at(j).GetPointZ(k2)); + } + for(int k2=0;k2<TrackZone.at(j).GetXPoints().size();k2++){ + // Get rid of the remnant of this TrackZone j integrated to the afore mentioned TrackZone i. + TrackZone.at(j).GetXPoints().pop_back(); + TrackZone.at(j).GetYPoints().pop_back(); + TrackZone.at(j).GetZPoints().pop_back(); } + RemoveTrack(TrackZone,j); + } } + } } - for(int i=0;i<TrackZone.size();i++){ - if(TrackZone.at(i).GetXPoints().size()<2){ - RemoveTrack(TrackZone,i); - i--; - } + } + for(int i=0;i<TrackZone.size();i++){ + if(TrackZone.at(i).GetXPoints().size()<2){ + RemoveTrack(TrackZone,i); + i--; } + } } @@ -320,38 +323,38 @@ void Cluster::RegroupZone() ////////////////////////////////////////////////// void Cluster::FitZone() { - double p1xy, p1xz, p1yz, p0xy, p0xz, p0yz; - for(unsigned int i = 0; i<TrackZone.size(); i++){ - 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); - myTrack.SetPointQ(q); - } - } - if(myTrack.GetXPoints().size()>0){ - Fit3D(myTrack); - - TrackZone[i].SetXm(myTrack.GetXm()); - TrackZone[i].SetYm(myTrack.GetYm()); - TrackZone[i].SetZm(myTrack.GetZm()); - TrackZone[i].SetXh(myTrack.GetXh()); - TrackZone[i].SetYh(myTrack.GetYh()); - TrackZone[i].SetZh(myTrack.GetZh()); - TrackZone[i].SetSlopesAndOffsets(); - } - myTrack.Clear(); + double p1xy, p1xz, p1yz, p0xy, p0xz, p0yz; + for(unsigned int i = 0; i<TrackZone.size(); i++){ + 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); + myTrack.SetPointQ(q); + } } + if(myTrack.GetXPoints().size()>0){ + Fit3D(myTrack); + + TrackZone[i].SetXm(myTrack.GetXm()); + TrackZone[i].SetYm(myTrack.GetYm()); + TrackZone[i].SetZm(myTrack.GetZm()); + TrackZone[i].SetXh(myTrack.GetXh()); + TrackZone[i].SetYh(myTrack.GetYh()); + TrackZone[i].SetZh(myTrack.GetZh()); + TrackZone[i].SetSlopesAndOffsets(); + } + myTrack.Clear(); + } } ////////////////////////////////////////////////// @@ -359,277 +362,277 @@ 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; - int zT, zTZ; - int element1 = TrackZone.at(id).GetPointX(0) + fNumberOfPadsX/4*TrackZone.at(id).GetPointY(0) + fNumberOfPadsX/4*fNumberOfPadsY/4*TrackZone.at(id).GetPointZ(0); - ExtractCoordinate(IDBigVoxels[element1].at(0),xT,yT,zT); - for(unsigned int i=0; i<TrackZoneSize; i++){ - for(unsigned int j=0; j<TrackZone[i].GetXPoints().size(); j++){ - 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(); - for(unsigned int k=0; k<IDBigVoxels[element2].size(); k++){ - if(IDBigVoxels[element2].at(k)>0){ - 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); - newTrack.SetPointZ(z_track); - newTrack.SetPointQ(q_track); - } - } - else if(i==id && IDBigVoxels[element2].at(k)==-1) ite++; - } - } + 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; + int zT, zTZ; + int element1 = TrackZone.at(id).GetPointX(0) + fNumberOfPadsX/4*TrackZone.at(id).GetPointY(0) + fNumberOfPadsX/4*fNumberOfPadsY/4*TrackZone.at(id).GetPointZ(0); + ExtractCoordinate(IDBigVoxels[element1].at(0),xT,yT,zT); + for(unsigned int i=0; i<TrackZoneSize; i++){ + for(unsigned int j=0; j<TrackZone[i].GetXPoints().size(); j++){ + 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(); + for(unsigned int k=0; k<IDBigVoxels[element2].size(); k++){ + if(IDBigVoxels[element2].at(k)>0){ + 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); + newTrack.SetPointZ(z_track); + newTrack.SetPointQ(q_track); } + } + else if(i==id && IDBigVoxels[element2].at(k)==-1) ite++; } - if((float)ite>=TrackZoneNumberOfPoints*0.9){ - vTrack.push_back(newTrack); - } - - newTrack.Clear(); + } } + } + if((float)ite>=TrackZoneNumberOfPoints*0.9){ + vTrack.push_back(newTrack); + } + + newTrack.Clear(); } + } } ////////////////////////////////////////////////// void Cluster::FitTrack() { - vector<int> to_erase; - for(unsigned int i=0; i<vTrack.size(); i++){ - if(vTrack[i].GetXPoints().size()>fClusterThreshold){ - Fit3D(vTrack[i]); - } - else to_erase.push_back(i); + vector<int> to_erase; + for(unsigned int i=0; i<vTrack.size(); i++){ + if(vTrack[i].GetXPoints().size()>fClusterThreshold){ + Fit3D(vTrack[i]); } + else to_erase.push_back(i); + } - for(int i=to_erase.size()-1; i>=0; i--) vTrack.erase(vTrack.begin()+to_erase[i]); + for(int i=to_erase.size()-1; i>=0; i--) vTrack.erase(vTrack.begin()+to_erase[i]); } ////////////////////////////////////////////////// double Cluster::Distance3D(NPL::Track aTrack, int ID) { - int x, y, z; - double xf, yf, zf; - double distance; - ExtractCoordinate(ID,x,y,z); + int x, y, z; + double xf, yf, zf; + double distance; + ExtractCoordinate(ID,x,y,z); - TVector3 V1; - TVector3 V2; + 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); + 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); + TVector3 V2CV1 = V2.Cross(V1); - distance = V2CV1.Mag()/V1.Mag(); + distance = V2CV1.Mag()/V1.Mag(); - return distance; + return distance; } ////////////////////////////////////////////////// void Cluster::ReorderTracks(vector<NPL::Track> &someTracks) { - NPL::Track MoveTrack; - if(someTracks.size()>1){ - for(unsigned int i=0; i<someTracks.size(); i++){ - for(int j=i+1; j<someTracks.size(); j++){ - if(someTracks.at(j).GetXPoints().size()>someTracks.at(i).GetXPoints().size()){ - MoveTrack = someTracks[j]; - someTracks[j] = someTracks[i]; - someTracks[i] = MoveTrack; - } - } + NPL::Track MoveTrack; + if(someTracks.size()>1){ + for(unsigned int i=0; i<someTracks.size(); i++){ + for(int j=i+1; j<someTracks.size(); j++){ + if(someTracks.at(j).GetXPoints().size()>someTracks.at(i).GetXPoints().size()){ + MoveTrack = someTracks[j]; + someTracks[j] = someTracks[i]; + someTracks[i] = MoveTrack; } + } } + } } ////////////////////////////////////////////////// void Cluster::ExtractCoordinate(int ID, int& x, int &y, int& z) { - x = ID%fNumberOfPadsX; - y = (ID%(fNumberOfPadsX*fNumberOfPadsY)-x)/fNumberOfPadsX; - z = ID/fNumberOfPadsX/fNumberOfPadsY; + x = ID%fNumberOfPadsX; + y = (ID%(fNumberOfPadsX*fNumberOfPadsY)-x)/fNumberOfPadsX; + z = ID/fNumberOfPadsX/fNumberOfPadsY; } ////////////////////////////////////////////////// void Cluster::RemoveTrack(vector<NPL::Track> &Ephem, int id) { - if(Ephem.size()>id) - for(int i=id;i<Ephem.size()-1;i++){ - Ephem.at(i) = Ephem.at(i+1); - } - Ephem.pop_back(); + if(Ephem.size()>id) + for(int i=id;i<Ephem.size()-1;i++){ + Ephem.at(i) = Ephem.at(i+1); + } + Ephem.pop_back(); } ////////////////////////////////////////////////// void Cluster::Init(vector<int> v1, vector<int> v2, vector<double> v3, vector<double> v4) { - Reset(); + Reset(); - vX = v1; - vY = v2; - vZ = v3; - vQ = v4; + vX = v1; + vY = v2; + vZ = v3; + vQ = v4; - fOriginalCloudSize = vX.size(); - double TotalCharge=0; - for(unsigned int i=0; i< vQ.size(); i++){ - TotalCharge += vQ[i]; - } - fTotalCharge = TotalCharge; + fOriginalCloudSize = vX.size(); + double TotalCharge=0; + for(unsigned int i=0; i< vQ.size(); i++){ + TotalCharge += vQ[i]; + } + fTotalCharge = TotalCharge; } ////////////////////////////////////////////////// void Cluster::Reset() { - fTotalCharge = 0; - vX.clear(); - vY.clear(); - vZ.clear(); - vQ.clear(); - TrackZone.clear(); - - for(int i=0; i<BigVoxelSize; i++){ - BigVoxels[i].clear(); - IDBigVoxels[i].clear(); - LocBigVoxels[i] = -1; - } - it = -1; + fTotalCharge = 0; + vX.clear(); + vY.clear(); + vZ.clear(); + vQ.clear(); + TrackZone.clear(); + + for(int i=0; i<BigVoxelSize; i++){ + BigVoxels[i].clear(); + IDBigVoxels[i].clear(); + LocBigVoxels[i] = -1; + } + it = -1; } ////////////////////////////////////////////////// void Cluster::Clear() { - vTrack.clear(); + vTrack.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; - double Xh,Yh,Zh; - double a,b; - double Sxx,Sxy,Syy,Sxz,Szz,Syz; - double theta; - double K11,K22,K12,K10,K01,K00; - 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++; - Xm+=X[i]; - Ym+=Y[i]; - Zm+=Z[i]; - Sxx+=X[i]*X[i]; - Syy+=Y[i]*Y[i]; - Szz+=Z[i]*Z[i]; - Sxy+=X[i]*Y[i]; - Sxz+=X[i]*Z[i]; - Syz+=Y[i]*Z[i]; - } - - Xm/=Q; - Ym/=Q; - Zm/=Q; - Sxx/=Q; - Syy/=Q; - Szz/=Q; - Sxy/=Q; - Sxz/=Q; - Syz/=Q; - Sxx-=(Xm*Xm); - Syy-=(Ym*Ym); - Szz-=(Zm*Zm); - 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) - { - rho=sqrt(-pow(p,3)/27.); - 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); - aTrack.SetXh(Xh); - aTrack.SetYh(Yh); - aTrack.SetZh(Zh); - aTrack.SetSlopesAndOffsets(); - - - return(fabs(dm2/Q)); + // 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; + double Xh,Yh,Zh; + double a,b; + double Sxx,Sxy,Syy,Sxz,Szz,Syz; + double theta; + double K11,K22,K12,K10,K01,K00; + 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++; + Xm+=X[i]; + Ym+=Y[i]; + Zm+=Z[i]; + Sxx+=X[i]*X[i]; + Syy+=Y[i]*Y[i]; + Szz+=Z[i]*Z[i]; + Sxy+=X[i]*Y[i]; + Sxz+=X[i]*Z[i]; + Syz+=Y[i]*Z[i]; + } + + Xm/=Q; + Ym/=Q; + Zm/=Q; + Sxx/=Q; + Syy/=Q; + Szz/=Q; + Sxy/=Q; + Sxz/=Q; + Syz/=Q; + Sxx-=(Xm*Xm); + Syy-=(Ym*Ym); + Szz-=(Zm*Zm); + 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) + { + rho=sqrt(-pow(p,3)/27.); + 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); + aTrack.SetXh(Xh); + aTrack.SetYh(Yh); + aTrack.SetZh(Zh); + aTrack.SetSlopesAndOffsets(); + + + return(fabs(dm2/Q)); } diff --git a/NPSimulation/Core/MaterialManager.cc b/NPSimulation/Core/MaterialManager.cc index e111e5f16c7d0f03596d55eaa27f13830fa5109d..776b66f51d583f861b76e425a35087a48e36c763 100644 --- a/NPSimulation/Core/MaterialManager.cc +++ b/NPSimulation/Core/MaterialManager.cc @@ -973,6 +973,17 @@ G4Material* MaterialManager::GetMaterialFromLibrary(string Name, G4Material* material = new G4Material("NPS_" + Name, density, 2); material->AddElement(GetElementFromLibrary("C"), 5); material->AddElement(GetElementFromLibrary("H"), 6); + + m_Material[Name] = material; + return material; + } + + else if (Name == "NE213_optical") { + if (!density) + density = 0.874 * g / cm3; + G4Material* material = new G4Material("NPS_" + Name, density, 2); + material->AddElement(GetElementFromLibrary("C"), 5); + material->AddElement(GetElementFromLibrary("H"), 6); //--------------------- Optical Properties ---------------------// const G4int NUMENTRIES = 15; diff --git a/NPSimulation/Detectors/Actar/Actar.cc b/NPSimulation/Detectors/Actar/Actar.cc index b45b90f881f75cc4842112e0a6a94a2c7a16b55e..3691d550d311eea3036da6ab4e4e122a8ce8e284 100644 --- a/NPSimulation/Detectors/Actar/Actar.cc +++ b/NPSimulation/Detectors/Actar/Actar.cc @@ -152,7 +152,11 @@ Actar::Actar(){ m_ReactionRegion=NULL; // Lookup table // // bool ReadingLookupTable = false; - string LT_FileName = "./Detectors/Actar/LT.dat"; + // Opening the LookUp Table LT.dat + G4String GlobalPath = getenv("NPTOOL"); + G4String LT_FileName = GlobalPath + "/NPSimulation/Detectors/Actar/LT.dat"; + //string LT_FileName = "./Detectors/Actar/LT.dat"; + //string LT_FileName = "./configs/LT.dat"; ifstream LTConfigFile; LTConfigFile.open(LT_FileName.c_str()); if(!LTConfigFile.is_open()){ @@ -274,7 +278,7 @@ G4LogicalVolume* Actar::BuildDetector(){ MPT->AddConstProperty("DE_ABSLENGTH",1*pc); MPT->AddConstProperty("DE_DRIFTSPEED",0.8*cm/microsecond); MPT->AddConstProperty("DE_TRANSVERSALSPREAD",2e-5*mm2/ns); - MPT->AddConstProperty("DE_LONGITUDINALSPREAD",7e-5*mm2/ns); + MPT->AddConstProperty("DE_LONGITUDINALSPREAD",4e-6*mm2/ns); DriftGasMaterial->SetMaterialPropertiesTable(MPT); @@ -655,6 +659,7 @@ void Actar::ReadSensitive(const G4Event*){ DataReduced.clear(); double Count = RandGauss::shoot(PadScorer->GetCharge(i),Actar_NS::ResoCharge*PadScorer->GetCharge(i)); double Time = PadScorer->GetTime(i);//RandGauss::shoot(Info[1],Actar_NS::ResoTime); + //cout << "Time= " << Time << endl; //int iTime = ((int) Time*20/512)+1; //int PadNbr = Info[7]; int Pad_X = PadScorer->GetPadX(i);//m_PadToXRow[PadNbr]; @@ -664,7 +669,7 @@ void Actar::ReadSensitive(const G4Event*){ int as = m_PadToAsad[Pad_X][Pad_Y]; int ag = m_PadToAGET[Pad_X][Pad_Y]; int ch = m_PadToChannel[Pad_X][Pad_Y]; - + if(Count>Actar_NS::ChargeThreshold){ DataReduced.globalchannelid = ch+(ag<<7)+(as<<9)+(co<<11); DataReduced.peakheight.push_back(Count); @@ -672,7 +677,7 @@ void Actar::ReadSensitive(const G4Event*){ } m_EventReduced->CoboAsad.push_back(DataReduced); } - + /*vector<double> Q, T; for(int i=0; i<h->GetNbinsX(); i++){ double count = h->GetBinContent(i); diff --git a/NPSimulation/Detectors/ChiNu/ChiNu.cc b/NPSimulation/Detectors/ChiNu/ChiNu.cc index 2f2eaa5b7a1486b8d4945836fd615a5f9a940a79..f8a2b108afc5759a27428894be0eacdb898d7afe 100644 --- a/NPSimulation/Detectors/ChiNu/ChiNu.cc +++ b/NPSimulation/Detectors/ChiNu/ChiNu.cc @@ -88,12 +88,6 @@ namespace ChiNu_NS{ // Lead shield const double Lead_Radius = 9*cm; const double Lead_Thickness = 2*mm; - - // Fission Chamber - const string FCWall_Material = "CH2"; - const double Cu_Thickness = 17*micrometer; - // const double Al_Thickness = 12*micrometer; - const double Kapton_Thickness = 50*micrometer; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... @@ -110,8 +104,6 @@ ChiNu::ChiNu(){ m_BuildLeadShield = 0; - m_FissionChamberWall = 0; - m_FissionChamberVolume = 0; // RGB Color + Transparency m_VisCylinder = new G4VisAttributes(G4Colour(0.0, 0.5, 1, 1)); @@ -143,178 +135,7 @@ void ChiNu::AddDetector(double R, double Theta, double Phi){ m_Theta.push_back(Theta); m_Phi.push_back(Phi); } - -//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... -G4AssemblyVolume* ChiNu::BuildFissionChamber(){ - if(!m_FissionChamberVolume){ - m_FissionChamberVolume = new G4AssemblyVolume(); - - G4RotationMatrix *Rv=new G4RotationMatrix(0,0,0); - G4ThreeVector Tv; - Tv.setX(0); Tv.setY(0); Tv.setZ(0); - - // Bottom PCB plate // - double PCB_width = 18.*cm; - double PCB_length = 33.*cm; - double PCB_Rogers_height = 1.6*mm; - double PCB_Cu_height = 6*35.*um; - double PCB_PosY = -8.5*cm; - // Cu layers - G4Box* PCB_Cu_solid = new G4Box("PCB_Cu_solid",0.5*PCB_width,0.5*PCB_Cu_height,0.5*PCB_length); - G4Material* Cu_material = MaterialManager::getInstance()->GetMaterialFromLibrary("Cu"); - G4LogicalVolume* PCB_Cu_vol = new G4LogicalVolume(PCB_Cu_solid, Cu_material,"PCB_Cu_logic",0,0,0); - PCB_Cu_vol->SetVisAttributes(m_VisCu); - Tv.setY(PCB_PosY); - m_FissionChamberVolume->AddPlacedVolume(PCB_Cu_vol, Tv, Rv); - - // Rogers 4003C layers - G4Box* PCB_Rogers_solid = new G4Box("PCB_Rogers_solid",0.5*PCB_width,0.5*PCB_Rogers_height,0.5*PCB_length); - G4Material* Rogers_material = MaterialManager::getInstance()->GetMaterialFromLibrary("Rogers4003C"); - G4LogicalVolume* PCB_Rogers_vol = new G4LogicalVolume(PCB_Rogers_solid, Rogers_material,"PCB_Rogers_logic",0,0,0); - PCB_Rogers_vol->SetVisAttributes(m_VisRogers4003C); - Tv.setY(PCB_PosY + 0.5*PCB_Cu_height + 0.5*PCB_Rogers_height); - m_FissionChamberVolume->AddPlacedVolume(PCB_Rogers_vol, Tv, Rv); - - //Al frame // - double frame1_width = 18.*cm; - double frame1_height = 5.*mm; - double frame1_length = 33.*cm; - double frame2_width = 15.2*cm; - double frame2_height = 5.1*mm; - double frame2_length = 30.2*cm; - - G4Box* frame1 = new G4Box("frame1", 0.5*frame1_width, 0.5*frame1_height, 0.5*frame1_length); - G4Box* frame2 = new G4Box("frame2", 0.5*frame2_width, 0.5*frame2_height, 0.5*frame2_length); - G4VSolid* Al_frame = (G4VSolid*) new G4SubtractionSolid("Al_frame",frame1,frame2,0,G4ThreeVector(0,0,0)); - G4Material* Al_material = MaterialManager::getInstance()->GetMaterialFromLibrary("Al"); - G4LogicalVolume* Al_frame_vol = new G4LogicalVolume(Al_frame,Al_material,"Al_frame_logic",0,0,0); - Al_frame_vol->SetVisAttributes(m_VisAl); - Tv.setY(PCB_PosY+ 0.5*PCB_Cu_height + PCB_Rogers_height + 0.5*frame1_height); - m_FissionChamberVolume->AddPlacedVolume(Al_frame_vol, Tv, Rv); - Tv.setY(PCB_PosY- 0.5*PCB_Cu_height - 0.5*frame1_height); - m_FissionChamberVolume->AddPlacedVolume(Al_frame_vol, Tv, Rv); - - double box1_width = 15.*cm; - double box1_height = 16.*cm; - double box1_length = 30.*cm; - double box2_width = 14.8*cm; - double box2_height = 15.8*cm; - double box2_length = 29.8*cm; - double box3_width = 12.5*cm; - double box3_height = 11.7*cm; - double box3_length = 30.1*cm; - double box4_width = 15.1*cm; - double box4_height = 11.8*cm; - double box4_length = 27.4*cm; - double box5_width = 12.4*cm; - double box5_height = 16.1*cm; - double box5_length = 27.4*cm; - - G4Box* box1 = new G4Box("box1", 0.5*box1_width, 0.5*box1_height, 0.5*box1_length); - G4Box* box2 = new G4Box("box2", 0.5*box2_width, 0.5*box2_height, 0.5*box2_length); - G4Box* box3 = new G4Box("box3", 0.5*box3_width, 0.5*box3_height, 0.5*box3_length); - G4Box* box4 = new G4Box("box4", 0.5*box4_width, 0.5*box4_height, 0.5*box4_length); - G4Box* box5 = new G4Box("box5", 0.5*box5_width, 0.5*box5_height, 0.5*box5_length); - - G4VSolid* box_int1 = (G4VSolid*) new G4SubtractionSolid("box_int1",box1,box2,0,G4ThreeVector(0,0,0)); - G4VSolid* box_int2 = (G4VSolid*) new G4SubtractionSolid("box_int2",box_int1,box3,0,G4ThreeVector(0,0,0)); - G4VSolid* box_int3 = (G4VSolid*) new G4SubtractionSolid("box_int3",box_int2,box4,0,G4ThreeVector(0,0,0)); - G4VSolid* box_int4 = (G4VSolid*) new G4SubtractionSolid("box_int4",box_int3,box5,0,G4ThreeVector(0,0,0)); - - G4LogicalVolume* full_box_vol = new G4LogicalVolume(box_int4, Al_material, "full_box_logic", 0,0,0); - full_box_vol->SetVisAttributes(m_VisAl); - Tv.setY(0); - m_FissionChamberVolume->AddPlacedVolume(full_box_vol, Tv, Rv); - - // Ti foils // - double foil1_width = 13*cm; - double foil1_length = 29*cm; - double foil1_thickness = 100*um; - double foil2_width = 13*cm; - double foil2_height = 14*cm; - double foil2_thickness = 50*um; - - G4Material* Ti_material = MaterialManager::getInstance()->GetMaterialFromLibrary("Ti"); - G4Box* foil1_solid = new G4Box("foil1", 0.5*foil1_width, 0.5*foil1_thickness, 0.5*foil1_length); - G4LogicalVolume* foil1_vol = new G4LogicalVolume(foil1_solid, Ti_material, "foil1_logic", 0, 0, 0); - foil1_vol->SetVisAttributes(m_VisTi); - Tv.setY(0.5*box2_height); - m_FissionChamberVolume->AddPlacedVolume(foil1_vol, Tv, Rv); - Tv.setY(0); - Tv.setX(-0.5*box2_width); - Rv->rotateZ(90*deg); - m_FissionChamberVolume->AddPlacedVolume(foil1_vol, Tv, Rv); - Tv.setX(0.5*box2_width); - m_FissionChamberVolume->AddPlacedVolume(foil1_vol, Tv, Rv); - - G4Box* foil2_solid = new G4Box("foil2", 0.5*foil2_width, 0.5*foil2_height, 0.5*foil2_thickness); - G4LogicalVolume* foil2_vol = new G4LogicalVolume(foil2_solid, Ti_material, "foil2_logic", 0, 0, 0); - foil2_vol->SetVisAttributes(m_VisTi); - Tv.setX(0);Tv.setY(0);Tv.setZ(-0.5*box2_length); - m_FissionChamberVolume->AddPlacedVolume(foil2_vol, Tv, Rv); - Tv.setZ(0.5*box2_length); - m_FissionChamberVolume->AddPlacedVolume(foil2_vol, Tv, Rv); - - // Cathode and Anode // - BuildCathode(-27.5); - double origine_anode = -25*mm; - double origine_cathode = -22.5*mm; - for(int i=0; i<11; i++){ - BuildAnode(origine_anode+i*5*mm); - } - for(int i=0; i<10; i++){ - BuildCathode(origine_cathode+i*5*mm); - } - BuildCathode(27.5); - - - } - return m_FissionChamberVolume; -} -//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... -void ChiNu::BuildCathode(double Zpos){ - // Al plate: 12 um - G4Tubs* Al_plate_solid = new G4Tubs("Al_plate",0,40*mm,12*micrometer,0,360*deg); - G4Material* Al_material = MaterialManager::getInstance()->GetMaterialFromLibrary("Al"); - G4LogicalVolume* Al_vol = new G4LogicalVolume(Al_plate_solid, Al_material,"logic_Al",0,0,0); - Al_vol->SetVisAttributes(m_VisAl); - - G4RotationMatrix *Rv=new G4RotationMatrix(0,0,0); - G4ThreeVector Tv; - Tv.setX(0); Tv.setY(0); Tv.setZ(0); - - Tv.setZ(Zpos); - m_FissionChamberVolume->AddPlacedVolume(Al_vol, Tv, Rv); - -} - -//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... -void ChiNu::BuildAnode(double Zpos){ - // Cu plate: 17 um - G4Tubs* Cu_plate_solid = new G4Tubs("Cu_plate",0,40*mm,0.5*ChiNu_NS::Cu_Thickness,0,360*deg); - G4Material* Cu_material = MaterialManager::getInstance()->GetMaterialFromLibrary("Cu"); - G4LogicalVolume* Cu_vol = new G4LogicalVolume(Cu_plate_solid, Cu_material,"logic_Cu",0,0,0); - Cu_vol->SetVisAttributes(m_VisCu); - - // Kapton: 50 um - G4Tubs* Kapton_solid = new G4Tubs("Kapton",0,40*mm,0.5*ChiNu_NS::Kapton_Thickness,0,360*deg); - G4Material* Kapton_material = MaterialManager::getInstance()->GetMaterialFromLibrary("Kapton"); - G4LogicalVolume* Kapton_vol = new G4LogicalVolume(Kapton_solid, Kapton_material,"logic_Kapton",0,0,0); - Kapton_vol->SetVisAttributes(m_VisFCWall); - - G4RotationMatrix *Rv=new G4RotationMatrix(0,0,0); - G4ThreeVector Tv; - Tv.setX(0); Tv.setY(0); Tv.setZ(0); - - Tv.setZ(Zpos); - m_FissionChamberVolume->AddPlacedVolume(Kapton_vol, Tv, Rv); - Tv.setZ(Zpos-0.5*ChiNu_NS::Kapton_Thickness-0.5*ChiNu_NS::Cu_Thickness); - m_FissionChamberVolume->AddPlacedVolume(Cu_vol, Tv, Rv); - Tv.setZ(Zpos+0.5*ChiNu_NS::Kapton_Thickness+0.5*ChiNu_NS::Cu_Thickness); - m_FissionChamberVolume->AddPlacedVolume(Cu_vol, Tv, Rv); - - -} + //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... G4AssemblyVolume* ChiNu::BuildDetector(){ if(!m_CylindricalDetector){ @@ -422,7 +243,7 @@ void ChiNu::ReadConfiguration(NPL::InputParser parser){ void ChiNu::ConstructDetector(G4LogicalVolume* world){ G4Material* Air = MaterialManager::getInstance()->GetMaterialFromLibrary("Air"); - //world->SetMaterial(Air); + world->SetMaterial(Air); for (unsigned short i = 0 ; i < m_R.size() ; i++) { @@ -446,12 +267,6 @@ void ChiNu::ConstructDetector(G4LogicalVolume* world){ G4RotationMatrix* Rot = new G4RotationMatrix(u,v,w); BuildDetector()->MakeImprint(world,Det_pos,Rot,i); } - - G4RotationMatrix* Rot_FC = new G4RotationMatrix(0,0,0); - G4ThreeVector Pos_FC = G4ThreeVector(0,0,0) ; - //BuildFissionChamber()->MakeImprint(world,Pos_FC,Rot_FC,0); - - } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... // Add Detector branch to the EventTree. diff --git a/NPSimulation/Detectors/ChiNu/ChiNu.hh b/NPSimulation/Detectors/ChiNu/ChiNu.hh index 879f101ae910fae1745f4d8625a7f3f303830ade..521f0f578b952290125ae9f2cf10665ee4d5f540 100644 --- a/NPSimulation/Detectors/ChiNu/ChiNu.hh +++ b/NPSimulation/Detectors/ChiNu/ChiNu.hh @@ -57,9 +57,6 @@ class ChiNu : public NPS::VDetector{ G4AssemblyVolume* BuildDetector(); - G4AssemblyVolume* BuildFissionChamber(); - void BuildAnode(double PosZ); - void BuildCathode(double PosZ); private: G4LogicalVolume* m_CylindricalDetector; @@ -68,8 +65,6 @@ class ChiNu : public NPS::VDetector{ G4LogicalVolume* m_LeadShield; G4AssemblyVolume* m_AssemblyVolume; - G4LogicalVolume* m_FissionChamberWall; - G4AssemblyVolume* m_FissionChamberVolume; //////////////////////////////////////////////////// ////// Inherite from NPS::VDetector class ///////// diff --git a/NPSimulation/Detectors/FissionChamber/CMakeLists.txt b/NPSimulation/Detectors/FissionChamber/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..673ee162588914dfc1d0c81b1cde5574f336f5ce --- /dev/null +++ b/NPSimulation/Detectors/FissionChamber/CMakeLists.txt @@ -0,0 +1,2 @@ +add_library(NPSFissionChamber SHARED FissionChamber.cc) +target_link_libraries(NPSFissionChamber NPSCore ${ROOT_LIBRARIES} ${Geant4_LIBRARIES} ${NPLib_LIBRARIES} -lNPFissionChamber) diff --git a/NPSimulation/Detectors/FissionChamber/FissionChamber.cc b/NPSimulation/Detectors/FissionChamber/FissionChamber.cc new file mode 100644 index 0000000000000000000000000000000000000000..1f2b491c78ccc9a2f9652c1142492dcb467080d5 --- /dev/null +++ b/NPSimulation/Detectors/FissionChamber/FissionChamber.cc @@ -0,0 +1,457 @@ +/***************************************************************************** + * Copyright (C) 2009-2020 this file is part of the NPTool Project * + * * + * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * + * For the list of contributors see $NPTOOL/Licence/Contributors * + *****************************************************************************/ + +/***************************************************************************** + * Original Author: Pierre Morfouace contact address: pierre.morfouace2@cea.fr * + * * + * Creation Date : September 2020 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class describe FissionChamber simulation * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + *****************************************************************************/ + +// C++ headers +#include <sstream> +#include <cmath> +#include <limits> +//G4 Geometry object +#include "G4Tubs.hh" +#include "G4Box.hh" + +//G4 sensitive +#include "G4SDManager.hh" +#include "G4MultiFunctionalDetector.hh" + +//G4 various object +#include "G4Material.hh" +#include "G4Transform3D.hh" +#include "G4PVPlacement.hh" +#include "G4VisAttributes.hh" +#include "G4Colour.hh" +#include "G4SubtractionSolid.hh" + +// NPTool header +#include "FissionChamber.hh" +#include "CalorimeterScorers.hh" +#include "InteractionScorers.hh" +#include "RootOutput.h" +#include "MaterialManager.hh" +#include "NPSDetectorFactory.hh" +#include "NPOptionManager.h" +#include "NPSHitsMap.hh" +// CLHEP header +#include "CLHEP/Random/RandGauss.h" + +using namespace std; +using namespace CLHEP; + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +namespace FissionChamber_NS{ + // Energy and time Resolution + const double EnergyThreshold = 0.1*MeV; + const double ResoTime = 4.5*ns ; + const double ResoEnergy = 1.0*MeV ; + const double Radius = 50*mm ; + const double Width = 100*mm ; + const double Thickness = 300*mm ; + const string Material = "BC400"; + + // Fission Chamber + const double Cu_Thickness = 17*micrometer; + // const double Al_Thickness = 12*micrometer; + const double Kapton_Thickness = 50*micrometer; + + + +} +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// FissionChamber Specific Method +FissionChamber::FissionChamber(){ + m_Event = new TFissionChamberData() ; + m_FissionChamberScorer = 0; + m_AssemblyVolume = 0; + m_FissionChamberVolume = 0; + + + // RGB Color + Transparency + m_VisFCWall = new G4VisAttributes(G4Colour(0.1,0.5,0.7,1)); + m_VisAl = new G4VisAttributes(G4Colour(0.839,0.803,0.803,1)); + m_VisTi = new G4VisAttributes(G4Colour(0.776,0.662,0.662,0.5)); + m_VisGas = new G4VisAttributes(G4Colour(0.576,0.662,0.662,0.3)); + m_VisCu = new G4VisAttributes(G4Colour(0.70, 0.40, 0. ,1)); + m_VisRogers4003C = new G4VisAttributes(G4Colour(0.60, 0.60, 0.2 ,1)); + + m_GasMaterial = "CF4"; + m_Pressure = 1.*bar; + +} + +FissionChamber::~FissionChamber(){ +} +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +void FissionChamber::AddDetector(G4ThreeVector POS){ + // Convert the POS value to R theta Phi as Spherical coordinate is easier in G4 + m_R.push_back(POS.mag()); + m_Theta.push_back(POS.theta()); + m_Phi.push_back(POS.phi()); +} + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +void FissionChamber::AddDetector(double R, double Theta, double Phi){ + m_R.push_back(R); + m_Theta.push_back(Theta); + m_Phi.push_back(Phi); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// Virtual Method of NPS::VDetector class + +// Read stream at Configfile to pick-up parameters of detector (Position,...) +// Called in DetecorConstruction::ReadDetextorConfiguration Method +void FissionChamber::ReadConfiguration(NPL::InputParser parser){ + vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("FissionChamber"); + if(NPOptionManager::getInstance()->GetVerboseLevel()) + cout << "//// " << blocks.size() << " detectors found " << endl; + + vector<string> cart = {"POS","GasMaterial","Pressure"}; + vector<string> sphe = {"R","Theta","Phi","GasMaterial","Pressure"}; + + for(unsigned int i = 0 ; i < blocks.size() ; i++){ + if(blocks[i]->HasTokenList(cart)){ + if(NPOptionManager::getInstance()->GetVerboseLevel()) + cout << endl << "//// FissionChamber " << i+1 << endl; + + G4ThreeVector Pos = NPS::ConvertVector(blocks[i]->GetTVector3("POS","mm")); + string GasMaterial = blocks[i]->GetString("GasMaterial"); + double Pressure = blocks[i]->GetDouble("Pressure","bar"); + AddDetector(Pos); + m_GasMaterial = GasMaterial; + m_Pressure = Pressure; + } + else if(blocks[i]->HasTokenList(sphe)){ + if(NPOptionManager::getInstance()->GetVerboseLevel()) + cout << endl << "//// FissionChamber " << 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 GasMaterial = blocks[i]->GetString("GasMaterial"); + double Pressure = blocks[i]->GetDouble("Pressure","bar"); + AddDetector(R,Theta,Phi); + m_GasMaterial = GasMaterial; + m_Pressure = Pressure; + } + else{ + cout << "ERROR: check your input file formatting " << endl; + exit(1); + } + } +} + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// Construct detector and inialise sensitive part. +// Called After DetecorConstruction::AddDetector Method +void FissionChamber::ConstructDetector(G4LogicalVolume* world){ + for (unsigned short i = 0 ; i < m_R.size() ; i++) { + + G4double wX = m_R[i] * sin(m_Theta[i] ) * cos(m_Phi[i] ) ; + G4double wY = m_R[i] * sin(m_Theta[i] ) * sin(m_Phi[i] ) ; + G4double wZ = m_R[i] * cos(m_Theta[i] ) ; + G4ThreeVector Det_pos = G4ThreeVector(wX, wY, wZ) ; + // So the face of the detector is at R instead of the middle + Det_pos+=Det_pos.unit()*FissionChamber_NS::Thickness*0.5; + // Building Detector reference frame + G4double ii = cos(m_Theta[i]) * cos(m_Phi[i]); + G4double jj = cos(m_Theta[i]) * sin(m_Phi[i]); + G4double kk = -sin(m_Theta[i]); + G4ThreeVector Y(ii,jj,kk); + G4ThreeVector w = Det_pos.unit(); + G4ThreeVector u = w.cross(Y); + G4ThreeVector v = w.cross(u); + v = v.unit(); + u = u.unit(); + + G4RotationMatrix* Rot = new G4RotationMatrix(0,0,0); + BuildFissionChamber()->MakeImprint(world, Det_pos, Rot, i); + } +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +G4AssemblyVolume* FissionChamber::BuildFissionChamber(){ + m_FissionChamberVolume = new G4AssemblyVolume(); + + G4RotationMatrix *Rv=new G4RotationMatrix(0,0,0); + G4ThreeVector Tv; + Tv.setX(0); Tv.setY(0); Tv.setZ(0); + + // Gas Volume // + double gas_height = 15.8*cm; + double gas_width = 14.8*cm; + double gas_length = 29.8*cm; + G4Box* gas_solid = new G4Box("Gas_solid", 0.5*gas_width, 0.5*gas_height, 0.5*gas_length); + + G4Material* gas_material = MaterialManager::getInstance()->GetGasFromLibrary(m_GasMaterial, m_Pressure, 300*kelvin); + G4LogicalVolume* gas_volume = new G4LogicalVolume(gas_solid, gas_material, "gas_logic", 0, 0, 0); + gas_volume->SetSensitiveDetector(m_FissionChamberScorer); + //m_VisGas->SetForceWireframe(true); + gas_volume->SetVisAttributes(m_VisGas); + m_FissionChamberVolume->AddPlacedVolume(gas_volume, Tv, Rv); + + + // Bottom PCB plate // + double PCB_width = 18.*cm; + double PCB_length = 33.*cm; + double PCB_Rogers_height = 1.6*mm; + double PCB_Cu_height = 6*35.*um; + double PCB_PosY = -8.5*cm; + // Cu layers + G4Box* PCB_Cu_solid = new G4Box("PCB_Cu_solid",0.5*PCB_width,0.5*PCB_Cu_height,0.5*PCB_length); + G4Material* Cu_material = MaterialManager::getInstance()->GetMaterialFromLibrary("Cu"); + G4LogicalVolume* PCB_Cu_vol = new G4LogicalVolume(PCB_Cu_solid, Cu_material,"PCB_Cu_logic",0,0,0); + PCB_Cu_vol->SetVisAttributes(m_VisCu); + Tv.setY(PCB_PosY); + m_FissionChamberVolume->AddPlacedVolume(PCB_Cu_vol, Tv, Rv); + + // Rogers 4003C layers + G4Box* PCB_Rogers_solid = new G4Box("PCB_Rogers_solid",0.5*PCB_width,0.5*PCB_Rogers_height,0.5*PCB_length); + G4Material* Rogers_material = MaterialManager::getInstance()->GetMaterialFromLibrary("Rogers4003C"); + G4LogicalVolume* PCB_Rogers_vol = new G4LogicalVolume(PCB_Rogers_solid, Rogers_material,"PCB_Rogers_logic",0,0,0); + PCB_Rogers_vol->SetVisAttributes(m_VisRogers4003C); + Tv.setY(PCB_PosY + 0.5*PCB_Cu_height + 0.5*PCB_Rogers_height); + m_FissionChamberVolume->AddPlacedVolume(PCB_Rogers_vol, Tv, Rv); + + // Al frame // + double frame1_width = 18.*cm; + double frame1_height = 5.*mm; + double frame1_length = 33.*cm; + double frame2_width = 15.2*cm; + double frame2_height = 5.1*mm; + double frame2_length = 30.2*cm; + + G4Box* frame1 = new G4Box("frame1", 0.5*frame1_width, 0.5*frame1_height, 0.5*frame1_length); + G4Box* frame2 = new G4Box("frame2", 0.5*frame2_width, 0.5*frame2_height, 0.5*frame2_length); + G4VSolid* Al_frame = (G4VSolid*) new G4SubtractionSolid("Al_frame",frame1,frame2,0,G4ThreeVector(0,0,0)); + G4Material* Al_material = MaterialManager::getInstance()->GetMaterialFromLibrary("Al"); + G4LogicalVolume* Al_frame_vol = new G4LogicalVolume(Al_frame,Al_material,"Al_frame_logic",0,0,0); + Al_frame_vol->SetVisAttributes(m_VisAl); + Tv.setY(PCB_PosY+ 0.5*PCB_Cu_height + PCB_Rogers_height + 0.5*frame1_height); + m_FissionChamberVolume->AddPlacedVolume(Al_frame_vol, Tv, Rv); + Tv.setY(PCB_PosY- 0.5*PCB_Cu_height - 0.5*frame1_height); + m_FissionChamberVolume->AddPlacedVolume(Al_frame_vol, Tv, Rv); + + double box1_width = 15.*cm; + double box1_height = 16.*cm; + double box1_length = 30.*cm; + double box2_width = 14.8*cm; + double box2_height = 15.8*cm; + double box2_length = 29.8*cm; + double box3_width = 12.5*cm; + double box3_height = 11.7*cm; + double box3_length = 30.1*cm; + double box4_width = 15.1*cm; + double box4_height = 11.8*cm; + double box4_length = 27.4*cm; + double box5_width = 12.4*cm; + double box5_height = 16.1*cm; + double box5_length = 27.4*cm; + + G4Box* box1 = new G4Box("box1", 0.5*box1_width, 0.5*box1_height, 0.5*box1_length); + G4Box* box2 = new G4Box("box2", 0.5*box2_width, 0.5*box2_height, 0.5*box2_length); + G4Box* box3 = new G4Box("box3", 0.5*box3_width, 0.5*box3_height, 0.5*box3_length); + G4Box* box4 = new G4Box("box4", 0.5*box4_width, 0.5*box4_height, 0.5*box4_length); + G4Box* box5 = new G4Box("box5", 0.5*box5_width, 0.5*box5_height, 0.5*box5_length); + + G4VSolid* box_int1 = (G4VSolid*) new G4SubtractionSolid("box_int1",box1,box2,0,G4ThreeVector(0,0,0)); + G4VSolid* box_int2 = (G4VSolid*) new G4SubtractionSolid("box_int2",box_int1,box3,0,G4ThreeVector(0,0,0)); + G4VSolid* box_int3 = (G4VSolid*) new G4SubtractionSolid("box_int3",box_int2,box4,0,G4ThreeVector(0,0,0)); + G4VSolid* box_int4 = (G4VSolid*) new G4SubtractionSolid("box_int4",box_int3,box5,0,G4ThreeVector(0,0,0)); + + G4LogicalVolume* full_box_vol = new G4LogicalVolume(box_int4, Al_material, "full_box_logic", 0,0,0); + full_box_vol->SetVisAttributes(m_VisAl); + Tv.setY(0); + m_FissionChamberVolume->AddPlacedVolume(full_box_vol, Tv, Rv); + + // Ti foils // + double foil1_width = 13*cm; + double foil1_length = 29*cm; + double foil1_thickness = 100*um; + double foil2_width = 13*cm; + double foil2_height = 14*cm; + double foil2_thickness = 50*um; + + G4Material* Ti_material = MaterialManager::getInstance()->GetMaterialFromLibrary("Ti"); + G4Box* foil1_solid = new G4Box("foil1", 0.5*foil1_width, 0.5*foil1_thickness, 0.5*foil1_length); + G4LogicalVolume* foil1_vol = new G4LogicalVolume(foil1_solid, Ti_material, "foil1_logic", 0, 0, 0); + foil1_vol->SetVisAttributes(m_VisTi); + Tv.setY(0.5*box2_height); + m_FissionChamberVolume->AddPlacedVolume(foil1_vol, Tv, Rv); + Tv.setY(0); + Tv.setX(-0.5*box2_width); + Rv->rotateZ(90*deg); + m_FissionChamberVolume->AddPlacedVolume(foil1_vol, Tv, Rv); + Tv.setX(0.5*box2_width); + m_FissionChamberVolume->AddPlacedVolume(foil1_vol, Tv, Rv); + + G4Box* foil2_solid = new G4Box("foil2", 0.5*foil2_width, 0.5*foil2_height, 0.5*foil2_thickness); + G4LogicalVolume* foil2_vol = new G4LogicalVolume(foil2_solid, Ti_material, "foil2_logic", 0, 0, 0); + foil2_vol->SetVisAttributes(m_VisTi); + Tv.setX(0);Tv.setY(0);Tv.setZ(-0.5*box2_length); + m_FissionChamberVolume->AddPlacedVolume(foil2_vol, Tv, Rv); + Tv.setZ(0.5*box2_length); + m_FissionChamberVolume->AddPlacedVolume(foil2_vol, Tv, Rv); + + // Cathode and Anode // + BuildCathode(-27.5); + double origine_anode = -25*mm; + double origine_cathode = -22.5*mm; + for(int i=0; i<11; i++){ + BuildAnode(origine_anode+i*5*mm); + } + for(int i=0; i<10; i++){ + BuildCathode(origine_cathode+i*5*mm); + } + + BuildCathode(27.5); + + + return m_FissionChamberVolume; +} +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +void FissionChamber::BuildCathode(double Zpos){ + // Al plate: 12 um + G4Tubs* Al_plate_solid = new G4Tubs("Al_plate",0,40*mm,12*micrometer,0,360*deg); + G4Material* Al_material = MaterialManager::getInstance()->GetMaterialFromLibrary("Al"); + G4LogicalVolume* Al_vol = new G4LogicalVolume(Al_plate_solid, Al_material,"logic_Al",0,0,0); + Al_vol->SetVisAttributes(m_VisAl); + + G4RotationMatrix *Rv=new G4RotationMatrix(0,0,0); + G4ThreeVector Tv; + Tv.setX(0); Tv.setY(0); Tv.setZ(0); + + Tv.setZ(Zpos); + m_FissionChamberVolume->AddPlacedVolume(Al_vol, Tv, Rv); + +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +void FissionChamber::BuildAnode(double Zpos){ + // Cu plate: 17 um + G4Tubs* Cu_plate_solid = new G4Tubs("Cu_plate",0,40*mm,0.5*FissionChamber_NS::Cu_Thickness,0,360*deg); + G4Material* Cu_material = MaterialManager::getInstance()->GetMaterialFromLibrary("Cu"); + G4LogicalVolume* Cu_vol = new G4LogicalVolume(Cu_plate_solid, Cu_material,"logic_Cu",0,0,0); + Cu_vol->SetVisAttributes(m_VisCu); + + // Kapton: 50 um + G4Tubs* Kapton_solid = new G4Tubs("Kapton",0,40*mm,0.5*FissionChamber_NS::Kapton_Thickness,0,360*deg); + G4Material* Kapton_material = MaterialManager::getInstance()->GetMaterialFromLibrary("Kapton"); + G4LogicalVolume* Kapton_vol = new G4LogicalVolume(Kapton_solid, Kapton_material,"logic_Kapton",0,0,0); + Kapton_vol->SetVisAttributes(m_VisFCWall); + + G4RotationMatrix *Rv=new G4RotationMatrix(0,0,0); + G4ThreeVector Tv; + Tv.setX(0); Tv.setY(0); Tv.setZ(0); + + Tv.setZ(Zpos); + m_FissionChamberVolume->AddPlacedVolume(Kapton_vol, Tv, Rv); + Tv.setZ(Zpos-0.5*FissionChamber_NS::Kapton_Thickness-0.5*FissionChamber_NS::Cu_Thickness); + m_FissionChamberVolume->AddPlacedVolume(Cu_vol, Tv, Rv); + Tv.setZ(Zpos+0.5*FissionChamber_NS::Kapton_Thickness+0.5*FissionChamber_NS::Cu_Thickness); + m_FissionChamberVolume->AddPlacedVolume(Cu_vol, Tv, Rv); + + +} +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// Add Detector branch to the EventTree. +// Called After DetecorConstruction::AddDetector Method +void FissionChamber::InitializeRootOutput(){ + RootOutput *pAnalysis = RootOutput::getInstance(); + TTree *pTree = pAnalysis->GetTree(); + if(!pTree->FindBranch("FissionChamber")){ + pTree->Branch("FissionChamber", "TFissionChamberData", &m_Event) ; + } + pTree->SetBranchAddress("FissionChamber", &m_Event) ; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// Read sensitive part and fill the Root tree. +// Called at in the EventAction::EndOfEventAvtion +void FissionChamber::ReadSensitive(const G4Event* ){ + m_Event->Clear(); + + /////////// + // Calorimeter scorer + CalorimeterScorers::PS_Calorimeter* Scorer= (CalorimeterScorers::PS_Calorimeter*) m_FissionChamberScorer->GetPrimitive(0); + + unsigned int size = Scorer->GetMult(); + for(unsigned int i = 0 ; i < size ; i++){ + vector<unsigned int> level = Scorer->GetLevel(i); + double Energy = RandGauss::shoot(Scorer->GetEnergy(i),FissionChamber_NS::ResoEnergy); + if(Energy>FissionChamber_NS::EnergyThreshold){ + double Time = RandGauss::shoot(Scorer->GetTime(i),FissionChamber_NS::ResoTime); + int DetectorNbr = level[0]; + m_Event->SetEnergy(DetectorNbr,Energy); + m_Event->SetTime(DetectorNbr,Time); + } + } +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//////////////////////////////////////////////////////////////// +void FissionChamber::InitializeScorers() { + // This check is necessary in case the geometry is reloaded + bool already_exist = false; + m_FissionChamberScorer = CheckScorer("FissionChamberScorer",already_exist) ; + + if(already_exist) + return ; + + // Otherwise the scorer is initialised + vector<int> level; level.push_back(0); + G4VPrimitiveScorer* Calorimeter= new CalorimeterScorers::PS_Calorimeter("Calorimeter",level, 0) ; + G4VPrimitiveScorer* Interaction= new InteractionScorers::PS_Interactions("Interaction",ms_InterCoord, 0) ; + //and register it to the multifunctionnal detector + m_FissionChamberScorer->RegisterPrimitive(Calorimeter); + m_FissionChamberScorer->RegisterPrimitive(Interaction); + G4SDManager::GetSDMpointer()->AddNewDetector(m_FissionChamberScorer) ; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//////////////////////////////////////////////////////////////////////////////// +// Construct Method to be pass to the DetectorFactory // +//////////////////////////////////////////////////////////////////////////////// +NPS::VDetector* FissionChamber::Construct(){ + return (NPS::VDetector*) new FissionChamber(); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//////////////////////////////////////////////////////////////////////////////// +// Registering the construct method to the factory // +//////////////////////////////////////////////////////////////////////////////// +extern"C" { + class proxy_nps_FissionChamber{ + public: + proxy_nps_FissionChamber(){ + NPS::DetectorFactory::getInstance()->AddToken("FissionChamber","FissionChamber"); + NPS::DetectorFactory::getInstance()->AddDetector("FissionChamber",FissionChamber::Construct); + } + }; + + proxy_nps_FissionChamber p_nps_FissionChamber; +} diff --git a/NPSimulation/Detectors/FissionChamber/FissionChamber.hh b/NPSimulation/Detectors/FissionChamber/FissionChamber.hh new file mode 100644 index 0000000000000000000000000000000000000000..713d0015e1a49ee1ce00d556a24d874d78bd1c7f --- /dev/null +++ b/NPSimulation/Detectors/FissionChamber/FissionChamber.hh @@ -0,0 +1,125 @@ +#ifndef FissionChamber_h +#define FissionChamber_h 1 +/***************************************************************************** + * Copyright (C) 2009-2020 this file is part of the NPTool Project * + * * + * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * + * For the list of contributors see $NPTOOL/Licence/Contributors * + *****************************************************************************/ + +/***************************************************************************** + * Original Author: Pierre Morfouace contact address: pierre.morfouace2@cea.fr * + * * + * Creation Date : September 2020 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class describe FissionChamber simulation * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + *****************************************************************************/ + +// C++ header +#include <string> +#include <vector> +using namespace std; + +// G4 headers +#include "G4ThreeVector.hh" +#include "G4RotationMatrix.hh" +#include "G4LogicalVolume.hh" +#include "G4MultiFunctionalDetector.hh" +#include "G4AssemblyVolume.hh" + +// NPTool header +#include "NPSVDetector.hh" +#include "TFissionChamberData.h" +#include "NPInputParser.h" + +class FissionChamber : public NPS::VDetector{ + //////////////////////////////////////////////////// + /////// Default Constructor and Destructor ///////// + //////////////////////////////////////////////////// + public: + FissionChamber() ; + virtual ~FissionChamber() ; + + //////////////////////////////////////////////////// + /////// Specific Function of this Class /////////// + //////////////////////////////////////////////////// + public: + // Cartesian + void AddDetector(G4ThreeVector POS); + // Spherical + void AddDetector(double R,double Theta,double Phi); + + + G4AssemblyVolume* BuildDetector(); + G4AssemblyVolume* BuildFissionChamber(); + void BuildAnode(double PosZ); + void BuildCathode(double PosZ); + + private: + G4LogicalVolume* m_Detector; + G4AssemblyVolume* m_AssemblyVolume; + + G4AssemblyVolume* m_FissionChamberVolume; + + //////////////////////////////////////////////////// + ////// Inherite from NPS::VDetector class ///////// + //////////////////////////////////////////////////// + public: + // Read stream at Configfile to pick-up parameters of detector (Position,...) + // Called in DetecorConstruction::ReadDetextorConfiguration Method + void ReadConfiguration(NPL::InputParser) ; + + // Construct detector and inialise sensitive part. + // Called After DetecorConstruction::AddDetector Method + void ConstructDetector(G4LogicalVolume* world) ; + + // Add Detector branch to the EventTree. + // Called After DetecorConstruction::AddDetector Method + void InitializeRootOutput() ; + + // Read sensitive part and fill the Root tree. + // Called at in the EventAction::EndOfEventAvtion + void ReadSensitive(const G4Event* event) ; + + public: // Scorer + void InitializeScorers() ; + + // Associated Scorer + G4MultiFunctionalDetector* m_FissionChamberScorer ; + //////////////////////////////////////////////////// + ///////////Event class to store Data//////////////// + //////////////////////////////////////////////////// + private: + TFissionChamberData* m_Event ; + + //////////////////////////////////////////////////// + ///////////////Private intern Data////////////////// + //////////////////////////////////////////////////// + private: // Geometry + // Detector Coordinate + vector<double> m_R; + vector<double> m_Theta; + vector<double> m_Phi; + + string m_GasMaterial; + double m_Pressure; + + // Visualisation Attribute + G4VisAttributes* m_VisFCWall; + G4VisAttributes* m_VisAl; + G4VisAttributes* m_VisCu; + G4VisAttributes* m_VisGas; + G4VisAttributes* m_VisTi; + G4VisAttributes* m_VisRogers4003C; + + // Needed for dynamic loading of the library + public: + static NPS::VDetector* Construct(); +}; +#endif diff --git a/NPSimulation/Detectors/PISTA/PISTA.cc b/NPSimulation/Detectors/PISTA/PISTA.cc index 73bec673c2b2113b8126421aa65a5d8f0f20c547..34c3f4dc24bdd70f5b9595cd9ebb02f9b460ebb4 100644 --- a/NPSimulation/Detectors/PISTA/PISTA.cc +++ b/NPSimulation/Detectors/PISTA/PISTA.cc @@ -121,7 +121,7 @@ G4LogicalVolume* PISTA::BuildTrapezoidDetector(){ "PISTA", 0,0,0); - G4VisAttributes* TrapezoidVisAtt = new G4VisAttributes(G4Colour(0.2, 0.80, 0.50)); + G4VisAttributes* TrapezoidVisAtt = new G4VisAttributes(G4Colour(0,0,0,0.5)); TrapezoidVisAtt->SetForceWireframe(true); logicTrapezoid->SetVisAttributes(TrapezoidVisAtt); @@ -147,7 +147,8 @@ G4LogicalVolume* PISTA::BuildTrapezoidDetector(){ logicFirstStage->SetSensitiveDetector(m_FirstStageScorer); // Visualisation of First Stage strips - G4VisAttributes* FirstStageVisAtt = new G4VisAttributes(G4Colour(0.2,0.8,0.5)); + //G4VisAttributes* FirstStageVisAtt = new G4VisAttributes(G4Colour(0.2,0.8,0.5)); + G4VisAttributes* FirstStageVisAtt = new G4VisAttributes(G4Colour(0.5,0.5,0.55,0.8)); logicFirstStage->SetVisAttributes(FirstStageVisAtt); ////// @@ -173,7 +174,7 @@ G4LogicalVolume* PISTA::BuildTrapezoidDetector(){ logicSecondStage->SetSensitiveDetector(m_SecondStageScorer); // Visualisation of Second Stage strips - G4VisAttributes* SecondStageVisAtt = new G4VisAttributes(G4Colour(0.2,0.8,0.5)); + G4VisAttributes* SecondStageVisAtt = new G4VisAttributes(G4Colour(0.5,0.5,0.5)); logicSecondStage->SetVisAttributes(SecondStageVisAtt); diff --git a/NPSimulation/Detectors/Plastic/Plastic.cc b/NPSimulation/Detectors/Plastic/Plastic.cc index 94f2013d2edff0bc81ad6e3d9e67a98a1f837122..fe2a6008e813d1f645e329a0cce7dd81ca5c61f3 100644 --- a/NPSimulation/Detectors/Plastic/Plastic.cc +++ b/NPSimulation/Detectors/Plastic/Plastic.cc @@ -42,6 +42,7 @@ // NPTool header #include "Plastic.hh" +#include "CalorimeterScorers.hh" #include "ObsoleteGeneralScorers.hh" #include "RootOutput.h" #include "MaterialManager.hh" @@ -59,7 +60,8 @@ using namespace CLHEP; namespace PLASTIC{ // Energy and time Resolution const G4double ResoTime = 1. ;// Resolution in ns // - const G4double ResoEnergy = 0.1 ;// Resolution in % + //const G4double ResoEnergy = 0.1 ;// Resolution in % + const G4double ResoEnergy = 1*keV; // Resolution } using namespace PLASTIC ; @@ -341,6 +343,17 @@ void Plastic::ReadSensitive(const G4Event* event){ G4String DetectorNumber; m_Event->Clear(); + CalorimeterScorers::PS_Calorimeter* Scorer = (CalorimeterScorers::PS_Calorimeter*) m_PlasticScorer->GetPrimitive(0); + + unsigned int size = Scorer->GetMult(); + for(unsigned int i=0; i<size; i++){ + double Energy = RandGauss::shoot(Scorer->GetEnergy(i), ResoEnergy); + double Time = RandGauss::shoot(Scorer->GetTime(i), ResoTime); + int DetectorNbr = Scorer->GetLevel(i)[0]; + m_Event->SetEnergyAndTime(DetectorNbr,Energy,Time); + } + + /* ////////////////////////////////////////////////////////////////////////////////////// //////////////////////// Used to Read Event Map of detector ////////////////////////// ////////////////////////////////////////////////////////////////////////////////////// @@ -386,45 +399,46 @@ void Plastic::ReadSensitive(const G4Event* event){ // Loop on Plastic Number for (G4int l = 0 ; l < sizeN ; l++) { - G4int N = *(DetectorNumber_itr->second); - G4int NTrackID = DetectorNumber_itr->first - N; - if (N > 0) { - det.push_back(N); - // Energy - Energy_itr = EnergyHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < sizeE ; h++) { - G4int ETrackID = Energy_itr->first - N; - G4double E = *(Energy_itr->second); - if (ETrackID == NTrackID) { - energy.push_back(RandGauss::shoot(E, E*ResoEnergy/100./2.35)) ; - } - Energy_itr++; - } + G4int N = *(DetectorNumber_itr->second); + G4int NTrackID = DetectorNumber_itr->first - N; + if (N > 0) { + det.push_back(N); + // Energy + Energy_itr = EnergyHitMap->GetMap()->begin(); + for (G4int h = 0 ; h < sizeE ; h++) { + G4int ETrackID = Energy_itr->first - N; + G4double E = *(Energy_itr->second); + if (ETrackID == NTrackID) { + //energy.push_back(RandGauss::shoot(E, E*ResoEnergy/100./2.35)) ; + energy.push_back(RandGauss::shoot(E, ResoEnergy)) ; + } + Energy_itr++; + } - // Time - Time_itr = TimeHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < sizeT ; h++) { - G4int TTrackID = Time_itr->first - N; - G4double T = *(Time_itr->second); - if (TTrackID == NTrackID) { - time.push_back(RandGauss::shoot(T, ResoTime)); - } - Time_itr++; - } - } - DetectorNumber_itr++; - } - unsigned int size=energy.size(); - for(unsigned int i = 0 ; i < size ; i++){ - m_Event->SetEnergyAndTime(det[i],energy[i],time[i]); - } + // Time + Time_itr = TimeHitMap->GetMap()->begin(); + for (G4int h = 0 ; h < sizeT ; h++) { + G4int TTrackID = Time_itr->first - N; + G4double T = *(Time_itr->second); + if (TTrackID == NTrackID) { + time.push_back(RandGauss::shoot(T, ResoTime)); + } + Time_itr++; +} +} +DetectorNumber_itr++; +} +unsigned int size=energy.size(); +for(unsigned int i = 0 ; i < size ; i++){ + m_Event->SetEnergyAndTime(det[i],energy[i],time[i]); +} - // clear map for next event - TimeHitMap->clear(); - DetectorNumberHitMap->clear(); - EnergyHitMap->clear(); +// clear map for next event +TimeHitMap->clear(); +DetectorNumberHitMap->clear(); +EnergyHitMap->clear();*/ } @@ -434,13 +448,19 @@ void Plastic::InitializeScorers() { m_PlasticScorer = CheckScorer("PlasticScorer",already_exist) ; if(already_exist) return ; - G4VPrimitiveScorer* DetNbr = new PSDetectorNumber("PlasticNumber","Plastic", 0); - G4VPrimitiveScorer* Energy = new PSEnergy("Energy","Plastic", 0); - G4VPrimitiveScorer* Time = new PSTOF("Time","Plastic", 0); + + vector<int> level; level.push_back(1); + G4VPrimitiveScorer* Calorimeter = new CalorimeterScorers::PS_Calorimeter("Plastic",level, 0); + m_PlasticScorer->RegisterPrimitive(Calorimeter); + + /*G4VPrimitiveScorer* DetNbr = new PSDetectorNumber("PlasticNumber","Plastic", 0); + G4VPrimitiveScorer* Energy = new PSEnergy("Energy","Plastic", 0); + G4VPrimitiveScorer* Time = new PSTOF("Time","Plastic", 0); //and register it to the multifunctionnal detector m_PlasticScorer->RegisterPrimitive(DetNbr); m_PlasticScorer->RegisterPrimitive(Energy); - m_PlasticScorer->RegisterPrimitive(Time); + m_PlasticScorer->RegisterPrimitive(Time);*/ + G4SDManager::GetSDMpointer()->AddNewDetector(m_PlasticScorer); } //////////////////////////////////////////////////////////////// diff --git a/Projects/ChiNu/Analysis.cxx b/Projects/ChiNu/Analysis.cxx index 05863c3778f662b54f55076ceb5fedc7fa8a559a..85e60dd976c99e5bc9e2606454656024fa01a69c 100644 --- a/Projects/ChiNu/Analysis.cxx +++ b/Projects/ChiNu/Analysis.cxx @@ -82,16 +82,15 @@ void Analysis::TreatEvent(){ double DeltaTheta = atan(89.0/Rdet); - double random_ThetaLab = ra.Uniform(init_ThetaLab-DeltaTheta, init_ThetaLab+DeltaTheta); + double exp_ThetaLab = m_ChiNu->GetVectorDetectorPosition(m_ChiNu->DetectorNumber[i]).Theta(); + double random_ThetaLab = ra.Uniform(exp_ThetaLab-DeltaTheta, exp_ThetaLab+DeltaTheta); double dEx = my_Reaction->ReconstructRelativistic(Elab[i], random_ThetaLab); ThetaLab.push_back(random_ThetaLab/deg); + //ThetaLab.push_back(exp_ThetaLab/deg); Ex.push_back(dEx); } } - - - } /////////////////////////////////////////////////////////////////////////////// diff --git a/Projects/ChiNu/ChiNu.detector b/Projects/ChiNu/ChiNu.detector deleted file mode 100644 index e3dc317cfd3c66051de7c31a390462c5bbfef330..0000000000000000000000000000000000000000 --- a/Projects/ChiNu/ChiNu.detector +++ /dev/null @@ -1,353 +0,0 @@ -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -%Target -% THICKNESS= 10 micrometer -% RADIUS= 20 mm -% MATERIAL= CD2 -% ANGLE= 0 deg -% X= 0 mm -% Y= 0 mm -% Z= 0 mm - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -%% LI -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 150 deg - PHI= 0 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 135 deg - PHI= 0 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 120 deg - PHI= 0 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 105 deg - PHI= 0 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 90 deg - PHI= 0 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 75 deg - PHI= 0 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 60 deg - PHI= 0 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 45 deg - PHI= 0 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 30 deg - PHI= 0 deg - LeadShield= 1 - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -%% LII -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 150 deg - PHI= 30 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 135 deg - PHI= 30 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 120 deg - PHI= 30 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 105 deg - PHI= 30 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 90 deg - PHI= 30 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 75 deg - PHI= 30 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 60 deg - PHI= 30 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 45 deg - PHI= 30 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 30 deg - PHI= 30 deg - LeadShield= 1 - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -%% LIII -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 150 deg - PHI= 60 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 135 deg - PHI= 60 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 120 deg - PHI= 60 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 105 deg - PHI= 60 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 90 deg - PHI= 60 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 75 deg - PHI= 60 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 60 deg - PHI= 60 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 45 deg - PHI= 60 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 30 deg - PHI= 60 deg - LeadShield= 1 - - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -%% RI -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 150 deg - PHI= 180 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 135 deg - PHI= 180 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 120 deg - PHI= 180 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 105 deg - PHI= 180 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 90 deg - PHI= 180 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 75 deg - PHI= 180 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 60 deg - PHI= 180 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 45 deg - PHI= 180 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 30 deg - PHI= 180 deg - LeadShield= 1 - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -%% RII -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 150 deg - PHI= 150 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 135 deg - PHI= 150 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 120 deg - PHI= 150 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 105 deg - PHI= 150 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 90 deg - PHI= 150 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 75 deg - PHI= 150 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 60 deg - PHI= 150 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 45 deg - PHI= 150 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 30 deg - PHI= 150 deg - LeadShield= 1 - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -%% RIII -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 150 deg - PHI= 120 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 135 deg - PHI= 120 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 120 deg - PHI= 120 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 105 deg - PHI= 120 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 90 deg - PHI= 120 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 75 deg - PHI= 120 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 60 deg - PHI= 120 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 45 deg - PHI= 120 deg - LeadShield= 1 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ChiNu - R= 1000 mm - THETA= 30 deg - PHI= 120 deg - LeadShield= 1 - diff --git a/Projects/ChiNu/chinu.detector b/Projects/ChiNu/chinu.detector new file mode 100644 index 0000000000000000000000000000000000000000..152b2dd9c346903cae81b13c8691653676de823f --- /dev/null +++ b/Projects/ChiNu/chinu.detector @@ -0,0 +1,20 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +Alias Theta + Action= Copy + Value= 30 45 60 75 90 105 120 135 150 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +Alias Phi + Action= Copy + Value= 0 30 60 120 150 180 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +ChiNu + R= 1000 mm + THETA= @Theta deg + PHI= @Phi deg + LeadShield= 1 + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +FissionChamber + POS= 0 0 0 mm + GasMaterial= CF4 + Pressure= 1 bar diff --git a/Projects/FissionChamber/Analysis.cxx b/Projects/FissionChamber/Analysis.cxx new file mode 100644 index 0000000000000000000000000000000000000000..7d5673e3d8df03c250cd40aa8921351cfb8210b6 --- /dev/null +++ b/Projects/FissionChamber/Analysis.cxx @@ -0,0 +1,68 @@ +/***************************************************************************** + * Copyright (C) 2009-2016 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: XAUTHORX contact address: XMAILX * + * * + * Creation Date : XMONTHX XYEARX * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class describe FissionChamber analysis project * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + *****************************************************************************/ + +#include<iostream> +using namespace std; +#include"Analysis.h" +#include"NPAnalysisFactory.h" +#include"NPDetectorManager.h" +//////////////////////////////////////////////////////////////////////////////// +Analysis::Analysis(){ +} +//////////////////////////////////////////////////////////////////////////////// +Analysis::~Analysis(){ +} + +//////////////////////////////////////////////////////////////////////////////// +void Analysis::Init(){ + FissionChamber= (TFissionChamberPhysicsPhysics*) m_DetectorManager->GetDetector("FissionChamber"); +} + +//////////////////////////////////////////////////////////////////////////////// +void Analysis::TreatEvent(){ +} + +//////////////////////////////////////////////////////////////////////////////// +void Analysis::End(){ +} + + +//////////////////////////////////////////////////////////////////////////////// +// Construct Method to be pass to the DetectorFactory // +//////////////////////////////////////////////////////////////////////////////// +NPL::VAnalysis* Analysis::Construct(){ + return (NPL::VAnalysis*) new Analysis(); +} + +//////////////////////////////////////////////////////////////////////////////// +// Registering the construct method to the factory // +//////////////////////////////////////////////////////////////////////////////// +extern "C"{ +class proxy{ + public: + proxy(){ + NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct); + } +}; + +proxy p; +} + diff --git a/Projects/FissionChamber/Analysis.h b/Projects/FissionChamber/Analysis.h new file mode 100644 index 0000000000000000000000000000000000000000..62f22889e258ebd9eb9d3f243142690d89dbe4e5 --- /dev/null +++ b/Projects/FissionChamber/Analysis.h @@ -0,0 +1,42 @@ +#ifndef Analysis_h +#define Analysis_h +/***************************************************************************** + * Copyright (C) 2009-2016 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: XAUTHORX contact address: XMAILX * + * * + * Creation Date : XMONTHX XYEARX * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class describe FissionChamber analysis project * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + *****************************************************************************/ + +#include"NPVAnalysis.h" +#include"TFissionChamberPhysics.h" +class Analysis: public NPL::VAnalysis{ + public: + Analysis(); + ~Analysis(); + + public: + void Init(); + void TreatEvent(); + void End(); + + static NPL::VAnalysis* Construct(); + + private: + TFissionChamberPhysics* FissionChamber; + +}; +#endif diff --git a/Projects/FissionChamber/CMakeLists.txt b/Projects/FissionChamber/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..22c74affdfc45019bdda2594f8439c52d4ab97ec --- /dev/null +++ b/Projects/FissionChamber/CMakeLists.txt @@ -0,0 +1,5 @@ +cmake_minimum_required (VERSION 2.8) +# Setting the policy to match Cmake version +cmake_policy(VERSION ${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION}) +# include the default NPAnalysis cmake file +include("../../NPLib/ressources/CMake/NPAnalysis.cmake") diff --git a/Projects/FissionChamber/FissionChamber.detector b/Projects/FissionChamber/FissionChamber.detector new file mode 100644 index 0000000000000000000000000000000000000000..cf8c9a2614d12120051da22fe34388b42de43834 --- /dev/null +++ b/Projects/FissionChamber/FissionChamber.detector @@ -0,0 +1,5 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +FissionChamber + POS= 0 0 0 mm + GasMaterial= CF4 + Pressure= 1 bar