From 1fe79370b6d12563e08a8731101c7549af9e7449 Mon Sep 17 00:00:00 2001
From: deserevi <deserevi@nptool>
Date: Tue, 12 Oct 2010 11:36:59 +0000
Subject: [PATCH] * Add support for annular shape in Gaspard + add reliable
scorer functionality for theta and phi strips in NPS + add theta and phi
strips position determination in NPL
* Status:
+ Analysis is working well, except that excitation energy, when using
the strips, is reconstructed -20 keV off.
---
.../gaspardV2Annular.detector | 3 +-
NPAnalysis/Gaspard/src/AnalysisNew.cc | 197 ++++++++++++
NPLib/GASPARD/GaspardTrackerAnnular.cxx | 303 ++++++++++++++++++
NPLib/GASPARD/GaspardTrackerAnnular.h | 77 +++++
NPLib/GASPARD/GaspardTrackerDummyShape.cxx | 28 +-
NPLib/GASPARD/TGaspardTrackerData.h | 4 +-
NPSimulation/src/GaspardScorers.cc | 7 +
NPSimulation/src/GaspardTrackerAnnular.cc | 38 +--
8 files changed, 609 insertions(+), 48 deletions(-)
create mode 100644 NPAnalysis/Gaspard/src/AnalysisNew.cc
create mode 100644 NPLib/GASPARD/GaspardTrackerAnnular.cxx
create mode 100644 NPLib/GASPARD/GaspardTrackerAnnular.h
diff --git a/Inputs/DetectorConfiguration/gaspardV2Annular.detector b/Inputs/DetectorConfiguration/gaspardV2Annular.detector
index fb193fcb3..39b7bf71f 100644
--- a/Inputs/DetectorConfiguration/gaspardV2Annular.detector
+++ b/Inputs/DetectorConfiguration/gaspardV2Annular.detector
@@ -22,6 +22,7 @@ Target
THICKNESS= 0.001
RADIUS= 7.5
MATERIAL= CD2
+ NBLAYERS= 0
X= 0
Y= 0
Z= 0
@@ -37,7 +38,7 @@ GPDAnnular
THIRDSTAGE= 1
VIS= all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1 Annular Front
-GPDAnnular
+%GPDAnnular
Z= 156.5
RMIN= 16
RMAX= 52
diff --git a/NPAnalysis/Gaspard/src/AnalysisNew.cc b/NPAnalysis/Gaspard/src/AnalysisNew.cc
new file mode 100644
index 000000000..b1c7a07a4
--- /dev/null
+++ b/NPAnalysis/Gaspard/src/AnalysisNew.cc
@@ -0,0 +1,197 @@
+#include "ObjectManager.hh"
+
+using namespace std;
+
+
+int main(int argc,char** argv)
+{
+ // test if number of arguments is correct
+ if (argc != 4) {
+ cout <<
+ "you need to specify both a Reaction file and a Detector file such as : Analysis myReaction.reaction myDetector.detector runToRead.run"
+ << endl;
+ return 0;
+ }
+
+ // get arguments
+ string reactionfileName = argv[1];
+ string detectorfileName = argv[2];
+ string runToReadfileName = argv[3];
+
+ // Instantiate RootInput and RootOutput singleton classes
+ RootInput:: getInstance(runToReadfileName);
+ RootOutput::getInstance("Analysis/Gaspard_AnalyzedData", "AnalyzedTree");
+
+ // Initialize the reaction
+ NPL::Reaction* myReaction = new Reaction();
+ myReaction->ReadConfigurationFile(reactionfileName);
+
+ // Initialize the detector
+ NPA::DetectorManager* myDetector = new DetectorManager;
+ myDetector->ReadConfigurationFile(detectorfileName);
+
+ // Calculate beam energy at target middle
+ // Get nominal beam energy
+ Double_t BeamEnergyNominal = myReaction->GetBeamEnergy() * MeV;
+ cout << BeamEnergyNominal << endl;
+ // Slow beam at target middle
+ Double_t BeamEnergy = BeamEnergyNominal - BeamTarget.Slow(BeamEnergyNominal, myDetector->GetTargetThickness()/2 * micrometer, 0);
+ cout << BeamEnergy << endl;
+ // Set energy beam at target middle
+ myReaction->SetBeamEnergy(BeamEnergy);
+
+ // Print target thickness
+ cout << myDetector->GetTargetThickness() << endl;
+
+ // Attach more branch to the output
+ double Ex = 0 ; double ExNoStrips = 0 ; double EE = 0 ; double TT = 0 ; double X = 0 ; double Y = 0 ; int det ;
+ RootOutput::getInstance()->GetTree()->Branch("ExcitationEnergy",&Ex,"Ex/D") ;
+ RootOutput::getInstance()->GetTree()->Branch("ExcitationEnergyNoStrips",&ExNoStrips,"ExNoStrips/D") ;
+ RootOutput::getInstance()->GetTree()->Branch("E",&EE,"EE/D") ;
+ RootOutput::getInstance()->GetTree()->Branch("A",&TT,"TT/D") ;
+ RootOutput::getInstance()->GetTree()->Branch("X",&X,"X/D") ;
+ RootOutput::getInstance()->GetTree()->Branch("Y",&Y,"Y/D") ;
+
+ // Get GaspardTracker pointer
+ GaspardTrackerNew* GPDTrack = (GaspardTrackerNew*) myDetector->m_Detector["GASPARD"];
+
+ // Get the input TChain and treat it
+ TChain* chain = RootInput:: getInstance() -> GetChain();
+
+ // Connect TInitialConditions branch
+ TInitialConditions *initCond = 0;
+ chain->SetBranchAddress("InitialConditions", &initCond);
+ chain->SetBranchStatus("InitialConditions", 1);
+
+ // Connect TInteractionCoordinates branch
+ TInteractionCoordinates *interCoord = 0;
+ chain->SetBranchAddress("InteractionCoordinates", &interCoord);
+ chain->SetBranchStatus("InteractionCoordinates", 0);
+
+ // Analysis is here!
+ int nentries = chain->GetEntries();
+ cout << "Number of entries to be analysed: " << nentries << endl;
+
+ // Default initialization
+ double XTarget = 0;
+ double YTarget = 0;
+ double BeamTheta = 0;
+ double BeamPhi = 0;
+
+ // random generator
+ TRandom3 *gene = new TRandom3();
+
+ // Loop on all events
+ for (int i = 0; i < nentries; i ++) {
+ if (i%10000 == 0 && i!=0) cout << "\r" << i << " analyzed events" << flush;
+ chain -> GetEntry(i);
+
+ // Treat Gaspard event
+ myDetector->ClearEventPhysics();
+ myDetector->BuildPhysicalEvent();
+
+ // Get total energy
+ double E = GPDTrack->GetEnergyDeposit();
+// cout << i << " " << E << endl;
+
+ // if there is a hit in the detector array, treat it.
+ double Theta, ThetaStrip, angle, ThetaCM;
+ double DetecX, DetecY, DetecZ;
+ double r;
+ TVector3 A;
+ if (E > -1000) {
+ // Get c.m. angle
+ ThetaCM = initCond->GetICEmittedAngleThetaCM(0) * deg;
+
+ // Get exact scattering angle from TInteractionCoordinates object
+// Theta = interCoord->GetDetectedAngleTheta(0) * deg;
+// cout << interCoord << endl;
+// interCoord->Dump();
+// cout << i << " mult: " << interCoord->GetDetectedMultiplicity() << endl;
+ DetecX = interCoord->GetDetectedPositionX(0);
+ DetecY = interCoord->GetDetectedPositionY(0);
+ DetecZ = interCoord->GetDetectedPositionZ(0);
+// cout << DetecX << " " << DetecY << " " << DetecZ << endl;
+ TVector3 Detec(DetecX, DetecY, DetecZ);
+
+ // Get interaction position in detector
+ // This takes into account the strips
+ A = GPDTrack->GetPositionOfInteraction();
+
+ // Get beam interaction coordinates on target (from initial condition)
+ XTarget = initCond->GetICPositionX(0);
+ YTarget = initCond->GetICPositionY(0);
+// cout << XTarget << " " << YTarget << endl;
+ BeamTheta = initCond->GetICIncidentAngleTheta(0)*deg;
+ BeamPhi = initCond->GetICIncidentAnglePhi(0)*deg;
+ TVector3 BeamDirection = TVector3(cos(BeamPhi)*sin(BeamTheta), sin(BeamPhi)*sin(BeamTheta), cos(BeamTheta));
+// cout << BeamDirection.X() << " " << BeamDirection.Y() << " " << BeamDirection.Z() << endl;
+
+ // Hit direction taking into account beam position on target
+ TVector3 HitDirection = A - TVector3(XTarget, YTarget, 0);
+// cout << "A: " << A.X() << " " << A.Y() << " " << A.Z() << endl;
+// cout << "HitDirection: " << HitDirection.X() << " " << HitDirection.Y() << " " << HitDirection.Z() << endl;
+
+ // Calculate scattering angle w.r.t. optical beam axis (do not take into account beam position on target)
+ ThetaStrip = ThetaCalculation(A, TVector3(0,0,1));
+ Theta = ThetaCalculation(Detec, TVector3(0, 0, 1));
+ // Calculate scattering angle w.r.t. beam (ideal case)
+// ThetaStrip = ThetaCalculation(HitDirection, BeamDirection);
+// Theta = ThetaCalculation(Detec - TVector3(XTarget, YTarget, 0), BeamDirection);
+ // Calculate scattering angle w.r.t. beam (finite spatial resolution)
+// double resol = 800; // in micrometer
+// angle = gene->Rndm() * 2*3.14;
+// r = fabs(gene->Gaus(0, resol)) * micrometer;
+// ThetaStrip = ThetaCalculation(A - TVector3(XTarget + r*cos(angle), YTarget + r*sin(angle), 0), BeamDirection);
+// Theta = ThetaCalculation(Detec - TVector3(XTarget + r*cos(angle), YTarget + r*sin(angle), 0), BeamDirection);
+//
+ // Correct for energy loss in the target
+ E = LightTarget.EvaluateInitialEnergy(E, myDetector->GetTargetThickness()/2 * micrometer, ThetaStrip);
+
+ // Calculate excitation energy
+// if (Theta/deg > 150 && Theta/deg < 180) {
+// if (Theta/deg < 60 && ThetaCM/deg < 90) {
+// if (Theta/deg > 35 && Theta/deg < 45 && E/MeV < 17) {
+// if (Theta/deg < 45) {
+// if (E/MeV < 38) { // for (p,t) reaction
+ if (Theta/deg > 90) { // for (d,p) reaction
+ ExNoStrips = myReaction->ReconstructRelativistic(E, Theta / rad);
+ Ex = myReaction->ReconstructRelativistic(E, ThetaStrip);
+ }
+ else {
+ Ex = -200;
+ ExNoStrips = -200;
+ }
+ }
+ else {
+ Ex = -100;
+ ExNoStrips = -100;
+ }
+
+ EE = E ; TT = ThetaStrip/deg;
+ if (E>-1000) {
+ X = A . X();
+ Y = A . Y();
+ }
+ else {
+ X = -1000 ; Y = -1000;
+ }
+
+ // Fill output tree
+ RootOutput::getInstance()->GetTree()->Fill();
+ }
+
+ // delete singleton classes
+ RootOutput::getInstance()->Destroy();
+ RootInput::getInstance()->Destroy();
+
+ return 0;
+}
+
+
+double ThetaCalculation (TVector3 A , TVector3 B)
+{
+ double Theta = acos( (A.Dot(B)) / (A.Mag()*B.Mag()) );
+ return Theta ;
+}
+
diff --git a/NPLib/GASPARD/GaspardTrackerAnnular.cxx b/NPLib/GASPARD/GaspardTrackerAnnular.cxx
new file mode 100644
index 000000000..e5a665fd9
--- /dev/null
+++ b/NPLib/GASPARD/GaspardTrackerAnnular.cxx
@@ -0,0 +1,303 @@
+#include "GaspardTrackerAnnular.h"
+
+// C++ headers
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <cmath>
+
+// Gaspard
+#include "TGaspardTrackerPhysicsNew.h"
+
+using namespace std;
+
+
+GaspardTrackerAnnular::GaspardTrackerAnnular(map<int, GaspardTrackerModule*> &Module,
+ TGaspardTrackerPhysicsNew* &EventPhysics)
+ : m_ModuleTest(Module),
+ m_EventPhysics(EventPhysics),
+ m_EventData(0),
+ m_PreTreatData(new TGaspardTrackerData),
+ m_NumberOfModule(0)
+{
+}
+
+
+
+GaspardTrackerAnnular::~GaspardTrackerAnnular()
+{
+ delete m_PreTreatData;
+}
+
+
+
+void GaspardTrackerAnnular::ReadConfiguration(string Path)
+{
+ // open config file
+ ifstream ConfigFile;
+ ConfigFile.open(Path.c_str());
+ string LineBuffer;
+ string DataBuffer;
+
+ double Z = 0, Rmin = 0, Rmax = 0;
+
+ // initialize flags
+ bool ReadingStatus = false;
+ bool check_Z = false;
+ bool check_Rmin = false;
+ bool check_Rmax = false;
+
+ // read config file
+ while (!ConfigFile.eof()) {
+ getline(ConfigFile, LineBuffer);
+
+ // If line is a GaspardXXX bloc, reading toggle to true
+ // and toggle to true flags indicating which shape is treated.
+ if (LineBuffer.compare(0, 10, "GPDAnnular") == 0) {
+ cout << "///////////////////////" << endl;
+ cout << "Annular module found:" << endl;
+ ReadingStatus = true;
+ }
+
+ // Reading Block
+ while (ReadingStatus) {
+ ConfigFile >> DataBuffer ;
+ // Comment Line
+ if (DataBuffer.compare(0, 1, "%") == 0) {
+ ConfigFile.ignore(std::numeric_limits<std::streamsize>::max(), '\n' );
+ }
+ // Finding another telescope (safety), toggle out
+ else if (DataBuffer.compare(0, 10, "GPDAnnular") == 0) {
+ cout << "WARNING: Another Module is find before standard sequence of Token, Error may occured in Telecope definition" << endl;
+ ReadingStatus = false;
+ }
+
+ //Position method
+ else if (DataBuffer.compare(0, 2, "Z=") == 0) {
+ check_Z = true;
+ ConfigFile >> DataBuffer ;
+ Z = atof(DataBuffer.c_str()) ;
+ cout << "Z: " << Z << endl;
+ }
+ else if (DataBuffer.compare(0, 5, "RMIN=") == 0) {
+ check_Rmin = true;
+ ConfigFile >> DataBuffer ;
+ Rmin = atof(DataBuffer.c_str()) ;
+ cout << "Rmin: " << Rmin << endl;
+ }
+ else if (DataBuffer.compare(0, 5, "RMAX=") == 0) {
+ check_Rmax = true;
+ ConfigFile >> DataBuffer ;
+ Rmax = atof(DataBuffer.c_str()) ;
+ cout << "Rmax: " << Rmax << endl;
+ }
+
+ /////////////////////////////////////////////////
+ // If All necessary information there, toggle out
+ if (check_Z && check_Rmin && check_Rmax) {
+ ReadingStatus = false;
+
+ // Add imodule
+ AddModule(Z, Rmin, Rmax);
+ m_ModuleTest[m_index["Annular"] + m_NumberOfModule] = this;
+
+ // reset boolean flag for point positioning
+ check_Z = false;
+ check_Rmin = false;
+ check_Rmax = false;
+ } // end test for adding a module
+ } // end while for reading block
+ } // end while for reading file
+
+ cout << endl << "/////////////////////////////" << endl<<endl;
+}
+
+
+void GaspardTrackerAnnular::PreTreat()
+{
+}
+
+
+
+void GaspardTrackerAnnular::BuildPhysicalEvent()
+{
+ // Check flags
+ bool Check_FirstStage = false;
+ bool Check_SecondStage = false;
+ bool Check_ThirdStage = false;
+
+ // Thresholds
+/*
+ double FirstStage_Front_E_Threshold = 0; double FirstStage_Front_T_Threshold = 0;
+ double FirstStage_Back_E_Threshold = 0; double FirstStage_Back_T_Threshold = 0;
+ double SecondStage_E_Threshold = 0; double SecondStage_T_Threshold = 0;
+ double ThirdStage_E_Threshold = 0; double ThirdStage_T_Threshold = 0;
+*/
+ // calculate multipicity in the first stage
+ int multXE = m_EventData->GetGPDTrkFirstStageFrontEMult();
+ int multYE = m_EventData->GetGPDTrkFirstStageBackEMult();
+ int multXT = m_EventData->GetGPDTrkFirstStageFrontTMult();
+ int multYT = m_EventData->GetGPDTrkFirstStageBackTMult();
+ // calculate multiplicity of 2nd and third stages
+ int mult2E = m_EventData->GetGPDTrkSecondStageEMult();
+ int mult2T = m_EventData->GetGPDTrkSecondStageTMult();
+ int mult3E = m_EventData->GetGPDTrkThirdStageEMult();
+ int mult3T = m_EventData->GetGPDTrkThirdStageTMult();
+
+ // Deal with multiplicity 1 for the first layer
+ if (multXE==1 && multYE==1 && multXT==1 && multYT==1) {
+ // calculate detector number
+ int det_ref = m_EventData->GetGPDTrkFirstStageFrontEDetectorNbr(0);
+ int detecXE = m_EventData->GetGPDTrkFirstStageFrontEDetectorNbr(0) / det_ref;
+ int detecXT = m_EventData->GetGPDTrkFirstStageFrontTDetectorNbr(0) / det_ref;
+ int detecYE = m_EventData->GetGPDTrkFirstStageBackEDetectorNbr(0) / det_ref;
+ int detecYT = m_EventData->GetGPDTrkFirstStageBackTDetectorNbr(0) / det_ref;
+
+ // module number starting from 0
+ det_ref -= m_index["Annular"];
+
+ // case of same detector
+ if (detecXE*detecXT*detecYE*detecYT == 1) {
+ // store module number
+ m_EventPhysics->SetModuleNumber(det_ref + m_index["Annular"]);
+ // calculate strip number
+ int stripXE = m_EventData->GetGPDTrkFirstStageFrontEStripNbr(0);
+ int stripXT = m_EventData->GetGPDTrkFirstStageFrontTStripNbr(0);
+ int stripYE = m_EventData->GetGPDTrkFirstStageBackEStripNbr(0);
+ int stripYT = m_EventData->GetGPDTrkFirstStageBackTStripNbr(0);
+
+ // case of same strips on X and Y
+ if (stripXE == stripXT && stripYE == stripYT) { // here we have a good strip event
+ // various
+ Check_FirstStage = true;
+ // store strip ID
+ m_EventPhysics->SetFirstStageFrontPosition(stripXE);
+ m_EventPhysics->SetFirstStageBackPosition(stripYE);
+ // get energy from strips and store it
+ double EnergyStripFront = m_EventData->GetGPDTrkFirstStageFrontEEnergy(0);
+ m_EventPhysics->SetFirstStageEnergy(EnergyStripFront);
+ double EnergyTot = EnergyStripFront;
+ // get time from strips and store it
+ double TimeStripBack = m_EventData->GetGPDTrkFirstStageBackTTime(0);
+ m_EventPhysics->SetFirstStageTime(TimeStripBack);
+
+ // check if we have a 2nd stage event
+ if (mult2E==1 && mult2T==1) {
+ Check_SecondStage = true;
+ double EnergySecond = m_EventData->GetGPDTrkSecondStageEEnergy(0);
+ m_EventPhysics->SetSecondStageEnergy(EnergySecond);
+ EnergyTot += EnergySecond;
+ }
+ else if (mult2E>1 || mult2T>1) {
+ cout << "Warning: multiplicity in second stage greater than in firststage" << endl;
+ }
+
+ // check if we have a third stage event
+ if (mult3E==1 && mult3T==1) {
+ Check_ThirdStage = true;
+ double EnergyThird = m_EventData->GetGPDTrkThirdStageEEnergy(0);
+ m_EventPhysics->SetThirdStageEnergy(EnergyThird);
+ EnergyTot += EnergyThird;
+ }
+ else if (mult3E>1 || mult3T>1) {
+ cout << "Warning: multiplicity in third stage greater than in firststage" << endl;
+ }
+
+ // Fill total energy
+ m_EventPhysics->SetTotalEnergy(EnergyTot);
+
+ // Fill default values for second an third stages
+ if (!Check_SecondStage) {
+ m_EventPhysics->SetSecondStageEnergy(-1000);
+ m_EventPhysics->SetSecondStageTime(-1000);
+ m_EventPhysics->SetSecondStagePosition(-1000);
+ }
+ if (!Check_ThirdStage) {
+ m_EventPhysics->SetThirdStageEnergy(-1000);
+ m_EventPhysics->SetThirdStageTime(-1000);
+ m_EventPhysics->SetThirdStagePosition(-1000);
+ }
+ }
+ else {
+ cout << "Not same strips" << endl;
+ }
+ }
+ else {
+ cout << "Not same detector" << endl;
+ }
+ }
+ else {
+/* cout << "Multiplicity is not one, it is: " << endl;
+ cout << "\tmultXE: " << multXE << endl;
+ cout << "\tmultXT: " << multXT << endl;
+ cout << "\tmultYE: " << multYE << endl;
+ cout << "\tmultYT: " << multYT << endl;*/
+ }
+}
+
+
+
+void GaspardTrackerAnnular::BuildSimplePhysicalEvent()
+{
+}
+
+
+
+
+void GaspardTrackerAnnular::AddModule(double zpos, double rmin, double rmax)
+{
+ m_NumberOfModule++;
+
+ // Characteristics
+ int NbPhiStrips = 16;
+ int NbThetaStrips = 16;
+ int NbQuadrant = 4;
+
+ // Theta & phi strips pitch
+ double thetaPitch = (rmax - rmin) / NbThetaStrips;
+ double phiPitch = 2*M_PI / (NbQuadrant*NbThetaStrips);
+
+ // Buffer object to fill Position Array
+ vector<double> lineX;
+ vector<double> lineY;
+ vector<double> lineZ;
+
+ vector< vector< double > > OneModuleStripPositionX;
+ vector< vector< double > > OneModuleStripPositionY;
+ vector< vector< double > > OneModuleStripPositionZ;
+
+ // loop on theta strips
+ for (int i = 0; i < NbThetaStrips*NbQuadrant; i++) {
+ lineX.clear();
+ lineY.clear();
+ lineZ.clear();
+
+ // center of theta strip
+ double r = rmin + thetaPitch/2 + thetaPitch*(i % NbThetaStrips);
+ cout << i << " " << i%NbThetaStrips << " " << r << endl;
+
+ // loop on phi strips
+ for (int j = 0; j < NbPhiStrips*NbQuadrant; j++) {
+ // center of phi strips
+ double phi = phiPitch/2 + phiPitch*j;
+
+ // calculate x and y projections
+ double x = r * cos(phi);
+ double y = r * sin(phi);
+
+ cout << i << " " << j << " " << r << " " << phi*180/M_PI << " " << x << " " << y << endl;
+ // fill lineX,Y,Z
+ lineX.push_back(x);
+ lineY.push_back(y);
+ lineZ.push_back(zpos);
+ }
+
+ OneModuleStripPositionX.push_back(lineX);
+ OneModuleStripPositionY.push_back(lineY);
+ OneModuleStripPositionZ.push_back(lineZ);
+ }
+
+ m_StripPositionX.push_back( OneModuleStripPositionX );
+ m_StripPositionY.push_back( OneModuleStripPositionY );
+ m_StripPositionZ.push_back( OneModuleStripPositionZ );
+}
diff --git a/NPLib/GASPARD/GaspardTrackerAnnular.h b/NPLib/GASPARD/GaspardTrackerAnnular.h
new file mode 100644
index 000000000..05e0d7e48
--- /dev/null
+++ b/NPLib/GASPARD/GaspardTrackerAnnular.h
@@ -0,0 +1,77 @@
+#ifndef GaspardTrackerAnnular_h
+#define GaspardTrackerAnnular_h 1
+
+// C++ headers
+#include <iostream>
+#include <string>
+#include <map>
+#include <vector>
+
+// Gaspard headers
+#include "TGaspardTrackerData.h"
+#include "TGaspardTrackerPhysicsNew.h"
+#include "GaspardTrackerModule.h"
+
+using namespace std;
+
+
+
+class GaspardTrackerAnnular : public GaspardTrackerModule
+{
+public:
+ ////////////////////////////////////////////////////
+ /////// Default Constructor and Destructor /////////
+ ////////////////////////////////////////////////////
+ GaspardTrackerAnnular(map<int, GaspardTrackerModule*> &Module, TGaspardTrackerPhysicsNew* &EventPhysics);
+ virtual ~GaspardTrackerAnnular();
+
+public:
+ ////////////////////////////////////////////////////
+ //// Inherite from GaspardTrackerModule class /////
+ ////////////////////////////////////////////////////
+ // Read stream at Configfile to pick-up parameters of detector (Position,...)
+ void ReadConfiguration(string Path);
+
+ // The goal of this method is to extract physical parameters from raw data
+ // Method called at each event read from the Input Tree
+ void BuildPhysicalEvent();
+
+ // Same as before but with a simpler treatment
+ void BuildSimplePhysicalEvent();
+
+private:
+ map<int, GaspardTrackerModule*> &m_ModuleTest;
+ TGaspardTrackerPhysicsNew* &m_EventPhysics;
+
+public:
+ void SetGaspardDataPointer(TGaspardTrackerData* gaspardData) {m_EventData = gaspardData;};
+ void PreTreat();
+
+private:
+ // Gaspard data coming from TGaspardTrackerPhysics through the
+ // SetGaspardDataPointer method
+ TGaspardTrackerData* m_EventData;
+ TGaspardTrackerData* m_PreTreatData;
+
+public:
+ ////////////////////////////////
+ // Specific to GaspardTracker //
+ ////////////////////////////////
+ // Add a Module
+ void AddModule(double zpos, double rmin, double rmax);
+
+ // Getters to retrieve the (X,Y,Z) coordinates of a pixel defined by strips (X,Y)
+ double GetStripPositionX(int N ,int X ,int Y) { return m_StripPositionX[N-1-m_index["Annular"]][X-1][Y-1]; }
+ double GetStripPositionY(int N ,int X ,int Y) { return m_StripPositionY[N-1-m_index["Annular"]][X-1][Y-1]; }
+ double GetStripPositionZ(int N ,int X ,int Y) { return m_StripPositionZ[N-1-m_index["Annular"]][X-1][Y-1]; }
+ double GetNumberOfModule() { return m_NumberOfModule; }
+
+private:
+ // Spatial Position of Strip Calculated on basis of detector position
+ int m_NumberOfModule;
+ vector< vector < vector < double > > > m_StripPositionX;
+ vector< vector < vector < double > > > m_StripPositionY;
+ vector< vector < vector < double > > > m_StripPositionZ;
+};
+
+#endif
diff --git a/NPLib/GASPARD/GaspardTrackerDummyShape.cxx b/NPLib/GASPARD/GaspardTrackerDummyShape.cxx
index 1fd813be7..7b61d0e5b 100644
--- a/NPLib/GASPARD/GaspardTrackerDummyShape.cxx
+++ b/NPLib/GASPARD/GaspardTrackerDummyShape.cxx
@@ -273,25 +273,16 @@ void GaspardTrackerDummyShape::BuildPhysicalEvent()
if (stripXE == stripXT && stripYE == stripYT) { // here we have a good strip event
// various
Check_FirstStage = true;
-// EventMultiplicity = 1;
// store strip ID
m_EventPhysics->SetFirstStageFrontPosition(stripXE);
m_EventPhysics->SetFirstStageBackPosition(stripYE);
// get energy from strips and store it
double EnergyStripFront = m_EventData->GetGPDTrkFirstStageFrontEEnergy(0);
-// double EnergyStripBack = m_EventData->GetGPDTrkFirstStageBackEEnergy(0);
-// double EnergyStrip = 0.5 * (EnergyStripFront + EnergyStripBack);
- double EnergyStrip = EnergyStripFront;
-// if (EnergyStripBack > EnergyStrip) EnergyStrip = EnergyStripBack;
- m_EventPhysics->SetFirstStageEnergy(EnergyStrip);
- double EnergyTot = EnergyStrip;
+ m_EventPhysics->SetFirstStageEnergy(EnergyStripFront);
+ double EnergyTot = EnergyStripFront;
// get time from strips and store it
- double TimeStripFront = m_EventData->GetGPDTrkFirstStageFrontEEnergy(0);
double TimeStripBack = m_EventData->GetGPDTrkFirstStageBackEEnergy(0);
- double TimeStrip = 0.5 * (TimeStripFront + TimeStripBack);
-// double TimeStrip = TimeStripFront;
-// if (TimeStripBack > TimeStrip) TimeStrip = TimeStripBack;
- m_EventPhysics->SetFirstStageTime(TimeStrip);
+ m_EventPhysics->SetFirstStageTime(TimeStripBack);
// check if we have a 2nd stage event
if (mult2E==1 && mult2T==1) {
@@ -303,6 +294,7 @@ void GaspardTrackerDummyShape::BuildPhysicalEvent()
else if (mult2E>1 || mult2T>1) {
cout << "Warning: multiplicity in second stage greater than in firststage" << endl;
}
+
// check if we have a third stage event
if (mult3E==1 && mult3T==1) {
Check_ThirdStage = true;
@@ -316,6 +308,18 @@ void GaspardTrackerDummyShape::BuildPhysicalEvent()
// Fill total energy
m_EventPhysics->SetTotalEnergy(EnergyTot);
+
+ // Fill default values for second an third stages
+ if (!Check_SecondStage) {
+ m_EventPhysics->SetSecondStageEnergy(-1000);
+ m_EventPhysics->SetSecondStageTime(-1000);
+ m_EventPhysics->SetSecondStagePosition(-1000);
+ }
+ if (!Check_ThirdStage) {
+ m_EventPhysics->SetThirdStageEnergy(-1000);
+ m_EventPhysics->SetThirdStageTime(-1000);
+ m_EventPhysics->SetThirdStagePosition(-1000);
+ }
}
else {
cout << "Not same strips" << endl;
diff --git a/NPLib/GASPARD/TGaspardTrackerData.h b/NPLib/GASPARD/TGaspardTrackerData.h
index bc45afa51..3eb72ac57 100644
--- a/NPLib/GASPARD/TGaspardTrackerData.h
+++ b/NPLib/GASPARD/TGaspardTrackerData.h
@@ -199,7 +199,7 @@ public:
UShort_t GetGPDTrkFirstStageFrontTStripNbr(Int_t i) {
return fGPDTrk_FirstStage_FrontT_StripNbr.at(i);
}
- Double_t GetGPDTrkFirstStageFrontTTnergy(Int_t i) {
+ Double_t GetGPDTrkFirstStageFrontTTime(Int_t i) {
return fGPDTrk_FirstStage_FrontT_Time.at(i);
}
// (Back, E)
@@ -225,7 +225,7 @@ public:
UShort_t GetGPDTrkFirstStageBackTStripNbr(Int_t i) {
return fGPDTrk_FirstStage_BackT_StripNbr.at(i);
}
- Double_t GetGPDTrkFirstStageBackTTnergy(Int_t i) {
+ Double_t GetGPDTrkFirstStageBackTTime(Int_t i) {
return fGPDTrk_FirstStage_BackT_Time.at(i);
}
diff --git a/NPSimulation/src/GaspardScorers.cc b/NPSimulation/src/GaspardScorers.cc
index 22ae3d76b..c920187fc 100644
--- a/NPSimulation/src/GaspardScorers.cc
+++ b/NPSimulation/src/GaspardScorers.cc
@@ -650,6 +650,8 @@ G4bool GPDScorerFirstStageFrontStripAnnular::ProcessHits(G4Step* aStep, G4Toucha
// Hit position in the world frame
G4ThreeVector POS = aStep->GetPreStepPoint()->GetPosition();
+ G4cout << "world frame hit position" << G4endl;
+ G4cout << POS.x() << " " << POS.y() << " " << POS.z() << G4endl;
// Hit position in the detector frame
POS = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(POS);
@@ -673,6 +675,11 @@ G4bool GPDScorerFirstStageFrontStripAnnular::ProcessHits(G4Step* aStep, G4Toucha
G4double ThetaStripNumber = floor(dummy / ThetaStripPitch);
ThetaStripNumber += PhiQuadrantNumber * GPDANNULAR::NbThetaStrips;
+ G4cout << "POS: " << POS << G4endl;
+ G4cout << "r, phi " << r << " " << phi << G4endl;
+ G4cout << "PhiWidth, PhiQuadrantNumber " << PhiWidth << " " << PhiQuadrantNumber << G4endl;
+ G4cout << "ThetaStripPitch, ThetaStripNumber, dummy " << ThetaStripPitch << " " << ThetaStripNumber << " " << dummy << G4endl;
+
if (ThetaStripNumber < 1e-6) {
/* G4cout << "POS: " << POS << G4endl;
G4cout << "r, phi " << r << " " << phi << G4endl;
diff --git a/NPSimulation/src/GaspardTrackerAnnular.cc b/NPSimulation/src/GaspardTrackerAnnular.cc
index a9e6efb95..8fc455a41 100644
--- a/NPSimulation/src/GaspardTrackerAnnular.cc
+++ b/NPSimulation/src/GaspardTrackerAnnular.cc
@@ -221,36 +221,6 @@ void GaspardTrackerAnnular::VolumeMaker(G4int TelescopeNumber ,
logicVacBox->SetVisAttributes(G4VisAttributes::Invisible);
- // Add a degrader plate between Si and CsI:
- /*
- G4Box* Degrader = new G4Box("Degrader" , 50*mm , 50*mm , 0.1*mm );
- G4LogicalVolume* logicDegrader = new G4LogicalVolume( Degrader , Harvar, "logicDegrader",0,0,0);
- PVPBuffer = new G4PVPlacement(0,G4ThreeVector(0,0,0),logicDegrader,"Degrader",logicVacBox,false,0) ;
- */
-
- //Place two marker to identify the u and v axis on silicon face:
- //marker are placed a bit before the silicon itself so they don't perturbate simulation
- //Uncomment to help debugging or if you want to understand the way the code work.
- //I should recommand to Comment it during simulation to avoid perturbation of simulation
- //Remember G4 is limitationg step on geometry constraints.
- /*
- G4ThreeVector positionMarkerU = CT*0.98 + MMu*SiliconFace/4;
- G4Box* solidMarkerU = new G4Box( "solidMarkerU" , SiliconFace/4 , 1*mm , 1*mm ) ;
- G4LogicalVolume* logicMarkerU = new G4LogicalVolume( solidMarkerU , Vacuum , "logicMarkerU",0,0,0) ;
- PVPBuffer = new G4PVPlacement(G4Transform3D(*MMrot,positionMarkerU),logicMarkerU,"MarkerU",world,false,0) ;
-
- G4VisAttributes* MarkerUVisAtt= new G4VisAttributes(G4Colour(0.,0.,0.5));//blue
- logicMarkerU->SetVisAttributes(MarkerUVisAtt);
-
- G4ThreeVector positionMarkerV = CT*0.98 + MMv*SiliconFace/4;
- G4Box* solidMarkerV = new G4Box( "solidMarkerU" , 1*mm , SiliconFace/4 , 1*mm ) ;
- G4LogicalVolume* logicMarkerV = new G4LogicalVolume( solidMarkerV , Vacuum , "logicMarkerV",0,0,0) ;
- PVPBuffer = new G4PVPlacement(G4Transform3D(*MMrot,positionMarkerV),logicMarkerV,"MarkerV",world,false,0) ;
-
- G4VisAttributes* MarkerVVisAtt= new G4VisAttributes(G4Colour(0.,0.5,0.5));//green
- logicMarkerV->SetVisAttributes(MarkerVVisAtt);
- */
-
////////////////////////////////////////////////////////////////
/////////////////// First Stage Construction////////////////////
////////////////////////////////////////////////////////////////
@@ -276,6 +246,8 @@ void GaspardTrackerAnnular::VolumeMaker(G4int TelescopeNumber ,
// Silicon detector itself
G4ThreeVector positionSilicon = G4ThreeVector(0, 0, Silicon_PosZ);
+ G4cout << "position en z de l'annulaire " << MMpos.z() * mm << G4endl;
+ G4cout << "position en z de l'annulaire " << Silicon_PosZ * mm << G4endl;
G4Tubs* solidSilicon = new G4Tubs("solidSilicon",
FirstStageRmin,
@@ -291,8 +263,7 @@ void GaspardTrackerAnnular::VolumeMaker(G4int TelescopeNumber ,
logicSilicon->SetSensitiveDetector(m_FirstStageScorer);
///Visualisation of Silicon Strip
-// G4VisAttributes* SiliconVisAtt = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5)) ;
- G4VisAttributes* SiliconVisAtt = new G4VisAttributes(G4Colour(0.0, 0.0, 0.9)) ; // bleu
+ G4VisAttributes* SiliconVisAtt = new G4VisAttributes(G4Colour(0.0, 0.0, 0.9)) ; // blue
logicSilicon->SetVisAttributes(SiliconVisAtt) ;
}
@@ -323,7 +294,6 @@ void GaspardTrackerAnnular::VolumeMaker(G4int TelescopeNumber ,
///Visualisation of Third Stage
G4VisAttributes* ThirdStageVisAtt = new G4VisAttributes(G4Colour(0.0, 0.9, 0.)) ; // green
logicThirdStage->SetVisAttributes(ThirdStageVisAtt) ;
-// logicThirdStage->SetVisAttributes(G4VisAttributes::Invisible);
// Set Third Stage sensible
logicThirdStage->SetSensitiveDetector(m_ThirdStageScorer);
@@ -464,6 +434,8 @@ void GaspardTrackerAnnular::ConstructDetector(G4LogicalVolume* world)
for (G4int i = 0; i < NumberOfModule; i++) {
// translation to position the module
+ // test if module is in the forward or backward hemisphere
+ (m_PosZ[i] < 0) ? m_PosZ[i] -= 0.5*Length : m_PosZ[i] += 0.5*Length;
MMpos = G4ThreeVector(0, 0, m_PosZ[i]);
// Passage Matrix from Lab Referential to Module Referential
--
GitLab