diff --git a/NPAnalysis/TAMU/Analysis.cxx b/NPAnalysis/TAMU/Analysis.cxx
index 12ee9b6fbf18f349f4fbd1855f5d80f2a0ff705f..c349e028825c042a874cdc407e2b3dc445b54578 100644
--- a/NPAnalysis/TAMU/Analysis.cxx
+++ b/NPAnalysis/TAMU/Analysis.cxx
@@ -1,78 +1,545 @@
#include "Analysis.h"
-using namespace std;
-
-int main(int argc, char** argv)
-{
- // command line parsing
- NPOptionManager* myOptionManager = NPOptionManager::getInstance(argc,argv);
-
- // Instantiate RootInput
- string runToReadfileName = myOptionManager->GetRunToReadFile();
- RootInput:: getInstance(runToReadfileName);
-
- // if input files are not given, use those from TAsciiFile
- if (myOptionManager->IsDefault("DetectorConfiguration")) {
- string name = RootInput::getInstance()->DumpAsciiFile("DetectorConfiguration");
- myOptionManager->SetDetectorFile(name);
- }
-
- // get input files from NPOptionManager
- string detectorfileName = myOptionManager->GetDetectorFile();
- string calibrationfileName = myOptionManager->GetCalibrationFile();
- string OutputfileName = myOptionManager->GetOutputFile();
-
- // Instantiate RootOutput
- RootOutput::getInstance("Analysis/"+OutputfileName, "AnalysedTree");
-
- // Instantiate the detector using a file
- NPA::DetectorManager* myDetector = new DetectorManager();
- myDetector->ReadConfigurationFile(detectorfileName);
- double EGammaDC[4];
- RootOutput::getInstance()->GetTree()->Branch("EGammaDC", &EGammaDC, "EGammaDC[4]/D");
- float beta=0.17;
-
- // Get the formed Chained Tree and Treat it
- TChain* Chain = RootInput:: getInstance() -> GetChain();
-
- // Get number of events to treat
- cout << endl << "///////// Starting Analysis ///////// "<< endl;
- int nentries = Chain->GetEntries();
- cout << " Number of Event to be treated : " << nentries << endl;
- clock_t begin = clock();
- clock_t end = begin;
-
- // main loop on entries
- for (int i = 0; i < nentries; i++) {
- if (i%10000 == 0 && i!=0) {
- cout.precision(5);
- end = clock();
- double TimeElapsed = (end-begin) / CLOCKS_PER_SEC;
- double percent = (double)i/nentries;
- double TimeToWait = (TimeElapsed/percent) - TimeElapsed;
- cout << " "<< flush;
- cout << "\r Progression:" << percent*100 << " % \t | \t Remaining time : ~" << TimeToWait <<"s"<< flush;
- }
- else if (i == nentries-1) cout << "\r Progression:" << " 100% " <<endl;
-
- // get data
- Chain -> GetEntry(i);
-
- myDetector->ClearEventPhysics();
- myDetector->BuildPhysicalEvent();
-
- /************************************************
-
- Put your code here
-
- ************************************************/
- RootOutput::getInstance()->GetTree()->Fill();
- }
-
- cout << "A total of " << nentries << " event has been annalysed " << endl ;
-
- RootOutput::getInstance()->Destroy();
- RootInput::getInstance()->Destroy();
- NPOptionManager::getInstance()->Destroy();
-
- return 0 ;
+
+int main(int argc, char** argv){
+ // command line parsing
+ NPOptionManager* myOptionManager = NPOptionManager::getInstance(argc,argv);
+
+ // Instantiate RootInput
+ string runToReadfileName = myOptionManager->GetRunToReadFile();
+ RootInput:: getInstance("RunToTreat.txt");
+ // Get the formed Chained Tree and Treat it
+ TChain* Chain = RootInput:: getInstance()->GetChain();
+ // if input files are not given, use those from TAsciiFile
+ if (myOptionManager->IsDefault("DetectorConfiguration")) {
+ string name = RootInput::getInstance()->DumpAsciiFile("DetectorConfiguration");
+ myOptionManager->SetDetectorFile(name);
+ }
+
+ // get input files from NPOptionManager
+ string detectorfileName = myOptionManager->GetDetectorFile();
+ string reactionfileName = myOptionManager->GetReactionFile();
+ // string calibrationfileName = myOptionManager->GetCalibrationFile();
+ string OutputfileName = myOptionManager->GetOutputFile();
+
+ // Instantiate RootOutput
+ RootOutput::getInstance("Analysis/"+OutputfileName, "AnalysedTree");
+ RootOutput::getInstance()->GetFile()->SetCompressionLevel(0);
+ // Instantiate the detector using a file
+ NPA::DetectorManager* myDetector = new DetectorManager();
+ myDetector->ReadConfigurationFile(detectorfileName);
+ // Attach new branch
+ InitOutputBranch();
+ InitInputBranch();
+
+ // Instantiate the Reaction
+ NPL::Reaction* P30dpReaction = new Reaction ;
+ P30dpReaction -> ReadConfigurationFile(reactionfileName) ;
+
+ //////////////////////////////////////////////////////
+
+ // Get pointer to the different detector
+
+ TTiaraHyballPhysics* TH = (TTiaraHyballPhysics*) myDetector -> GetDetector("TiaraHyball");
+ TTiaraBarrelPhysics* TB = (TTiaraBarrelPhysics*) myDetector -> GetDetector("TiaraBarrel");
+
+ // intermediate variable
+ //TRandom3 Rand = TRandom3();
+ TRandom *Rand = new TRandom3();
+ double ThetaNormalTarget = 0 ;
+ double ThetaTHSurface = 0;
+ double ThetaTBSurface = 0;
+ double Si_E_TH = 0 ;
+ double Si_E_TB = 0 ;
+ double Energy = 0;
+ double TargetThickness = myDetector->GetTargetThickness()*micrometer;
+
+ double XTarget = 0;
+ double YTarget =0;
+ TVector3 BeamDirection = TVector3(0,0,1);
+
+ // Get number of events to treat
+ cout << endl << "///////// Starting Analysis ///////// "<< endl;
+ int nentries = Chain->GetEntries();
+ cout << " Number of Event to be treated : " << nentries << endl;
+ clock_t begin = clock();
+ clock_t end = begin;
+ cout.precision(5);
+ ////////////////////////////////////////////////////////////////////////////
+ // main loop on entries
+ for (int i = 0; i < nentries; i++) {
+ if (i%10000 == 0 && i!=0) {
+ cout.precision(5);
+ end = clock();
+ double TimeElapsed = (end-begin) / CLOCKS_PER_SEC;
+ double percent = (double)i/nentries;
+ double TimeToWait = (TimeElapsed/percent) - TimeElapsed;
+ cout << " "<< flush;
+ cout << "\r Progression:"
+ << percent*100 << " % \t | \t Remaining time : ~"
+ << TimeToWait <<"s | Analysis Rate :"
+ << (double) i/TimeElapsed << flush;
+ }
+ else if (i == nentries-1) cout << "\r Progression:" << " 100% " <<endl;
+
+ // Get the raw data
+ Chain -> GetEntry(i);
+ // Clear previous event
+ myDetector->ClearEventPhysics();
+ // Build the current event
+ myDetector->BuildPhysicalEvent();
+ // Reinitiate calculated variable
+ ReInitValue();
+
+
+
+ // Beam energy is measured using F3 and F2 plastic TOF (time of flight)
+ double BeamEnergy = Rand->Gaus(Init->GetIncidentInitialKineticEnergy(),4.5);
+ BeamEnergy = P30_CD2.Slow(BeamEnergy,TargetThickness/2.,0);
+ P30dpReaction->SetBeamEnergy(BeamEnergy);
+
+ //////////////////////////// LOOP on TiaraHyball + SSSD Hit //////////////////
+
+ for(unsigned int countTiaraHyball = 0 ; countTiaraHyball < TH->Strip_E.size() ; countTiaraHyball++){
+ /************************************************/
+ //Part 0 : Get the useful Data
+ // TiaraHyball
+
+ /************************************************/
+ // Part 1 : Impact Angle
+ ThetaTHSurface = 0;
+ ThetaNormalTarget = 0;
+ if(XTarget>-1000 && YTarget>-1000){
+ TVector3 BeamImpact(XTarget,YTarget,0);
+ TVector3 HitDirection = TH -> GetPositionOfInteraction(countTiaraHyball) - BeamImpact ;
+ ThetaLab = HitDirection.Angle( BeamDirection );
+
+ ThetaTHSurface = HitDirection.Angle(TVector3(0,0,-1) );
+ ThetaNormalTarget = HitDirection.Angle( TVector3(0,0,1) ) ;
+ }
+
+ else{
+ BeamDirection = TVector3(-1000,-1000,-1000);
+ ThetaTHSurface = -1000 ;
+ ThetaNormalTarget = -1000 ;
+ }
+
+ /************************************************/
+
+ /************************************************/
+
+ // Part 2 : Impact Energy
+ Energy = ELab = 0;
+ Si_E_TH = TH->Strip_E[countTiaraHyball];
+ Energy = Si_E_TH;
+
+ // Evaluate energy using the thickness
+ ELab = proton_Al.EvaluateInitialEnergy( Energy ,0.4*micrometer , ThetaTHSurface);
+ // Target Correction
+ ELab = proton_CD2.EvaluateInitialEnergy( ELab ,TargetThickness/2., ThetaNormalTarget);
+ /************************************************/
+
+ /************************************************/
+ // Part 3 : Excitation Energy Calculation
+ Ex = P30dpReaction -> ReconstructRelativistic( ELab , ThetaLab );
+ int RingNumber = TH -> Strip_Ring[countTiaraHyball];
+ if(RingNumber==1){
+ Ring1Ex = Ex;
+ }
+ if(RingNumber==2){
+ Ring2Ex = Ex;
+ }
+ if(RingNumber==3){
+ Ring3Ex = Ex;
+ }
+ if(RingNumber==4){
+ Ring4Ex = Ex;
+ }
+ if(RingNumber==5){
+ Ring5Ex = Ex;
+ }
+ if(RingNumber==6){
+ Ring6Ex = Ex;
+ }
+ if(RingNumber==7){
+ Ring7Ex = Ex;
+ }
+ if(RingNumber==8){
+ Ring8Ex = Ex;
+ }
+ if(RingNumber==9){
+ Ring9Ex = Ex;
+ }
+ if(RingNumber==10){
+ Ring10Ex = Ex;
+ }
+ if(RingNumber==11){
+ Ring11Ex = Ex;
+ }
+ if(RingNumber==12){
+ Ring12Ex = Ex;
+ }
+ if(RingNumber==13){
+ Ring13Ex = Ex;
+ }
+ if(RingNumber==14){
+ Ring14Ex = Ex;
+ }
+ if(RingNumber==15){
+ Ring15Ex = Ex;
+ }
+ if(RingNumber==16){
+ Ring16Ex = Ex;
+ }
+ /************************************************/
+
+ /************************************************/
+ // Part 4 : Theta CM Calculation
+ ThetaCM = P30dpReaction -> EnergyLabToThetaCM( ELab , ThetaLab)/deg;
+ ThetaLab=ThetaLab/deg;
+
+ /************************************************/
+
+ /************************************************/
+ // Part 5a : Implementing impact matrix for the entire Hyball (all 6 wedges)
+ TVector3 ImpactPosition = TH -> GetPositionOfInteraction(countTiaraHyball);
+ int DetectorNumber = TH -> DetectorNumber[countTiaraHyball];
+ ImpactMatrixCoordX = ImpactPosition.X();
+ ImpactMatrixCoordY = ImpactPosition.Y();
+
+ // Part 5b : Implementing impact matrix for each wedge in the Hyball individually
+ if(DetectorNumber==1){
+ TVector3 DetectorOnePOI = TH -> GetPositionOfInteraction(countTiaraHyball);
+ Detector1CoordX = DetectorOnePOI.X();
+ Detector1CoordY = DetectorOnePOI.Y();
+ }
+ if(DetectorNumber==2){
+ TVector3 DetectorTwoPOI = TH -> GetPositionOfInteraction(countTiaraHyball);
+ Detector2CoordX = DetectorTwoPOI.X();
+ Detector2CoordY = DetectorTwoPOI.Y();
+ }
+ if(DetectorNumber==3){
+ TVector3 DetectorThreePOI = TH -> GetPositionOfInteraction(countTiaraHyball);
+ Detector3CoordX = DetectorThreePOI.X();
+ Detector3CoordY = DetectorThreePOI.Y();
+ }
+ if(DetectorNumber==4){
+ TVector3 DetectorFourPOI = TH -> GetPositionOfInteraction(countTiaraHyball);
+ Detector4CoordX = DetectorFourPOI.X();
+ Detector4CoordY = DetectorFourPOI.Y();
+ }
+ if(DetectorNumber==5){
+ TVector3 DetectorFivePOI = TH -> GetPositionOfInteraction(countTiaraHyball);
+ Detector5CoordX = DetectorFivePOI.X();
+ Detector5CoordY = DetectorFivePOI.Y();
+ }
+ if(DetectorNumber==6){
+ TVector3 DetectorSixPOI = TH -> GetPositionOfInteraction(countTiaraHyball);
+ Detector6CoordX = DetectorSixPOI.X();
+ Detector6CoordY = DetectorSixPOI.Y();
+ }
+ /************************************************/
+
+ /************************************************/
+ // Part 6a : Implementing randomised position impact matrix for the entire Hyball (all 6 wedges)
+ TVector3 ImpactPositionR = TH -> GetPositionOfInteraction(countTiaraHyball); // gets the event's x, y and z coordinates and puts them into a vector entitled ImpactPositionR (for randomised)
+ ImpactPositionR.SetZ(0.0);
+ DetectorNumber = TH -> DetectorNumber[countTiaraHyball];
+ double R = ImpactPositionR.Mag();
+ double Theta = ImpactPositionR.Theta(); // defines Theta, the inclination coordinate in a spherical coordinate system
+ double Phi = ImpactPositionR.Phi(); // defines Phi, the azimuthal coordinate in a spherical coordinate system
+ double RandomNumber1 = Rand->Rndm();
+ double DeltaR = ((RandomNumber1 * 6.4)-3.2);
+ R = R + DeltaR; // randomises R within a given detector ring
+ double RandomNumber2 = Rand->Rndm();
+ double DeltaAngle = ((RandomNumber2 * 0.118682389)-0.0593411946);
+ Phi = Phi + DeltaAngle; // randomises Phi within a given detector sector
+ ImpactMatrixCoordRandomX = R*(sin(Theta))*(cos(Phi)); // defines the final randomised X coordinate by transforming the randomised spherical coordinates
+ ImpactMatrixCoordRandomY = R*(sin(Phi))*(sin(Theta)); // defines the final randomised Y coordinate by transforming the randomised spherical coordinates
+
+ // Part 6b : Implementing randomised position impact matrix for each wedge in the Hyball individually
+ if(DetectorNumber==1){
+ Detector1RandomCoordX = ImpactMatrixCoordRandomX;
+ Detector1RandomCoordY = ImpactMatrixCoordRandomY;
+ }
+ if(DetectorNumber==2){
+ Detector2RandomCoordX = ImpactMatrixCoordRandomX;
+ Detector2RandomCoordY = ImpactMatrixCoordRandomY;
+ }
+ if(DetectorNumber==3){
+ Detector3RandomCoordX = ImpactMatrixCoordRandomX;
+ Detector3RandomCoordY = ImpactMatrixCoordRandomY;
+ }
+ if(DetectorNumber==4){
+ Detector4RandomCoordX = ImpactMatrixCoordRandomX;
+ Detector4RandomCoordY = ImpactMatrixCoordRandomY;
+ }
+ if(DetectorNumber==5){
+ Detector5RandomCoordX = ImpactMatrixCoordRandomX;
+ Detector5RandomCoordY = ImpactMatrixCoordRandomY;
+ }
+ if(DetectorNumber==6){
+ Detector6RandomCoordX = ImpactMatrixCoordRandomX;
+ Detector6RandomCoordY = ImpactMatrixCoordRandomY;
+ }
+
+ /************************************************/
+ } // end loop TiaraHyball*/
+
+ //////////////////////////// LOOP on TiaraBarrel //////////////////
+
+ for(unsigned int countTiaraBarrel = 0 ; countTiaraBarrel < TB->Strip_E.size() ; countTiaraBarrel++){
+ /************************************************/
+ //Part 0 : Get the useful Data
+ //TiaraBarrel
+
+ /************************************************/
+ // Part 1 : Impact Angle
+ ThetaTBSurface = 0;
+ ThetaNormalTarget = 0;
+ if(XTarget>-1000 && YTarget>-1000){
+ TVector3 BeamImpact(XTarget,YTarget,0);
+ TVector3 HitDirection = TB -> GetPositionOfInteraction(countTiaraBarrel) - BeamImpact ;
+ ThetaLab = HitDirection.Angle( BeamDirection );
+
+ ThetaTBSurface = HitDirection.Angle(TVector3(0,0,-1) );
+ ThetaNormalTarget = HitDirection.Angle( TVector3(0,0,1) ) ;
+ }
+
+ else{
+ BeamDirection = TVector3(-1000,-1000,-1000);
+ ThetaTBSurface = -1000 ;
+ ThetaNormalTarget = -1000 ;
+ }
+
+ /************************************************/
+
+ /************************************************/
+
+ // Part 2 : Impact Energy
+ Energy = ELab = 0;
+ Si_E_TB = TB->Strip_E[countTiaraBarrel];
+ Energy = Si_E_TB;
+
+ // Evaluate energy using the thickness
+ ELab = proton_Al.EvaluateInitialEnergy( Energy ,0.4*micrometer , ThetaTBSurface);
+ // Target Correction
+ ELab = proton_CD2.EvaluateInitialEnergy( ELab ,TargetThickness/2., ThetaNormalTarget); // -> have I written these lines correctly? (they look different to the syntax used in Example1) Also, should this be the beam particle or the heavy projectile particle?
+ /************************************************/
+
+ /************************************************/
+ // Part 3 : Excitation Energy Calculation
+ Ex = P30dpReaction -> ReconstructRelativistic( ELab , ThetaLab );
+ /************************************************/
+
+ /************************************************/
+ // Part 4 : Theta CM Calculation
+ ThetaCM = P30dpReaction -> EnergyLabToThetaCM( ELab , ThetaLab)/deg;
+ ThetaLab=ThetaLab/deg;
+
+ /************************************************/
+
+ /************************************************/
+ // Part 5a : Implementing impact matrix for the Tiara Barrel (all 8 detecting strips)
+ TVector3 ImpactPositionB = TB -> GetPositionOfInteraction(countTiaraBarrel);
+// int xcheck = ImpactPositionB.X();
+// cout << "event x value is " << xcheck << endl;
+// int StripNumberB = TB -> Strip_N[countTiaraBarrel];
+// cout << "Barrel resistive strip number is " << StripNumberB << endl;
+ int DetectorNumberB = TB -> DetectorNumber[countTiaraBarrel];
+// cout << "Barrel detector number is " << DetectorNumberB << endl;
+// cout << "" << endl;
+ ImpactMatrixBCoordX = ImpactPositionB.X();
+ ImpactMatrixBCoordZ = ImpactPositionB.Z();
+
+ // Part 5b : Implementing impact matrix for each wedge in the Hyball individually
+ if(DetectorNumberB==1){
+// TVector3 BarrelStripOnePOI = TH -> GetPositionOfInteraction(countTiaraBarrel);
+ BarrelStrip1CoordX = ImpactPositionB.X();
+ BarrelStrip1CoordZ = ImpactPositionB.Z();
+ }
+ if(DetectorNumberB==2){
+// TVector3 BarrelStripTwoPOI = TH -> GetPositionOfInteraction(countTiaraBarrel);
+ BarrelStrip2CoordX = ImpactPositionB.X();
+ BarrelStrip2CoordZ = ImpactPositionB.Z();
+ }
+ if(DetectorNumberB==3){
+// TVector3 BarrelStripThreePOI = TH -> GetPositionOfInteraction(countTiaraBarrel);
+ BarrelStrip3CoordX = ImpactPositionB.X();
+ BarrelStrip3CoordZ = ImpactPositionB.Z();
+ }
+ if(DetectorNumberB==4){
+// TVector3 BarrelStripFourPOI = TH -> GetPositionOfInteraction(countTiaraBarrel);
+ BarrelStrip4CoordX = ImpactPositionB.X();
+ BarrelStrip4CoordZ = ImpactPositionB.Z();
+ }
+ if(DetectorNumberB==5){
+// TVector3 BarrelStripFivePOI = TH -> GetPositionOfInteraction(countTiaraBarrel);
+ BarrelStrip5CoordX = ImpactPositionB.X();
+ BarrelStrip5CoordZ = ImpactPositionB.Z();
+ }
+ if(DetectorNumberB==6){
+// TVector3 BarrelStripSixPOI = TH -> GetPositionOfInteraction(countTiaraBarrel);
+ BarrelStrip6CoordX = ImpactPositionB.X();
+ BarrelStrip6CoordZ = ImpactPositionB.Z();
+ }
+ if(DetectorNumberB==7){
+// TVector3 BarrelStripSevenPOI = TH -> GetPositionOfInteraction(countTiaraBarrel);
+ BarrelStrip7CoordX = ImpactPositionB.X();
+ BarrelStrip7CoordZ = ImpactPositionB.Z();
+ }
+ if(DetectorNumberB==8){
+// TVector3 BarrelStripEightPOI = TH -> GetPositionOfInteraction(countTiaraBarrel);
+ BarrelStrip8CoordX = ImpactPositionB.X();
+ BarrelStrip8CoordZ = ImpactPositionB.Z();
+ }
+
+ /************************************************/
+
+ /************************************************/
+ // Part 6a : Implementing randomised position impact matrix for the entire Hyball (all 6 wedges)
+ /*TVector3 ImpactPositionR = TH -> GetPositionOfInteraction(countTiaraHyball); // gets the event's x, y and z coordinates and puts them into a vector entitled ImpactPositionR (for randomised)
+ ImpactPositionR.SetZ(0.0);
+ DetectorNumber = TH -> DetectorNumber[countTiaraHyball];
+ double R = ImpactPositionR.Mag();
+ double Theta = ImpactPositionR.Theta(); // defines Theta, the inclination coordinate in a spherical coordinate system
+ double Phi = ImpactPositionR.Phi(); // defines Phi, the azimuthal coordinate in a spherical coordinate system
+ double RandomNumber1 = Rand->Rndm();
+ double DeltaR = ((RandomNumber1 * 6.4)-3.2);
+ R = R + DeltaR; // randomises R within a given detector ring
+ double RandomNumber2 = Rand->Rndm();
+ double DeltaAngle = ((RandomNumber2 * 0.118682389)-0.0593411946);
+ Phi = Phi + DeltaAngle; // randomises Phi within a given detector sector
+ ImpactMatrixCoordRandomX = R*(sin(Theta))*(cos(Phi)); // defines the final randomised X coordinate by transforming the randomised spherical coordinates
+ ImpactMatrixCoordRandomY = R*(sin(Phi))*(sin(Theta)); // defines the final randomised Y coordinate by transforming the randomised spherical coordinates
+
+ // Part 6b : Implementing randomised position impact matrix for each wedge in the Hyball individually
+ if(DetectorNumber==1){
+ Detector1RandomCoordX = ImpactMatrixCoordRandomX;
+ Detector1RandomCoordY = ImpactMatrixCoordRandomY;
+ }
+ if(DetectorNumber==2){
+ Detector2RandomCoordX = ImpactMatrixCoordRandomX;
+ Detector2RandomCoordY = ImpactMatrixCoordRandomY;
+ }
+ if(DetectorNumber==3){
+ Detector3RandomCoordX = ImpactMatrixCoordRandomX;
+ Detector3RandomCoordY = ImpactMatrixCoordRandomY;
+ }
+ if(DetectorNumber==4){
+ Detector4RandomCoordX = ImpactMatrixCoordRandomX;
+ Detector4RandomCoordY = ImpactMatrixCoordRandomY;
+ }
+ if(DetectorNumber==5){
+ Detector5RandomCoordX = ImpactMatrixCoordRandomX;
+ Detector5RandomCoordY = ImpactMatrixCoordRandomY;
+ }
+ if(DetectorNumber==6){
+ Detector6RandomCoordX = ImpactMatrixCoordRandomX;
+ Detector6RandomCoordY = ImpactMatrixCoordRandomY;
+ }
+ */
+ /************************************************/
+ } // end loop TiaraBarrel
+
+
+////////////////////////////////////////////////////////////////
+
+ if(ELab>0)
+ RootOutput::getInstance()->GetTree()->Fill();
+ } // loop over events
+
+ cout << "A total of " << nentries << " event has been analysed " <<endl ;
+
+ RootOutput::getInstance()->Destroy();
+ RootInput::getInstance()->Destroy();
+ NPOptionManager::getInstance()->Destroy();
+ ////////////////////////////////////////////////////////////////////////////
+ return 0 ;
}
+
+////////////////////////////////////////////////////////////////////////////
+void InitOutputBranch() {
+ RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex,"Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("Ring1Ex",&Ring1Ex,"Ring1Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("Ring2Ex",&Ring2Ex,"Ring2Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("Ring3Ex",&Ring3Ex,"Ring3Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("Ring4Ex",&Ring4Ex,"Ring4Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("Ring5Ex",&Ring5Ex,"Ring5Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("Ring6Ex",&Ring6Ex,"Ring6Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("Ring7Ex",&Ring7Ex,"Ring7Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("Ring8Ex",&Ring8Ex,"Ring8Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("Ring9Ex",&Ring9Ex,"Ring9Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("Ring10Ex",&Ring10Ex,"Ring10Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("Ring11Ex",&Ring11Ex,"Ring11Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("Ring12Ex",&Ring12Ex,"Ring12Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("Ring13Ex",&Ring13Ex,"Ring13Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("Ring14Ex",&Ring14Ex,"Ring14Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("Ring15Ex",&Ring15Ex,"Ring15Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("Ring16Ex",&Ring16Ex,"Ring16Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab,"ELab/D");
+ RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab,"ThetaLab/D");
+ RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM,"ThetaCM/D");
+ RootOutput::getInstance()->GetTree()->Branch("ImpactMatrixCoordX",&ImpactMatrixCoordX,"ImpactMatrixCoordX/D");
+ RootOutput::getInstance()->GetTree()->Branch("ImpactMatrixCoordY",&ImpactMatrixCoordY,"ImpactMatrixCoordY/D");
+ RootOutput::getInstance()->GetTree()->Branch("Detector1CoordX",&Detector1CoordX,"Detector1CoordX/D");
+ RootOutput::getInstance()->GetTree()->Branch("Detector1CoordY",&Detector1CoordY,"Detector1CoordY/D");
+ RootOutput::getInstance()->GetTree()->Branch("Detector2CoordX",&Detector2CoordX,"Detector2CoordX/D");
+ RootOutput::getInstance()->GetTree()->Branch("Detector2CoordY",&Detector2CoordY,"Detector2CoordY/D");
+ RootOutput::getInstance()->GetTree()->Branch("Detector3CoordX",&Detector3CoordX,"Detector3CoordX/D");
+ RootOutput::getInstance()->GetTree()->Branch("Detector3CoordY",&Detector3CoordY,"Detector3CoordY/D");
+ RootOutput::getInstance()->GetTree()->Branch("Detector4CoordX",&Detector4CoordX,"Detector4CoordX/D");
+ RootOutput::getInstance()->GetTree()->Branch("Detector4CoordY",&Detector4CoordY,"Detector4CoordY/D");
+ RootOutput::getInstance()->GetTree()->Branch("Detector5CoordX",&Detector5CoordX,"Detector5CoordX/D");
+ RootOutput::getInstance()->GetTree()->Branch("Detector5CoordY",&Detector5CoordY,"Detector5CoordY/D");
+ RootOutput::getInstance()->GetTree()->Branch("Detector6CoordX",&Detector6CoordX,"Detector6CoordX/D");
+ RootOutput::getInstance()->GetTree()->Branch("Detector6CoordY",&Detector6CoordY,"Detector6CoordY/D");
+ RootOutput::getInstance()->GetTree()->Branch("ImpactMatrixCoordRandomX",&ImpactMatrixCoordRandomX,"ImpactMatrixCoordRandomX/D");
+ RootOutput::getInstance()->GetTree()->Branch("ImpactMatrixCoordRandomY",&ImpactMatrixCoordRandomY,"ImpactMatrixCoordRandomY/D");
+ RootOutput::getInstance()->GetTree()->Branch("Detector1RandomCoordX",&Detector1RandomCoordX,"Detector1RandomCoordX/D");
+ RootOutput::getInstance()->GetTree()->Branch("Detector1RandomCoordY",&Detector1RandomCoordY,"Detector1RandomCoordY/D");
+ RootOutput::getInstance()->GetTree()->Branch("Detector2RandomCoordX",&Detector2RandomCoordX,"Detector2RandomCoordX/D");
+ RootOutput::getInstance()->GetTree()->Branch("Detector2RandomCoordY",&Detector2RandomCoordY,"Detector2RandomCoordY/D");
+ RootOutput::getInstance()->GetTree()->Branch("Detector3RandomCoordX",&Detector3RandomCoordX,"Detector3RandomCoordX/D");
+ RootOutput::getInstance()->GetTree()->Branch("Detector3RandomCoordY",&Detector3RandomCoordY,"Detector3RandomCoordY/D");
+ RootOutput::getInstance()->GetTree()->Branch("Detector4RandomCoordX",&Detector4RandomCoordX,"Detector4RandomCoordX/D");
+ RootOutput::getInstance()->GetTree()->Branch("Detector4RandomCoordY",&Detector4RandomCoordY,"Detector4RandomCoordY/D");
+ RootOutput::getInstance()->GetTree()->Branch("Detector5RandomCoordX",&Detector5RandomCoordX,"Detector5RandomCoordX/D");
+ RootOutput::getInstance()->GetTree()->Branch("Detector5RandomCoordY",&Detector5RandomCoordY,"Detector5RandomCoordY/D");
+ RootOutput::getInstance()->GetTree()->Branch("Detector6RandomCoordX",&Detector6RandomCoordX,"Detector6RandomCoordX/D");
+ RootOutput::getInstance()->GetTree()->Branch("Detector6RandomCoordY",&Detector6RandomCoordY,"Detector6RandomCoordY/D");
+ RootOutput::getInstance()->GetTree()->Branch("ImpactMatrixBCoordX",&ImpactMatrixBCoordX,"ImpactMatrixBCoordX/D");
+ RootOutput::getInstance()->GetTree()->Branch("ImpactMatrixBCoordZ",&ImpactMatrixBCoordZ,"ImpactMatrixBCoordZ/D");
+ RootOutput::getInstance()->GetTree()->Branch("BarrelStrip1CoordX",&BarrelStrip1CoordX,"BarrelStrip1CoordX/D");
+ RootOutput::getInstance()->GetTree()->Branch("BarrelStrip1CoordZ",&BarrelStrip1CoordZ,"BarrelStrip1CoordZ/D");
+ RootOutput::getInstance()->GetTree()->Branch("BarrelStrip2CoordX",&BarrelStrip2CoordX,"BarrelStrip2CoordX/D");
+ RootOutput::getInstance()->GetTree()->Branch("BarrelStrip2CoordZ",&BarrelStrip2CoordZ,"BarrelStrip2CoordZ/D");
+ RootOutput::getInstance()->GetTree()->Branch("BarrelStrip3CoordX",&BarrelStrip3CoordX,"BarrelStrip3CoordX/D");
+ RootOutput::getInstance()->GetTree()->Branch("BarrelStrip3CoordZ",&BarrelStrip3CoordZ,"BarrelStrip3CoordZ/D");
+ RootOutput::getInstance()->GetTree()->Branch("BarrelStrip4CoordX",&BarrelStrip4CoordX,"BarrelStrip4CoordX/D");
+ RootOutput::getInstance()->GetTree()->Branch("BarrelStrip4CoordZ",&BarrelStrip4CoordZ,"BarrelStrip4CoordZ/D");
+ RootOutput::getInstance()->GetTree()->Branch("BarrelStrip5CoordX",&BarrelStrip5CoordX,"BarrelStrip5CoordX/D");
+ RootOutput::getInstance()->GetTree()->Branch("BarrelStrip5CoordZ",&BarrelStrip5CoordZ,"BarrelStrip5CoordZ/D");
+ RootOutput::getInstance()->GetTree()->Branch("BarrelStrip6CoordX",&BarrelStrip6CoordX,"BarrelStrip6CoordX/D");
+ RootOutput::getInstance()->GetTree()->Branch("BarrelStrip6CoordZ",&BarrelStrip6CoordZ,"BarrelStrip6CoordZ/D");
+ RootOutput::getInstance()->GetTree()->Branch("BarrelStrip7CoordX",&BarrelStrip7CoordX,"BarrelStrip7CoordX/D");
+ RootOutput::getInstance()->GetTree()->Branch("BarrelStrip7CoordZ",&BarrelStrip7CoordZ,"BarrelStrip7CoordZ/D");
+ RootOutput::getInstance()->GetTree()->Branch("BarrelStrip8CoordX",&BarrelStrip8CoordX,"BarrelStrip8CoordX/D");
+ RootOutput::getInstance()->GetTree()->Branch("BarrelStrip8CoordZ",&BarrelStrip8CoordZ,"BarrelStrip8CoordZ/D");
+}
+
+/////////////////////////////////////////////////////////////////////////////
+void InitInputBranch(){
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("InitialConditions",&Init );
+ RootInput:: getInstance()->GetChain()->SetBranchStatus("InitialConditions",true );
+ RootInput:: getInstance()->GetChain()->SetBranchStatus("fIC_*",true );
+}
+/////////////////////////////////////////////////////////////////////////////
+void ReInitValue(){
+ Ex = -1000 ;
+ ELab = -1000;
+ ThetaLab = -1000;
+ ThetaCM = -1000;
+}
+
+
+/////////////////////////////////////////////////////////////////////////////
diff --git a/NPAnalysis/TAMU/Analysis.h b/NPAnalysis/TAMU/Analysis.h
index d5ca8c305a8f3b7d788651234860e30edb569327..8ab72955d7329753094dc6552968f8a59ca014aa 100644
--- a/NPAnalysis/TAMU/Analysis.h
+++ b/NPAnalysis/TAMU/Analysis.h
@@ -5,10 +5,17 @@
/////////////////////////////////////////////////////////////////////////////////////////////////
// -------------------------------------- VARIOUS INCLUDE ---------------------------------------
-// NPA
+// NPL
#include "DetectorManager.h"
#include "NPOptionManager.h"
-
+#include "NPReaction.h"
+#include "RootInput.h"
+#include "RootOutput.h"
+#include "TTiaraHyballPhysics.h"
+#include "TTiaraBarrelPhysics.h"
+#include "TInitialConditions.h"
+#include "NPEnergyLoss.h"
+using namespace NPL ;
// STL C++
#include <iostream>
#include <fstream>
@@ -16,107 +23,124 @@
#include <string>
#include <cmath>
#include <cstdlib>
-
+using namespace std;
// ROOT
#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
-#include <TLeaf.h>
#include <TVector3.h>
-#include <TRandom.h>
-
-// NPL
-#include "RootInput.h"
-#include "RootOutput.h"
-#include "NPReaction.h"
-#include "TInitialConditions.h"
-#include "TPlasticData.h"
-#include "TMust2Data.h"
-#include "TMust2Physics.h"
-#include "TExogamPhysics.h"
-#include "TSSSDPhysics.h"
-#include "TPlasticPhysics.h"
-#include "GaspardTracker.h"
-
-// Use CLHEP System of unit and Physical Constant
-#include "NPGlobalSystemOfUnits.h"
-#include "NPPhysicalConstants.h"
-
+#include <TRandom3.h>
+#include <TMath.h>
+#include <TObject.h>
// ----------------------------------------------------------------------------------------------
-double ThetaCalculation (TVector3 A , TVector3 B) ;
+void InitOutputBranch() ;
+void InitInputBranch() ;
+void ReInitValue() ;
/////////////////////////////////////////////////////////////////////////////////////////////////
// ----------------------------------- DOUBLE, INT, BOOL AND MORE -------------------------------
-namespace VARIABLE
- {
- // Declare your Variable here:
-
- double X1,Y1,Z1 ;
- int N1,N2 = 0 ;
- bool check= false ;
-
- // A Usefull Simple Random Generator
- TRandom Rand;
- }
-
+namespace VARIABLE{
+ double Ex;
+ double Ring1Ex;
+ double Ring2Ex;
+ double Ring3Ex;
+ double Ring4Ex;
+ double Ring5Ex;
+ double Ring6Ex;
+ double Ring7Ex;
+ double Ring8Ex;
+ double Ring9Ex;
+ double Ring10Ex;
+ double Ring11Ex;
+ double Ring12Ex;
+ double Ring13Ex;
+ double Ring14Ex;
+ double Ring15Ex;
+ double Ring16Ex;
+ double ELab;
+ double ThetaLab;
+ double ThetaCM;
+ double ImpactMatrixCoordX;
+ double ImpactMatrixCoordY;
+ double Detector1CoordX;
+ double Detector1CoordY;
+ double Detector2CoordX;
+ double Detector2CoordY;
+ double Detector3CoordX;
+ double Detector3CoordY;
+ double Detector4CoordX;
+ double Detector4CoordY;
+ double Detector5CoordX;
+ double Detector5CoordY;
+ double Detector6CoordX;
+ double Detector6CoordY;
+ TVector3 DetectorOnePOI;
+ TVector3 DetectorTwoPOI;
+ TVector3 DetectorThreePOI;
+ TVector3 DetectorFourPOI;
+ TVector3 DetectorFivePOI;
+ TVector3 DetectorSixPOI;
+ double BarrelStrip1CoordX;
+ double BarrelStrip1CoordZ;
+ double BarrelStrip2CoordX;
+ double BarrelStrip2CoordZ;
+ double BarrelStrip3CoordX;
+ double BarrelStrip3CoordZ;
+ double BarrelStrip4CoordX;
+ double BarrelStrip4CoordZ;
+ double BarrelStrip5CoordX;
+ double BarrelStrip5CoordZ;
+ double BarrelStrip6CoordX;
+ double BarrelStrip6CoordZ;
+ double BarrelStrip7CoordX;
+ double BarrelStrip7CoordZ;
+ double BarrelStrip8CoordX;
+ double BarrelStrip8CoordZ;
+ TVector3 BarrelStripOnePOI;
+ TVector3 BarrelStripTwoPOI;
+ TVector3 BarrelStripThreePOI;
+ TVector3 BarrelStripFourPOI;
+ TVector3 BarrelStripFivePOI;
+ TVector3 BarrelStripSixPOI;
+ TVector3 BarrelStripSevenPOI;
+ TVector3 BarrelStripEightPOI;
+ double ImpactMatrixCoordRandomX;
+ double ImpactMatrixCoordRandomY;
+ double Detector1RandomCoordX;
+ double Detector1RandomCoordY;
+ double Detector2RandomCoordX;
+ double Detector2RandomCoordY;
+ double Detector3RandomCoordX;
+ double Detector3RandomCoordY;
+ double Detector4RandomCoordX;
+ double Detector4RandomCoordY;
+ double Detector5RandomCoordX;
+ double Detector5RandomCoordY;
+ double Detector6RandomCoordX;
+ double Detector6RandomCoordY;
+ TVector3 DetectorOnePOIR;
+ TVector3 DetectorTwoPOIR;
+ TVector3 DetectorThreePOIR;
+ TVector3 DetectorFourPOIR;
+ TVector3 DetectorFivePOIR;
+ TVector3 DetectorSixPOIR;
+ double ImpactMatrixBCoordX;
+ double ImpactMatrixBCoordZ;
+ TInitialConditions* Init = new TInitialConditions();
+}
+
using namespace VARIABLE ;
// ----------------------------------------------------------------------------------------------
-
-
-
/////////////////////////////////////////////////////////////////////////////////////////////////
-// -----------------------------------GRAPH------------------------------------------------------
-#include <TObject.h>
-#include <TH1.h>
-#include <TH1F.h>
-#include <TH2.h>
-#include <TH2F.h>
-#include <TGraph2D.h>
-
-namespace GRAPH
- {
- // Declare your Spectra here:
-
- TH1F *myHist1D = new TH1F("Hist1D","Histogramm 1D ; x ; count", 1000 , -5 , 5 ) ;
-
- TH2F *myHist2D = new TH2F("Hist2D","Histogramm 2D ; x ; y ", 128 , 1 , 128 , 128 , 1 , 128 ) ;
-
- }
-
-using namespace GRAPH ;
-// --------------------------------------------------------------------------------------------
-
-
-
-///////////////////////////////////////////////////////////////////////////////////////////////
-// -----------------------------------CUT------------------------------------------------------
-#include <TCutG.h>
-namespace CUT
- {
- // Declare your Cut here:
-
- }
-
-using namespace CUT ;
-// --------------------------------------------------------------------------------------------
-
-
-
-////////////////////////////////////////////////////////////////////////////////////////////////
// -----------------------------------ENERGY LOSS----------------------------------------------
-#include "NPEnergyLoss.h"
-using namespace NPL ;
-namespace ENERGYLOSS
- {
-
- // Declare your Energy loss here :
- /* EnergyLoss ProtonTarget = EnergyLoss ( "CD2.txt" ,
- 100 ,
- 1,
- 1 );
- */
- }
-
+namespace ENERGYLOSS{
+ // Energy loss table: the G4Table are generated by the simulation
+ EnergyLoss proton_CD2 = EnergyLoss("proton_CD2.G4table","G4Table",100 );
+ EnergyLoss proton_Al = EnergyLoss("proton_Al.G4table","G4Table",10);
+ EnergyLoss proton_Si = EnergyLoss("proton_Si.G4table","G4Table",10);
+ EnergyLoss P30_CD2 = EnergyLoss("P30[0.0]_CD2.G4table","G4Table",100);
+}
+
using namespace ENERGYLOSS ;
// ----------------------------------------------------------------------------------------------
/////////////////////////////////////////////////////////////////////////////////////////////////
diff --git a/NPLib/Tiara/TTiaraHyballPhysics.cxx b/NPLib/Tiara/TTiaraHyballPhysics.cxx
index 453e8c9a1487a567ab9174ab2930ea2cb7683909..4c5bd042caf786a6a12e0918c26894a0640f2908 100644
--- a/NPLib/Tiara/TTiaraHyballPhysics.cxx
+++ b/NPLib/Tiara/TTiaraHyballPhysics.cxx
@@ -589,8 +589,8 @@ void TTiaraHyballPhysics::AddWedgeDetector( double R,double Phi,double Z){
double Wedge_R_Min = 32.6+R;
double Wedge_R_Max = 135.1+R;
- double Wedge_Phi_Min = -27.2*deg/rad ;
- double Wedge_Phi_Max = 27.2*deg/rad ;
+ double Wedge_Phi_Min = -27.2*deg/rad;
+ double Wedge_Phi_Max = 27.2*deg/rad;
Phi= Phi*deg/rad;
int Wedge_Ring_NumberOfStrip = 16 ;
@@ -643,10 +643,17 @@ TVector3 TTiaraHyballPhysics::GetDetectorNormal( const int i) const{
}
+//double TTiaraHyballPhysics::GetNumberOfDetector(const int i) const{
+// double DetectorNumber = double ( GetDetectorNumber() );
+
+// return(DetectorNumber) ;
+
+//}
+
TVector3 TTiaraHyballPhysics::GetPositionOfInteraction(const int i) const{
TVector3 Position = TVector3 ( GetStripPositionX(DetectorNumber[i],Strip_Ring[i],Strip_Sector[i] ) ,
GetStripPositionY( DetectorNumber[i],Strip_Ring[i],Strip_Sector[i] ) ,
- GetStripPositionZ( DetectorNumber[i],Strip_Ring[i],Strip_Sector[i] ) ) ;
+ GetStripPositionZ( DetectorNumber[i],Strip_Ring[i],Strip_Sector[i] ) ) ;
return(Position) ;
diff --git a/NPLib/Tiara/TTiaraHyballPhysics.h b/NPLib/Tiara/TTiaraHyballPhysics.h
index d8b0b9c2fcb529b63e8a750915f7d629a27dda59..9aa797e81f002f3db66d2bde8781fd30ffc8a0ad 100644
--- a/NPLib/Tiara/TTiaraHyballPhysics.h
+++ b/NPLib/Tiara/TTiaraHyballPhysics.h
@@ -151,7 +151,7 @@ class TTiaraHyballPhysics : public TObject, public NPA::VDetector{
double GetStripPositionY( const int N , const int Ring , const int Sector ) const{ return m_StripPositionY[N-1][Ring-1][Sector-1] ; } ;
double GetStripPositionZ( const int N , const int Ring , const int Sector ) const{ return m_StripPositionZ[N-1][Ring-1][Sector-1] ; } ;
- double GetNumberOfDetector() const { return m_NumberOfDetector; };
+// double GetNumberOfDetector() const { return m_NumberOfDetector; };
// To be called after a build Physical Event
int GetEventMultiplicity() const { return EventMultiplicity; };
diff --git a/NPSimulation/Tiara/Tiara.cc b/NPSimulation/Tiara/Tiara.cc
index 7970be9f452b25d68cbe060d4b63468d15eda8a1..8e6c4136309bd91031c4577abe185673eee8e1bf 100644
--- a/NPSimulation/Tiara/Tiara.cc
+++ b/NPSimulation/Tiara/Tiara.cc
@@ -115,13 +115,13 @@ void Tiara::ReadConfiguration(string Path){
}
// Inner Barrel case
- else if (DataBuffer=="TiaraOuterBarrel="){
+ else if (DataBuffer=="TiaraInnerBarrel="){
if(VerboseLevel==1) G4cout << "Inner Barrel found " << G4endl ;
ConfigFile >> m_boolInner;
}
// Outter Barrel case
- else if (DataBuffer=="TiaraInnerBarrel="){
+ else if (DataBuffer=="TiaraOuterBarrel="){
if(VerboseLevel==1) G4cout << "Outer Barrel found " << G4endl ;
ConfigFile >> m_boolOuter;
}
@@ -269,6 +269,9 @@ void Tiara::ReadSensitive(const G4Event* event){
// Front Energy
double EF = RandGauss::shoot(Info[0],ResoEnergy);
if(EF>EnergyThreshold){
+ int RingNumber=Info[8];
+ RingNumber=abs(RingNumber-17); //
+ Info[8]=RingNumber;
m_EventHyball->SetRingE(Info[7],Info[8],EF);
m_EventHyball->SetRingT(Info[7],Info[8],Info[1]);
}
diff --git a/NPSimulation/src/SiliconScorers.cc b/NPSimulation/src/SiliconScorers.cc
index d23e61dbda7cbef1a2fda63f8a7610767f63e906..f0f3604f1a2964890c86b280b11f4de954cf3621 100644
--- a/NPSimulation/src/SiliconScorers.cc
+++ b/NPSimulation/src/SiliconScorers.cc
@@ -278,12 +278,12 @@ G4bool PS_Silicon_Resistive::ProcessHits(G4Step* aStep, G4TouchableHistory*){
m_Position = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(m_Position);
- m_StripWidthNumber = (int)((m_Position.y() + m_StripPlaneWidth / 2.) / m_StripPitchWidth ) + 1 ;
- m_StripWidthNumber = m_NumberOfStripWidth - m_StripWidthNumber + 1 ;
+ m_StripWidthNumber = (int)((m_Position.x() + m_StripPlaneWidth / 2.) / m_StripPitchWidth ) + 1 ;
+ // m_StripWidthNumber = m_NumberOfStripWidth - m_StripWidthNumber + 1 ;
// The energy is divided in two depending on the position
// position along the resistive strip
- double P = (m_Position.x())/(0.5*m_StripPlaneLength);
+ double P = (m_Position.z())/(0.5*m_StripPlaneLength);
// Upstream Energy
EnergyAndTime[0] = aStep->GetTotalEnergyDeposit()*(1+P)*0.5;