diff --git a/NPLib/Detectors/Exogam/Geometry_Clover_Exogam.h b/NPLib/Detectors/Exogam/Geometry_Clover_Exogam.h new file mode 100644 index 0000000000000000000000000000000000000000..256d54f88b4deafe27f76d7608d7fca337093d67 --- /dev/null +++ b/NPLib/Detectors/Exogam/Geometry_Clover_Exogam.h @@ -0,0 +1,577 @@ +/************************************************************************************************ + + +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; +} diff --git a/NPLib/Detectors/Exogam/TExogamData.h b/NPLib/Detectors/Exogam/TExogamData.h index c5d377ba575788f0cb29d59061b2d07fd8d096f6..e29da96ce2d0a50799250b1f29ef8ba96f3d871a 100644 --- a/NPLib/Detectors/Exogam/TExogamData.h +++ b/NPLib/Detectors/Exogam/TExogamData.h @@ -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 }; diff --git a/NPLib/Detectors/Exogam/TExogamPhysics.cxx b/NPLib/Detectors/Exogam/TExogamPhysics.cxx index 7575a007d15a70e8031b473a55d409141e74973b..ddc738cffce85fce65f9afc9c618a618f4b29160 100644 --- a/NPLib/Detectors/Exogam/TExogamPhysics.cxx +++ b/NPLib/Detectors/Exogam/TExogamPhysics.cxx @@ -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); } }; diff --git a/NPLib/Detectors/Exogam/TExogamPhysics.h b/NPLib/Detectors/Exogam/TExogamPhysics.h index ced6f1e29e0f577fdd9795543b333e06d0626590..28c818dfcfb01366edeca5f6e3ddd65a66ab274b 100644 --- a/NPLib/Detectors/Exogam/TExogamPhysics.h +++ b/NPLib/Detectors/Exogam/TExogamPhysics.h @@ -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;