Skip to content
Snippets Groups Projects
Commit 8cf70b25 authored by Hugo Jacob's avatar Hugo Jacob
Browse files

Added features to Exogam data and physics, still needs some work to be fully operational

parent eb001ff6
No related branches found
No related tags found
1 merge request!27Draft: [Epic] Preparation of the environement for the new GaseousDetectorScorers...
Pipeline #286770 passed
Stage: build-NPLib
Stage: build-NPSimulation
Stage: test
Hide whitespace changes
Inline Side-by-side
NPLib/Detectors/Exogam/Geometry_Clover_Exogam.h 0 → 100644
+ 577
− 0
View file @ 8cf70b25
/************************************************************************************************
The first function return the structur "Clover_struc" for the flange number given in argument (from 1 to 16 in EXOGAM)
@ The convention :
theta = scattering angle by respect to the beam axis
phi = azimutal angle
phi convention : 0 degree for flange 12 (top-theta=90 deg) and we go in the clockwise looking in the beam direction
BEWARE that the output phi is the angle of the PROJECTION on the (0,x,y) plan of
the gamma vector !!!!!!!!!!!!!!! (see ROOT manual Chapter 18)
@ Output are in rad
*
@ The distance between the clover and the target has to be defined by the user (see #define)
@ The routine works only with the ROOT Lib
@ To print the result : switch COMMENTFLAG to 1
************************************************************************************************
The second function make the Doppler correction of an input energy according to the scattering
(see function at the end of the file for input)
*************************************************************************************************
E. Clment CEA-Saclay/SPhN
June 2006
E553 update E.Clement GANIL
Sept 2008
E530 update N. de Sereville IPNO (target can now be at any position)
April 2009
E553 update E.Clement GANIL
May 2009: *more correct crystal size and can thickness
*segment mean angle depends on the interaction depth due to the crystal shape
************************************************************************************************
************************************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <ctype.h>
#include "TROOT.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TH2F.h"
#include "TH1F.h"
#include "TH3F.h"
#include "TChain.h"
#include "TTree.h"
#include "TNtuple.h"
#include "TFile.h"
#include "TRandom.h"
#include "TVector3.h"
#include "TRotation.h"
#include "TLorentzVector.h"
#include "TGenPhaseSpace.h"
#include "TMath.h"
#include "TRint.h"
#include "TString.h"
#include "TCutG.h"
#include "TLatex.h"
#include "TGraphErrors.h"
#include "TArrow.h"
#include "TF1.h"
#include "TLatex.h"
#include "TDatime.h"
#include "TSystem.h"
struct Clover_struc {
float Theta_Clover; // clover center scattering angle by respect to beam axe
float Phi_Clover; // clover center phi angle of the PROJECTION on the (0,x,y) plan of the gamma vector
float X_Clover,Y_Clover,Z_Clover;
float Theta_Crystal[4]; //crystalA=0 crystaB=1 etc ....
float Phi_Crystal[4];
float X_Crystal[4],Y_Crystal[4],Z_Crystal[4];
float Theta_Crystal_Seg[4][4]; //CrystalA Segment1 =0,0 CrystalA Segment2=0,1 etc ...
float Phi_Crystal_Seg[4][4];
float X_Crystal_Seg[4][4],Y_Crystal_Seg[4][4],Z_Crystal_Seg[4][4];
};
#define COMMENTFLAG 0
#define D_CloveFlange_Targ1 (140.0 +7.0) //mm + distance capot-crystal (np18-22-03)
#define D_CloveFlange_Targ2 (140.0 +7.0) //mm
#define D_CloveFlange_Targ3 (140.0 +7.0) //mm
#define D_CloveFlange_Targ4 (140.0 +7.0) //mm
#define D_CloveFlange_Targ5 (140.0 +7.0) //mm
#define D_CloveFlange_Targ6 (140.0 +7.0) //mm
#define D_CloveFlange_Targ7 (140.0 +7.0) //mm
#define D_CloveFlange_Targ8 (140.0 +7.0) //mm
#define D_CloveFlange_Targ9 (140.0 +7.0) //mm
#define D_CloveFlange_Targ10 (140.0 +7.0) //mm
#define D_CloveFlange_Targ11 (140.0 +7.0) //mm
#define D_CloveFlange_Targ12 (140.0 +7.0) //mm
#define D_CloveFlange_Targ13 (140.0 +7.0) //mm
#define D_CloveFlange_Targ14 (140.0 +7.0) //mm
#define D_CloveFlange_Targ15 (140.0 +7.0) //mm
#define D_CloveFlange_Targ16 (140.0 +7.0) //mm
#define D_CloveFlange_Targ17 (140.0 +7.0) //mm
#define TARGET_POSITION_X 0. // mm
#define TARGET_POSITION_Y 0. // mm
#define TARGET_POSITION_Z 0. // mm
#define InteractionDepth 20. // mm
struct Clover_struc Ask_For_Angles(int flange, double InterDepth){
struct Clover_struc Result;
int i,j;
float Real_position[17] ={D_CloveFlange_Targ1, D_CloveFlange_Targ2, D_CloveFlange_Targ3, D_CloveFlange_Targ4, D_CloveFlange_Targ5, D_CloveFlange_Targ6,
D_CloveFlange_Targ7, D_CloveFlange_Targ8, D_CloveFlange_Targ9, D_CloveFlange_Targ10, D_CloveFlange_Targ11,
D_CloveFlange_Targ12, D_CloveFlange_Targ13, D_CloveFlange_Targ14, D_CloveFlange_Targ15, D_CloveFlange_Targ16,D_CloveFlange_Targ17};
//Elements definition
TVector3 flange12(D_CloveFlange_Targ12,0,0); //clover vector
TVector3 flange12Crist[4]; //crystals vectors
TVector3 flange12CristSeg[4][4]; //segments vectors
//Initialisation by default
TVector3 v1(1,0,0);
TVector3 v2(0,1,0);
TVector3 v3(0,0,1);
for(i=0;i<4;i++){
flange12Crist[i].SetX(1.);
flange12Crist[i].SetY(1.);
flange12Crist[i].SetZ(1.);
flange12Crist[i].SetTheta(1.);
flange12Crist[i].SetPhi(1.);
flange12Crist[i].SetMag(1.);
for(j=0;j<4;j++){
flange12CristSeg[i][j].SetX(1.);
flange12CristSeg[i][j].SetY(1.);
flange12CristSeg[i][j].SetZ(1.);
flange12CristSeg[i][j].SetTheta(1.);
flange12CristSeg[i][j].SetPhi(1.);
flange12CristSeg[i][j].SetMag(1.);
}
}
/*Initialisation of all the vectors to the flange 12 BUT with the correct distance for the asked clover
between target and germanium crystal*/
//Clover position
flange12.SetTheta(90.0*TMath::Pi()/180.0);
flange12.SetPhi(0.0*TMath::Pi()/180.0);
flange12.SetMag(Real_position[flange-1]);
//segment mean angle depends on the interaction depth due to the crystal shape; this creates a linear function with InteractionDepth as input
// see plan np18-22-05 && np18-22-04
TF1 *ShapeC = new TF1("ShapeC","0.132*x+20.54");
TF1 *ShapeS1 = new TF1("ShapeS1","0.273*x+30.81"); //gd
TF1 *ShapeS2 = new TF1("ShapeS2","0.066*x+10.27"); //pt
float EXOGAM_Crystal_Center;
float EXOGAM_Segment_Pos1,EXOGAM_Segment_Pos2 ;
// Original code by E. Clément using fixed interaction depth
//if(InteractionDepth>=30){EXOGAM_Crystal_Center= 24.5;EXOGAM_Segment_Pos1=39.; EXOGAM_Segment_Pos2=12.25; }
//else{EXOGAM_Crystal_Center=ShapeC->Eval(InteractionDepth);
// EXOGAM_Segment_Pos1=ShapeS1->Eval(InteractionDepth); //gd
// EXOGAM_Segment_Pos2=ShapeS2->Eval(InteractionDepth); //pt
//}
// Modified code by H. Jacob using input interaction depth
if(InterDepth>=30){EXOGAM_Crystal_Center= 24.5;EXOGAM_Segment_Pos1=39.; EXOGAM_Segment_Pos2=12.25; }
else{EXOGAM_Crystal_Center=ShapeC->Eval(InterDepth);
EXOGAM_Segment_Pos1=ShapeS1->Eval(InterDepth); //gd
EXOGAM_Segment_Pos2=ShapeS2->Eval(InterDepth); //pt
}
//Crystal A
flange12Crist[0].SetY(flange12.Y()+(EXOGAM_Crystal_Center));
flange12Crist[0].SetZ(flange12.Z()-(EXOGAM_Crystal_Center));
flange12Crist[0].SetX(Real_position[flange-1]);
//segment1
flange12CristSeg[0][0].SetY(flange12.Y()+(EXOGAM_Segment_Pos1)); //gd
flange12CristSeg[0][0].SetZ(flange12.Z()-(EXOGAM_Segment_Pos1)); //gd
flange12CristSeg[0][0].SetX(Real_position[flange-1]);
//segment2
flange12CristSeg[0][1].SetY(flange12.Y()+(EXOGAM_Segment_Pos2));
flange12CristSeg[0][1].SetZ(flange12.Z()-(EXOGAM_Segment_Pos1));
flange12CristSeg[0][1].SetX(Real_position[flange-1]);
//segment3
flange12CristSeg[0][2].SetY(flange12.Y()+(EXOGAM_Segment_Pos2));
flange12CristSeg[0][2].SetZ(flange12.Z()-(EXOGAM_Segment_Pos2));
flange12CristSeg[0][2].SetX(Real_position[flange-1]);
//segment4
flange12CristSeg[0][3].SetY(flange12.Y()+(EXOGAM_Segment_Pos1));
flange12CristSeg[0][3].SetZ(flange12.Z()-(EXOGAM_Segment_Pos2));
flange12CristSeg[0][3].SetX(Real_position[flange-1]);
//Crystal B
flange12Crist[1].SetY(flange12.Y()-(EXOGAM_Crystal_Center));
flange12Crist[1].SetZ(flange12.Z()-(EXOGAM_Crystal_Center));
flange12Crist[1].SetX(Real_position[flange-1]);
//segment1
flange12CristSeg[1][0].SetY(flange12.Y()-(EXOGAM_Segment_Pos1));
flange12CristSeg[1][0].SetZ(flange12.Z()-(EXOGAM_Segment_Pos1));
flange12CristSeg[1][0].SetX(Real_position[flange-1]);
//segment2
flange12CristSeg[1][1].SetY(flange12.Y()-(EXOGAM_Segment_Pos1));
flange12CristSeg[1][1].SetZ(flange12.Z()-(EXOGAM_Segment_Pos2));
flange12CristSeg[1][1].SetX(Real_position[flange-1]);
//segment3
flange12CristSeg[1][2].SetY(flange12.Y()-(EXOGAM_Segment_Pos2));
flange12CristSeg[1][2].SetZ(flange12.Z()-(EXOGAM_Segment_Pos2));
flange12CristSeg[1][2].SetX(Real_position[flange-1]);
//segment4
flange12CristSeg[1][3].SetY(flange12.Y()-(EXOGAM_Segment_Pos2));
flange12CristSeg[1][3].SetZ(flange12.Z()-(EXOGAM_Segment_Pos1));
flange12CristSeg[1][3].SetX(Real_position[flange-1]);
//Crystal C
flange12Crist[2].SetY(flange12.Y()-(EXOGAM_Crystal_Center));
flange12Crist[2].SetZ(flange12.Z()+(EXOGAM_Crystal_Center));
flange12Crist[2].SetX(Real_position[flange-1]);
//segment1
flange12CristSeg[2][0].SetY(flange12.Y()-(EXOGAM_Segment_Pos1));
flange12CristSeg[2][0].SetZ(flange12.Z()+(EXOGAM_Segment_Pos1));
flange12CristSeg[2][0].SetX(Real_position[flange-1]);
//segment2
flange12CristSeg[2][1].SetY(flange12.Y()-(EXOGAM_Segment_Pos2));
flange12CristSeg[2][1].SetZ(flange12.Z()+(EXOGAM_Segment_Pos1));
flange12CristSeg[2][1].SetX(Real_position[flange-1]);
//segment3
flange12CristSeg[2][2].SetY(flange12.Y()-(EXOGAM_Segment_Pos2));
flange12CristSeg[2][2].SetZ(flange12.Z()+(EXOGAM_Segment_Pos2));
flange12CristSeg[2][2].SetX(Real_position[flange-1]);
//segment4
flange12CristSeg[2][3].SetY(flange12.Y()-(EXOGAM_Segment_Pos1));
flange12CristSeg[2][3].SetZ(flange12.Z()+(EXOGAM_Segment_Pos2));
flange12CristSeg[2][3].SetX(Real_position[flange-1]);
//Crystal D
flange12Crist[3].SetY(flange12.Y()+(EXOGAM_Crystal_Center));
flange12Crist[3].SetZ(flange12.Z()+(EXOGAM_Crystal_Center));
flange12Crist[3].SetX(Real_position[flange-1]);
//segment1
flange12CristSeg[3][0].SetY(flange12.Y()+(EXOGAM_Segment_Pos1));
flange12CristSeg[3][0].SetZ(flange12.Z()+(EXOGAM_Segment_Pos1));
flange12CristSeg[3][0].SetX(Real_position[flange-1]);
//segment2
flange12CristSeg[3][1].SetY(flange12.Y()+(EXOGAM_Segment_Pos1));
flange12CristSeg[3][1].SetZ(flange12.Z()+(EXOGAM_Segment_Pos2));
flange12CristSeg[3][1].SetX(Real_position[flange-1]);
//segment3
flange12CristSeg[3][2].SetY(flange12.Y()+(EXOGAM_Segment_Pos2));
flange12CristSeg[3][2].SetZ(flange12.Z()+(EXOGAM_Segment_Pos2));
flange12CristSeg[3][2].SetX(Real_position[flange-1]);
//segment4
flange12CristSeg[3][3].SetY(flange12.Y()+(EXOGAM_Segment_Pos2));
flange12CristSeg[3][3].SetZ(flange12.Z()+(EXOGAM_Segment_Pos1));
flange12CristSeg[3][3].SetX(Real_position[flange-1]);
if (flange >=1 && flange <=17){
switch(flange){ //which flange ??
case 1:
flange12.RotateY(-45.0*TMath::Pi()/180.0);
for(i=0;i<4;i++){
flange12Crist[i].RotateY(-45.0*TMath::Pi()/180.0);
for(j=0;j<4;j++){
flange12CristSeg[i][j].RotateY(-45.0*TMath::Pi()/180.0);
}
}
break;
case 2:
flange12.RotateY(90.0*TMath::Pi()/180.0);
flange12.RotateX(-45.0*TMath::Pi()/180.0);
for(i=0;i<4;i++){
flange12Crist[i].RotateY(90.0*TMath::Pi()/180.0);
flange12Crist[i].RotateX(-45.0*TMath::Pi()/180.0);
for(j=0;j<4;j++){
flange12CristSeg[i][j].RotateY(90.0*TMath::Pi()/180.0);
flange12CristSeg[i][j].RotateX(-45.0*TMath::Pi()/180.0);
}
}
break;
case 3:
flange12.RotateY(90.0*TMath::Pi()/180.0);
flange12.RotateX(45.0*TMath::Pi()/180.0);
for(i=0;i<4;i++){
flange12Crist[i].RotateY(90.0*TMath::Pi()/180.0);
flange12Crist[i].RotateX(45.0*TMath::Pi()/180.0);
for(j=0;j<4;j++){
flange12CristSeg[i][j].RotateY(90.0*TMath::Pi()/180.0);
flange12CristSeg[i][j].RotateX(45.0*TMath::Pi()/180.0);
}
}
break;
case 4:
flange12.RotateY(45.0*TMath::Pi()/180.0);
for(i=0;i<4;i++){
flange12Crist[i].RotateY(45.0*TMath::Pi()/180.0);
for(j=0;j<4;j++){
flange12CristSeg[i][j].RotateY(45.0*TMath::Pi()/180.0);
}
}
break;
case 5:
flange12.RotateY(135.0*TMath::Pi()/180.0);
for(i=0;i<4;i++){
flange12Crist[i].RotateY(135.0*TMath::Pi()/180.0);
for(j=0;j<4;j++){
flange12CristSeg[i][j].RotateY(135.0*TMath::Pi()/180.0);
}
}
break;
case 6:
flange12.RotateY(90.0*TMath::Pi()/180.0);
flange12.RotateX(90.0*TMath::Pi()/180.0);
for(i=0;i<4;i++){
flange12Crist[i].RotateY(90.0*TMath::Pi()/180.0);
flange12Crist[i].RotateX(90.0*TMath::Pi()/180.0);
for(j=0;j<4;j++){
flange12CristSeg[i][j].RotateY(90.0*TMath::Pi()/180.0);
flange12CristSeg[i][j].RotateX(90.0*TMath::Pi()/180.0);
}
}
break;
case 7:
flange12.RotateY(135.0*TMath::Pi()/180.0);
flange12.RotateX(90.0*TMath::Pi()/180.0);
for(i=0;i<4;i++){
flange12Crist[i].RotateY(135.0*TMath::Pi()/180.0);
flange12Crist[i].RotateX(90.0*TMath::Pi()/180.0);
for(j=0;j<4;j++){
flange12CristSeg[i][j].RotateY(135.0*TMath::Pi()/180.0);
flange12CristSeg[i][j].RotateX(90.0*TMath::Pi()/180.0);
}
}
break;
case 8:
flange12.RotateY(180.0*TMath::Pi()/180.0);
for(i=0;i<4;i++){
flange12Crist[i].RotateY(180.0*TMath::Pi()/180.0);
for(j=0;j<4;j++){
flange12CristSeg[i][j].RotateY(180.0*TMath::Pi()/180.0);
}
}
break;
case 9:
flange12.RotateY(135.0*TMath::Pi()/180.0);
flange12.RotateX(-90.0*TMath::Pi()/180.0);
for(i=0;i<4;i++){
flange12Crist[i].RotateY(135.0*TMath::Pi()/180.0);
flange12Crist[i].RotateX(-90.0*TMath::Pi()/180.0);
for(j=0;j<4;j++){
flange12CristSeg[i][j].RotateY(135.0*TMath::Pi()/180.0);
flange12CristSeg[i][j].RotateX(-90.0*TMath::Pi()/180.0);
}
}
break;
case 10:
flange12.RotateY(90.0*TMath::Pi()/180.0);
flange12.RotateX(-90.0*TMath::Pi()/180.0);
for(i=0;i<4;i++){
flange12Crist[i].RotateY(90.0*TMath::Pi()/180.0);
flange12Crist[i].RotateX(-90.0*TMath::Pi()/180.0);
for(j=0;j<4;j++){
flange12CristSeg[i][j].RotateY(90.0*TMath::Pi()/180.0);
flange12CristSeg[i][j].RotateX(-90.0*TMath::Pi()/180.0);
}
}
break;
case 11:
flange12.RotateY(45.0*TMath::Pi()/180.0);
flange12.RotateX(-90.0*TMath::Pi()/180.0);
for(i=0;i<4;i++){
flange12Crist[i].RotateY(45.0*TMath::Pi()/180.0);
flange12Crist[i].RotateX(-90.0*TMath::Pi()/180.0);
for(j=0;j<4;j++){
flange12CristSeg[i][j].RotateY(45.0*TMath::Pi()/180.0);
flange12CristSeg[i][j].RotateX(-90.0*TMath::Pi()/180.0);
}
}
break;
case 12: //Clover of initialisation --> not move
break;
case 13:
flange12.RotateY(45.0*TMath::Pi()/180.0);
flange12.RotateX(90.0*TMath::Pi()/180.0);
for(i=0;i<4;i++){
flange12Crist[i].RotateY(45.0*TMath::Pi()/180.0);
flange12Crist[i].RotateX(90.0*TMath::Pi()/180.0);
for(j=0;j<4;j++){
flange12CristSeg[i][j].RotateY(45.0*TMath::Pi()/180.0);
flange12CristSeg[i][j].RotateX(90.0*TMath::Pi()/180.0);
}
}
break;
case 14:
flange12.RotateY(-135.0*TMath::Pi()/180.0);
for(i=0;i<4;i++){
flange12Crist[i].RotateY(-135.0*TMath::Pi()/180.0);
for(j=0;j<4;j++){
flange12CristSeg[i][j].RotateY(-135.0*TMath::Pi()/180.0);
}
}
break;
case 15:
flange12.RotateY(90.0*TMath::Pi()/180.0);
flange12.RotateX(-135.0*TMath::Pi()/180.0);
for(i=0;i<4;i++){
flange12Crist[i].RotateY(90.0*TMath::Pi()/180.0);
flange12Crist[i].RotateX(-135.0*TMath::Pi()/180.0);
for(j=0;j<4;j++){
flange12CristSeg[i][j].RotateY(90.0*TMath::Pi()/180.0);
flange12CristSeg[i][j].RotateX(-135.0*TMath::Pi()/180.0);
}
}
break;
case 16:
flange12.RotateY(90.0*TMath::Pi()/180.0);
flange12.RotateX(135.0*TMath::Pi()/180.0);
for(i=0;i<4;i++){
flange12Crist[i].RotateY(90.0*TMath::Pi()/180.0);
flange12Crist[i].RotateX(135.0*TMath::Pi()/180.0);
for(j=0;j<4;j++){
flange12CristSeg[i][j].RotateY(90.0*TMath::Pi()/180.0);
flange12CristSeg[i][j].RotateX(135.0*TMath::Pi()/180.0);
}
}
break;
case 17:
flange12.RotateY(90.0*TMath::Pi()/180.0);
for(i=0;i<4;i++){
flange12Crist[i].RotateY(90.0*TMath::Pi()/180.0);
for(j=0;j<4;j++){
flange12CristSeg[i][j].RotateY(90.0*TMath::Pi()/180.0);
}
}
break;
default :
break;
}
// Take into account the case of a target which is not at the center of the EXOGAM array
// Target position
TVector3 targetPos(TARGET_POSITION_X, TARGET_POSITION_Y, TARGET_POSITION_Z);
// case of the selected EXOGAM detector
/* for (Int_t k = 0; k < 3; k++) {
cout << "flange12 avant : " << flange12(k) << endl;
cout << "targetPos : " << targetPos(k) << endl;
}*/
flange12 -= targetPos;
/* for (Int_t k = 0; k < 3; k++) {
cout << "flange12 apres : " << flange12(k) << endl;
}*/
// loop on cristals
for (Int_t ii = 0; ii < 4; ii++) {
flange12Crist[ii] -= targetPos;
// loop on segments
for (Int_t jj = 0; jj < 4; jj++) {
flange12CristSeg[ii][jj] -= targetPos;
}
}
#if COMMENTFLAG
printf(" flange %d theta %f phi Proj %f \n",flange,flange12.Theta()*180.0/(TMath::Pi()),flange12.Phi()*180.0/(TMath::Pi()));
for(i=0;i<4;i++){
printf(" flange %d cristal %d theta %f phi Proj %f \n",
flange,i+1,flange12Crist[i].Theta()*180.0/(TMath::Pi()),flange12Crist[i].Phi()*180.0/(TMath::Pi()));
for(j=0;j<4;j++){
printf(" flange %d cristal %d seg %d theta %f phi Proj %f \n",
flange,i+1,j+1,flange12CristSeg[i][j].Theta()*180.0/(TMath::Pi()),flange12CristSeg[i][j].Phi()*180.0/(TMath::Pi()));
}
}
printf("#########################\n");
#endif
//Output result in rad in the "Clover_struc" structur
Result.Theta_Clover=flange12.Theta();
Result.Phi_Clover=flange12.Phi();
Result.X_Clover=flange12.X();
Result.Y_Clover=flange12.Y();
Result.Z_Clover=flange12.Z();
for(i=0;i<4;i++){
Result.Theta_Crystal[i]=flange12Crist[i].Theta();
Result.Phi_Crystal[i]=flange12Crist[i].Phi();
Result.X_Crystal[i]=flange12Crist[i].X();
Result.Y_Crystal[i]=flange12Crist[i].Y();
Result.Z_Crystal[i]=flange12Crist[i].Z();
for(j=0;j<4;j++){
Result.Theta_Crystal_Seg[i][j]=flange12CristSeg[i][j].Theta();
Result.Phi_Crystal_Seg[i][j]=flange12CristSeg[i][j].Phi();
Result.X_Crystal_Seg[i][j]=flange12CristSeg[i][j].X();
Result.Y_Crystal_Seg[i][j]=flange12CristSeg[i][j].Y();
Result.Z_Crystal_Seg[i][j]=flange12CristSeg[i][j].Z();
}
}
}
else {printf("Bad flange number, flange %d doesn't exist in EXOGAM !!!!!! \n",flange);}
return Result;
}
// Routine of doppler correction
float Doppler_Correction(float Theta_Gamma, float Phi_Gamma, float Theta_Part, float Phi_Part, float Beta_Part, float energie_Mes){ //rad, v/c
float energievraie,cosinusPSI;
cosinusPSI =TMath::Sin(Theta_Part)*TMath::Cos(Phi_Part)*TMath::Sin(Theta_Gamma)*TMath::Cos(Phi_Gamma)+
TMath::Sin(Theta_Part)*TMath::Sin(Phi_Part)*TMath::Sin(Theta_Gamma)*TMath::Sin(Phi_Gamma)+
TMath::Cos(Theta_Part)*TMath::Cos(Theta_Gamma);
energievraie = energie_Mes*(1.-Beta_Part*cosinusPSI)/sqrt(1.-Beta_Part*Beta_Part);
return energievraie;
}
NPLib/Detectors/Exogam/TExogamData.h
+ 16
− 16
View file @ 8cf70b25
......@@ -52,10 +52,10 @@ class TExogamData : public TObject {
///////////////////// SETTERS ////////////////////////
inline void SetExo(unsigned int& Crystal, unsigned int& Energy,
unsigned int& Energy_HG, unsigned long long& TS, unsigned int& TDC,
unsigned int& BGO, unsigned int& CsI, unsigned int& Outer1,
unsigned int& Outer2, unsigned int& Outer3,unsigned int& Outer4) {
inline void SetExo(const unsigned int& Crystal,const unsigned int& Energy,
const unsigned int& Energy_HG,const unsigned long long& TS,const unsigned int& TDC,
const unsigned int& BGO,const unsigned int& CsI,const unsigned int& Outer1,
const unsigned int& Outer2,const unsigned int& Outer3,const unsigned int& Outer4) {
fExo_Crystal.push_back(Crystal);
fExo_E.push_back(Energy);
fExo_E_HG.push_back(Energy_HG);
......@@ -69,18 +69,18 @@ class TExogamData : public TObject {
fExo_Outer4.push_back(Outer4);
}
///////////////////// GETTERS ////////////////////////
inline unsigned int GetExoMult(unsigned int& i) { return fExo_Crystal[i]; }
inline unsigned int GetExoCrystal(unsigned int& i) { return fExo_Crystal[i]; }
inline unsigned int GetExoE(unsigned int& i) { return fExo_E[i]; }
inline unsigned int GetExoEHG(unsigned int& i) { return fExo_E_HG[i]; }
inline unsigned long long GetExoTS(unsigned int& i) { return fExo_TS[i]; }
inline unsigned int GetExoTDC(unsigned int& i) { return fExo_TDC[i]; }
inline unsigned int GetExoBGO(unsigned int& i) { return fExo_BGO[i]; }
inline unsigned int GetExoCsI(unsigned int& i) { return fExo_CsI[i]; }
inline unsigned int GetExoOuter1(unsigned int& i) { return fExo_Outer1[i]; }
inline unsigned int GetExoOuter2(unsigned int& i) { return fExo_Outer2[i]; }
inline unsigned int GetExoOuter3(unsigned int& i) { return fExo_Outer3[i]; }
inline unsigned int GetExoOuter4(unsigned int& i) { return fExo_Outer4[i]; }
inline unsigned int GetExoMult() { return fExo_Crystal.size(); }
inline unsigned int GetExoCrystal(const unsigned int& i) const { return fExo_Crystal[i]; }
inline unsigned int GetExoE(const unsigned int& i) const { return fExo_E[i]; }
inline unsigned int GetExoEHG(const unsigned int& i) const { return fExo_E_HG[i]; }
inline unsigned long long GetExoTS(const unsigned int& i) const { return fExo_TS[i]; }
inline unsigned int GetExoTDC(const unsigned int& i) const { return fExo_TDC[i]; }
inline unsigned int GetExoBGO(const unsigned int& i) const { return fExo_BGO[i]; }
inline unsigned int GetExoCsI(const unsigned int& i) const { return fExo_CsI[i]; }
inline unsigned int GetExoOuter1(const unsigned int& i) const { return fExo_Outer1[i]; }
inline unsigned int GetExoOuter2(const unsigned int& i) const { return fExo_Outer2[i]; }
inline unsigned int GetExoOuter3(const unsigned int& i) const { return fExo_Outer3[i]; }
inline unsigned int GetExoOuter4(const unsigned int& i) const { return fExo_Outer4[i]; }
ClassDef(TExogamData, 1) // ExogamData structure
};
......
NPLib/Detectors/Exogam/TExogamPhysics.cxx
+ 206
− 96
View file @ 8cf70b25
......@@ -26,6 +26,7 @@ using namespace EXOGAM_LOCAL;
#include <iostream>
#include <sstream>
#include <stdlib.h>
#include <functional>
// NPL
#include "NPDetectorFactory.h"
......@@ -47,6 +48,11 @@ ClassImp(TExogamPhysics)
NumberOfHitCristal = 0;
m_Spectra = NULL;
NumberOfClover = 0;
m_EXO_E_RAW_Threshold = 0;
m_EXO_E_Threshold = 0;
m_EXO_EHG_RAW_Threshold = 0;
m_EXO_TDC_RAW_Threshold = 0;
m_EXO_TDC_RAW_Threshold = 0;
m_PreTreatedData = new TExogamData;
m_EventData = new TExogamData;
......@@ -59,112 +65,59 @@ ClassImp(TExogamPhysics)
void TExogamPhysics::BuildSimplePhysicalEvent() { BuildPhysicalEvent(); }
///////////////////////////////////////////////////////////////////////////
void TExogamPhysics::PreTreat() {
/*ClearPreTreatedData();
// Clearing PreTreat TExogamData
ClearPreTreatedData();
// Clearing local variables for pretreat
ResetPreTreatVariable();
//E
for(unsigned int i = 0 ; i < EventData -> GetECCEMult(); i++) {
UShort_t cristal_E = 10000 ; UShort_t cristal_T = 2000;
//if(IsValidChannel)
{
int clover = EventData -> GetECCEClover(i);
int cristal = EventData -> GetECCECristal(i);
if(EventData -> GetECCEEnergy(i) < 3000) cristal_E = CalibrationManager::getInstance()->
ApplyCalibration("EXOGAM/Cl"+ NPL::itoa(clover)+"_Cr"+ NPL::itoa(cristal)+"_Elow", EventData -> GetECCEEnergy(i));
else cristal_E = CalibrationManager::getInstance()->
ApplyCalibration("EXOGAM/Cl"+ NPL::itoa(clover)+"_Cr"+ NPL::itoa(cristal)+"_Ehigh", EventData -> GetECCEEnergy(i));
if(cristal_E > Threshold_ECC)
{
PreTreatedData->SetECCEClover ( clover ) ;
PreTreatedData->SetECCECristal( cristal ) ;
PreTreatedData->SetECCEEnergy ( cristal_E ) ;
bool checkT = false;
for(unsigned int k = 0; k < EventData -> GetECCTMult(); k++){
if(clover == EventData -> GetECCTClover(k) && cristal == EventData -> GetECCTCristal(k)){
// cout << EventData -> GetECCTTime(k) << endl;
if(EventData -> GetECCTTime(k) < 16383) cristal_T = CalibrationManager::getInstance()->
ApplyCalibration("EXOGAM/Cl"+ NPL::itoa(clover)+"_Cr"+ NPL::itoa(cristal)+"_T", EventData -> GetECCTTime(k)); else
cristal_T = 2500;
//if(cristal_T >5000 && cristal_T !=25000 ) cout << "PreTreat " << cristal_T << " " << EventData ->
GetECCTTime(k) << " " << clover << " " << cristal << " " << EventData->GetECCTMult() << endl;
checkT=true;
PreTreatedData->SetECCTClover (clover ) ;
PreTreatedData->SetECCTCristal( cristal ) ;
PreTreatedData->SetECCTTime ( cristal_T ) ;
ECC_Multiplicity ++;
GOCCE_Multiplicity++;
}
}
if(!checkT) {
PreTreatedData->SetECCTClover (clover ) ;
PreTreatedData->SetECCTCristal( cristal ) ;
PreTreatedData->SetECCTTime ( -1000 ) ;
}
}
}
}
//cout << PreTreatedData-> GetECCTMult() << " " << PreTreatedData-> GetECCEMult() << endl;
//GOCCE
//E
for(unsigned int i = 0 ; i < EventData -> GetGOCCEEMult(); i++) {
UShort_t segment_E = 25000;
//if(IsValidChannel)
{
int clover = EventData -> GetGOCCEEClover(i);
int cristal = EventData -> GetGOCCEECristal(i);
int segment = EventData -> GetGOCCEESegment(i);
if(EventData -> GetGOCCEEEnergy(i) > RawThreshold_GOCCE)
{
segment_E = CalibrationManager::getInstance()->ApplyCalibration("EXOGAM/Cl"+ NPL::itoa(clover)+"_Cr"+
NPL::itoa(cristal)+"_Seg"+ NPL::itoa(segment)+"_E", EventData -> GetGOCCEEEnergy(i));
if(segment_E > Threshold_GOCCE)
{
PreTreatedData->SetGOCCEEClover ( clover ) ;
PreTreatedData->SetGOCCEECristal( cristal ) ;
PreTreatedData->SetGOCCEESegment( segment ) ;
PreTreatedData->SetGOCCEEEnergy ( segment_E ) ;
}
}
else
{
m_EXO_Mult = m_EventData->GetExoMult();
for (unsigned int i = 0; i < m_EXO_Mult; ++i) {
if (m_EventData->GetExoE(i) > m_EXO_E_RAW_Threshold)
EXO_E = fEXO_E(m_EventData, i);
if (m_EventData->GetExoEHG(i) > m_EXO_EHG_RAW_Threshold)
EXO_EHG = fEXO_EHG(m_EventData, i);
if (m_EventData->GetExoTDC(i) > m_EXO_TDC_RAW_Threshold)
EXO_TDC = fEXO_T(m_EventData, i);
EXO_Outer1 = fEXO_Outer(m_EventData, i, 1);
EXO_Outer2 = fEXO_Outer(m_EventData, i, 2);
EXO_Outer3 = fEXO_Outer(m_EventData, i, 3);
EXO_Outer4 = fEXO_Outer(m_EventData, i, 4);
if(EXO_E > m_EXO_E_Threshold){
m_PreTreatedData->SetExo(m_EventData->GetExoCrystal(i), EXO_E,
EXO_EHG, m_EventData->GetExoTS(i), EXO_TDC,
m_EventData->GetExoBGO(i), m_EventData->GetExoCsI(i), EXO_Outer1,
EXO_Outer2, EXO_Outer3, EXO_Outer4);
}
}
}
}
//cout << "EXOGAM pretreat ok!" << endl;
return;
*/
}
///////////////////////////////////////////////////////////////////////////
void TExogamPhysics::ResetPreTreatVariable(){
EXO_E = -1000;
EXO_EHG = -1000;
EXO_TDC = -1000;
EXO_Outer1 = -1000;
EXO_Outer2 = -1000;
EXO_Outer3 = -1000;
EXO_Outer4 = -1000;
}
void TExogamPhysics::BuildPhysicalEvent() {
if (NPOptionManager::getInstance()->IsReader() == true) {
m_EventData = &(**r_ReaderEventData);
}
PreTreat();
//PreTreat();
//for(unsigned int i = 0; i < m_PreTreatedData->GetExoMult(); i++){
// mean_free_path = ComputeMeanFreePath(m_PreTreatedData->GetExoE(i));
//}
std::cout << ComputeMeanFreePath(7000) << std::endl;
/*
if(PreTreatedData -> GetECCEMult() != PreTreatedData -> GetECCTMult()) cout << PreTreatedData -> GetECCEMult() << " "
<< PreTreatedData -> GetECCTMult() << endl;
......@@ -396,6 +349,24 @@ void TExogamPhysics::BuildPhysicalEvent() {
*/
}
double TExogamPhysics::ComputeMeanFreePath(double Energy){
auto b = Map_PhotonCS.lower_bound(Energy);
auto a = prev(b);
if(b == Map_PhotonCS.begin()){
a = b;
b++;
}
else if(b == Map_PhotonCS.end()){
b--;
a = prev(b);
}
std::cout << a->first << " " << b->first << " " << a->second << " " << b->second << std::endl;
double coeff = (Energy - a->first)/(b->first - a->first);
double PhotonCrossSection = a->second + coeff*(b->second - a->second); // mm2/g
return 1./(GeDensity*PhotonCrossSection);
}
double TExogamPhysics::DopplerCorrection(double E, double Theta) {
double Pi = 3.141592654;
TString filename = "configs/beta.txt";
......@@ -475,8 +446,77 @@ void TExogamPhysics::ReadConfiguration(NPL::InputParser parser) {
exit(1);
}
}
ReadAnalysisConfig();
}
void TExogamPhysics::ReadAnalysisConfig() {
bool ReadingStatus = false;
// path to photon cross section
string CSFilename = "../../Inputs/PhotonCrossSection/CoherentGe.xcom";
string LineBuffer;
ifstream CSFile;
CSFile.open(CSFilename.c_str());
if (!CSFile.is_open()) {
cout << " No CS file found "
<< CSFilename << endl;
return;
}
while(CSFile.good()){
double gammaE, CrossSection;
getline(CSFile, LineBuffer);
istringstream ss(LineBuffer);
ss >> gammaE >> CrossSection; // E in MeV, converted to keV, CrossSection in cm2/g
gammaE *= 1000.; // Convertion to keV
CrossSection *= 100.;
Map_PhotonCS[gammaE] = CrossSection;
}
// path to file
string FileName = "./configs/ConfigExogam.dat";
// open analysis config file
ifstream AnalysisConfigFile;
AnalysisConfigFile.open(FileName.c_str());
if (!AnalysisConfigFile.is_open()) {
cout << " No ConfigExogam.dat found: Default parameters loaded for "
"Analysis "
<< FileName << endl;
return;
}
string DataBuffer, whatToDo;
while (!AnalysisConfigFile.eof()) {
// Pick-up next line
getline(AnalysisConfigFile, LineBuffer);
// search for "header"
if (LineBuffer.compare(0, 11, "ConfigExogam") == 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 == "EXO_Threshold") {
//AnalysisConfigFile >> DataBuffer;
//m_MaximumStripMultiplicityAllowed = atoi(DataBuffer.c_str());
//cout << "MAXIMUN STRIP MULTIPLICITY " << m_MaximumStripMultiplicityAllowed << endl;
}
}
}
}
///////////////////////////////////////////////////////////////////////////
void TExogamPhysics::InitSpectra() { m_Spectra = new TExogamSpectra(NumberOfClover); }
......@@ -587,6 +627,14 @@ void TExogamPhysics::AddParameterToCalibrationManager() {
// Activated associated Branches and link it to the private member DetectorData address
// In this method mother Branches (Detector) AND daughter leaf (fDetector_parameter) have to be activated
void TExogamPhysics::InitializeRootInputRaw() {
TChain* inputChain = RootInput::getInstance()->GetChain();
// Option to use the nptreereader anaysis
if (NPOptionManager::getInstance()->IsReader() == true) {
TTreeReader* inputTreeReader = RootInput::getInstance()->GetTreeReader();
inputTreeReader->SetTree(inputChain);
}
// Option to use the standard npanalysis
else{
TChain* inputChain = RootInput::getInstance()->GetChain();
inputChain->SetBranchStatus("EXOGAM", true);
inputChain->SetBranchStatus("fEXO_*", true);
......@@ -599,12 +647,21 @@ void TExogamPhysics::InitializeRootInputRaw() {
cristal_mult = new TH1F("cristal_mult","cristal_mult",20,0,20);
outputList->Add(cristal_mult);
*/
}
}
/////////////////////////////////////////////////////////////////////
// Activated associated Branches and link it to the private member DetectorPhysics address
// In this method mother Branches (Detector) AND daughter leaf (parameter) have to be activated
void TExogamPhysics::InitializeRootInputPhysics() {
TChain* inputChain = RootInput::getInstance()->GetChain();
// Option to use the nptreereader anaysis
if (NPOptionManager::getInstance()->IsReader() == true) {
TTreeReader* inputTreeReader = RootInput::getInstance()->GetTreeReader();
inputTreeReader->SetTree(inputChain);
}
// Option to use the standard npanalysis
else{
TChain* inputChain = RootInput::getInstance()->GetChain();
inputChain->SetBranchStatus("EventMultiplicty", true);
inputChain->SetBranchStatus("ECC_Multiplicity", true);
......@@ -627,6 +684,7 @@ void TExogamPhysics::InitializeRootInputPhysics() {
inputChain->SetBranchStatus("Position", true);
inputChain->SetBranchStatus("Theta", true);
inputChain->SetBranchAddress("EXOGAM", &m_EventPhysics);
}
}
/////////////////////////////////////////////////////////////////////
......@@ -644,9 +702,58 @@ void TExogamPhysics::InitializeRootOutput() {
*/
}
void TExogamPhysics::SetTreeReader(TTreeReader* TreeReader) {
TExogamPhysicsReader::r_SetTreeReader(TreeReader);
}
///////////////////////////////////////////////////////////////////////////
namespace EXOGAM_LOCAL {
// tranform an integer to a string
double fEXO_E(const TExogamData* m_EventData, const unsigned int& i) {
static string name;
name = "EXOGAM/Cr_";
name += NPL::itoa(m_EventData->GetExoE(i));
name += "_E";
return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoE(i));
}
double fEXO_EHG(const TExogamData* m_EventData, const unsigned int& i) {
static string name;
name = "EXOGAM/Cr_";
name += NPL::itoa(m_EventData->GetExoEHG(i));
name += "_EHG";
return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoEHG(i));
}
double fEXO_T(const TExogamData* m_EventData, const unsigned int& i) {
static string name;
name = "EXOGAM/Cr_";
name += NPL::itoa(m_EventData->GetExoTDC(i));
name += "_TDC";
return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoTDC(i));
}
double fEXO_Outer(const TExogamData* m_EventData, const unsigned int& i, const unsigned int OuterNumber) {
static string name;
name = "EXOGAM/Cr_";
name += NPL::itoa(m_EventData->GetExoE(i));
name += "_Outer";
name += NPL::itoa(OuterNumber);
name += "_E";
if(OuterNumber == 1)
return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoOuter1(i));
else if(OuterNumber == 2)
return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoOuter2(i));
else if(OuterNumber == 3)
return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoOuter3(i));
else if(OuterNumber == 4)
return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoOuter4(i));
else{
std::cout << "WARNING: Outer number != 1-4, something is wrong\n";
return 0;
};
}
string itoa(int value) {
std::ostringstream o;
......@@ -663,6 +770,8 @@ namespace EXOGAM_LOCAL {
////////////////////////////////////////////////////////////////////////////////
NPL::VDetector* TExogamPhysics::Construct() { return (NPL::VDetector*)new TExogamPhysics(); }
NPL::VTreeReader* TExogamPhysics::ConstructReader() { return (NPL::VTreeReader*)new TExogamPhysicsReader(); }
////////////////////////////////////////////////////////////////////////////////
// Registering the construct method to the factory //
////////////////////////////////////////////////////////////////////////////////
......@@ -672,6 +781,7 @@ class proxy_exogam {
proxy_exogam() {
NPL::DetectorFactory::getInstance()->AddToken("Exogam", "Exogam");
NPL::DetectorFactory::getInstance()->AddDetector("Exogam", TExogamPhysics::Construct);
NPL::DetectorFactory::getInstance()->AddDetectorReader("Exogam", TExogamPhysics::ConstructReader);
}
};
......
NPLib/Detectors/Exogam/TExogamPhysics.h
+ 31
− 1
View file @ 8cf70b25
......@@ -29,6 +29,7 @@
// NPL
#include "NPCalibrationManager.h"
#include "NPVDetector.h"
#include "NPVTreeReader.h"
#include "TExogamData.h"
#include "TExogamSpectra.h"
#include "TExogamPhysicsReader.h"
......@@ -144,6 +145,8 @@ class TExogamPhysics : public TObject, public NPL::VDetector, public TExogamPhys
void CheckSpectra();
// Used for Online only, clear all the spectra hold by the Spectra class
void ClearSpectra();
void SetTreeReader(TTreeReader* TreeReader);
private: // Root Input and Output tree classes
......@@ -170,11 +173,32 @@ class TExogamPhysics : public TObject, public NPL::VDetector, public TExogamPhys
// Retrieve raw and pre-treated data
TExogamData* GetRawData() const {return m_EventData;}
TExogamData* GetPreTreatedData() const {return m_PreTreatedData;}
void ResetPreTreatVariable();
void ReadAnalysisConfig();
double ComputeMeanFreePath(double GammaEnergy);
private: // Variables for analysis
unsigned int m_EXO_Mult;
double m_EXO_E_RAW_Threshold;
double m_EXO_E_Threshold;
double m_EXO_EHG_RAW_Threshold;
double m_EXO_TDC_RAW_Threshold;
double EXO_E;
double EXO_EHG;
double EXO_TDC;
double EXO_Outer1;
double EXO_Outer2;
double EXO_Outer3;
double EXO_Outer4;
double mean_free_path;
const double GeDensity = 0.005323; //! g/mm3
std::map<double, double> Map_PhotonCS;
private: // Spectra Class
TExogamSpectra* m_Spectra;//!
......@@ -184,11 +208,17 @@ class TExogamPhysics : public TObject, public NPL::VDetector, public TExogamPhys
public: // Static constructor to be passed to the Detector Factory
static NPL::VDetector* Construct();
static NPL::VTreeReader* ConstructReader();
ClassDef(TExogamPhysics,1) // ExogamPhysics structure
};
namespace EXOGAM_LOCAL
{
double fEXO_E(const TExogamData* m_EventData, const unsigned int& i);
double fEXO_EHG(const TExogamData* m_EventData, const unsigned int& i);
double fEXO_T(const TExogamData* m_EventData, const unsigned int& i);
double fEXO_Outer(const TExogamData* m_EventData, const unsigned int& i, const unsigned int OuterNumber);
const double Threshold_ECC = 50;
const double Threshold_GOCCE = 0;
const double RawThreshold_ECC = 0;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment