diff --git a/NPLib/Detectors/Exogam/TExogamPhysics.cxx b/NPLib/Detectors/Exogam/TExogamPhysics.cxx
index 86cb043e0f44dcec5f8f21c4fbb5c1f6e1108ee8..602bbed1432ac582970e90c4fc8a444442dd35db 100644
--- a/NPLib/Detectors/Exogam/TExogamPhysics.cxx
+++ b/NPLib/Detectors/Exogam/TExogamPhysics.cxx
@@ -48,6 +48,7 @@ ClassImp(TExogamPhysics)
NumberOfHitCristal = 0;
m_Spectra = NULL;
NumberOfClover = 0;
+ m_EXO_OuterUp_RAW_Threshold = 4294966296;
m_EXO_E_RAW_Threshold = 0;
m_EXO_E_Threshold = 0;
m_EXO_EHG_RAW_Threshold = 0;
@@ -67,34 +68,54 @@ void TExogamPhysics::BuildSimplePhysicalEvent() { BuildPhysicalEvent(); }
void TExogamPhysics::PreTreat() {
// Clearing PreTreat TExogamData
ClearPreTreatedData();
- // Clearing local variables for pretreat
//E
m_EXO_Mult = m_EventData->GetExoMult();
for (unsigned int i = 0; i < m_EXO_Mult; ++i) {
- ResetPreTreatVariable();
- // std::cout << "test1 " << EXO_E << std::endl;
+
+ ResetPreTreatVariable();
+
if (m_EventData->GetExoE(i) > m_EXO_E_RAW_Threshold)
EXO_E = fEXO_E(m_EventData, i);
- // std::cout << "test2 " << EXO_E << std::endl;
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);
+
+ if (m_EventData->GetExoOuter1(i) < m_EXO_OuterUp_RAW_Threshold)
+ EXO_Outer1 = fEXO_Outer(m_EventData, i, 0);
+ else
+ EXO_Outer1 = 0;
- EXO_Outer1 = fEXO_Outer(m_EventData, i, 0);
- EXO_Outer2 = fEXO_Outer(m_EventData, i, 1);
- EXO_Outer3 = fEXO_Outer(m_EventData, i, 2);
- EXO_Outer4 = fEXO_Outer(m_EventData, i, 3);
-
+ if (m_EventData->GetExoOuter2(i) < m_EXO_OuterUp_RAW_Threshold)
+ EXO_Outer2 = fEXO_Outer(m_EventData, i, 1);
+ else
+ EXO_Outer2 = 0;
+
+ if (m_EventData->GetExoOuter3(i) < m_EXO_OuterUp_RAW_Threshold)
+ EXO_Outer3 = fEXO_Outer(m_EventData, i, 2);
+ else
+ EXO_Outer3 = 0;
+
+ if (m_EventData->GetExoOuter4(i) < m_EXO_OuterUp_RAW_Threshold)
+ EXO_Outer4 = fEXO_Outer(m_EventData, i, 3);
+ else
+ EXO_Outer4 = 0;
+ // std::cout << m_EventData->GetExoOuter4(i) << std::endl;
+ // std::cout << EXO_E << " " << EXO_EHG << " " << EXO_Outer1 << " " << EXO_Outer2 << " " << EXO_Outer3 << " " << EXO_Outer4 << std::endl;
+ //EXO_Outer2 = fEXO_Outer(m_EventData, i, 1);
+ //EXO_Outer3 = fEXO_Outer(m_EventData, i, 2);
+ //EXO_Outer4 = fEXO_Outer(m_EventData, i, 3);
+
+ // *1000 to convert MeV into keV
if(EXO_E > m_EXO_E_Threshold){
m_PreTreatedData->SetExo(m_EventData->GetExoCrystal(i), EXO_E*1000,
EXO_EHG*1000, m_EventData->GetExoTS(i), EXO_TDC,
- m_EventData->GetExoBGO(i), m_EventData->GetExoCsI(i), EXO_Outer1,
- EXO_Outer2, EXO_Outer3, EXO_Outer4);
+ m_EventData->GetExoBGO(i), m_EventData->GetExoCsI(i), EXO_Outer1*1000,
+ EXO_Outer2*1000, EXO_Outer3*1000, EXO_Outer4*1000);
}
}
}
@@ -113,274 +134,96 @@ void TExogamPhysics::BuildPhysicalEvent() {
if (NPOptionManager::getInstance()->IsReader() == true) {
m_EventData = &(**r_ReaderEventData);
}
+ // std::cout << m_EventData << std::endl;
PreTreat();
+ // This maps stores ID of events sorted by flange number. Map key is flange nbr, vector should contain ID of events
+ std::map<unsigned int,std::vector<unsigned int>> HitsID;
+
for(unsigned int i = 0; i < m_PreTreatedData->GetExoMult(); i++){
-
- mean_free_path = ComputeMeanFreePath(m_PreTreatedData->GetExoE(i));
+
+ // Doing flange and crystal matching
flange_nbr = MapCrystalFlangeCLover[m_PreTreatedData->GetExoCrystal(i)].first;
crystal_nbr = MapCrystalFlangeCLover[m_PreTreatedData->GetExoCrystal(i)].second;
-
- EnergyDoppler = 0;
- double MaxOuter = 0;
- unsigned int MaxOuter_Id;
- if(m_PreTreatedData->GetExoOuter1(i) > MaxOuter){
- MaxOuter = m_PreTreatedData->GetExoOuter1(i);
- MaxOuter_Id = 0;
- }
- if(m_PreTreatedData->GetExoOuter2(i) > MaxOuter){
- MaxOuter = m_PreTreatedData->GetExoOuter2(i);
- MaxOuter_Id = 1;
- }
- if(m_PreTreatedData->GetExoOuter3(i) > MaxOuter){
- MaxOuter = m_PreTreatedData->GetExoOuter3(i);
- MaxOuter_Id = 2;
- }
- if(m_PreTreatedData->GetExoOuter4(i) > MaxOuter){
- MaxOuter = m_PreTreatedData->GetExoOuter4(i);
- MaxOuter_Id = 3;
- }
- if(MaxOuter > 0){
- Exogam_struc = Ask_For_Angles(flange_nbr, mean_free_path);
- double Theta_seg = Exogam_struc.Theta_Crystal_Seg[crystal_nbr][MaxOuter_Id];
- double Phi_seg = Exogam_struc.Phi_Crystal_Seg[crystal_nbr][MaxOuter_Id];
- EnergyDoppler = Doppler_Correction(Theta_seg,Phi_seg,0,0,Beta,m_PreTreatedData->GetExoE(i));
- }
E_Cal.push_back(m_PreTreatedData->GetExoE(i));
- Flange_N.push_back(flange_nbr);
- Crystal_N.push_back(crystal_nbr);
- E_Doppler.push_back(EnergyDoppler);
+ EHG_Cal.push_back(m_PreTreatedData->GetExoEHG(i));
+ Outer1_Cal.push_back(m_PreTreatedData->GetExoOuter1(i));
+ Outer2_Cal.push_back(m_PreTreatedData->GetExoOuter2(i));
+ Outer3_Cal.push_back(m_PreTreatedData->GetExoOuter3(i));
+ Outer4_Cal.push_back(m_PreTreatedData->GetExoOuter4(i));
+ FlangeN.push_back(flange_nbr);
+ CrystalN.push_back(crystal_nbr);
+
+ // Filling HitsID
+ // std::cout << i << " " << flange_nbr << " " << crystal_nbr << std::endl;
+
+ HitsID[flange_nbr].push_back(i);
}
-/*
- if(PreTreatedData -> GetECCEMult() != PreTreatedData -> GetECCTMult()) cout << PreTreatedData -> GetECCEMult() << " "
- << PreTreatedData -> GetECCTMult() << endl;
-
- for(unsigned int i = 0 ; i < PreTreatedData -> GetECCEMult(); i++) {
-
- // cout << i << " " << cristal_E << endl;
- // if(PreTreatedData->GetECCTTime(i) > 0)
- {
- ECC_E.push_back(PreTreatedData->GetECCEEnergy(i));
- ECC_T.push_back(PreTreatedData->GetECCTTime(i));
- ECC_CloverNumber.push_back(PreTreatedData->GetECCEClover(i));
- ECC_CristalNumber.push_back(PreTreatedData->GetECCECristal(i));
-
- // cout << "BuildPhys " << PreTreatedData->GetECCEClover(i) << " " << PreTreatedData->GetECCECristal(i)<< " " <<
- PreTreatedData->GetECCTTime(i) << " " << endl;
+ // Now that HitsID is full, we use it to process simple AddBack of events in the same flange
+ // Basically looping on all flanges, then on al events ID in each flange
+ for(auto it = HitsID.begin(); it != HitsID.end(); it++){
+ double E_AddBack = 0;
+ double E_Max = 0;
+ unsigned int Id_Max = 0;
+ for(auto itvec = (*it).second.begin(); itvec !=(*it).second.end(); itvec++){
+ E_AddBack+= m_PreTreatedData->GetExoE(*itvec);
+ if(E_Max < m_PreTreatedData->GetExoE(*itvec)){
+ E_Max = m_PreTreatedData->GetExoE(*itvec);
+ Id_Max = *itvec;
}
+ }
+ // Doing it again for this loop, it's a bit unhappy but didnt find a better way to do it yet
+ flange_nbr = (*it).first;
+ crystal_nbr = MapCrystalFlangeCLover[m_PreTreatedData->GetExoCrystal(Id_Max)].second;
+
+ // Adding all AddBack (AB) related stuff
+ E_AB.push_back(E_AddBack);
+ FlangeN_AB.push_back(flange_nbr);
+ Size_AB.push_back((*it).second.size());
+
+ // Adding these parameters for Doppler correction purposes (D)
+ CrystalN_ABD.push_back(crystal_nbr);
+ int MaxOuterId = GetMaxOuter(Id_Max);
+ OuterN_ABD.push_back(GetMaxOuter(Id_Max));
+
+ // If a max Outer is found, Do doppler correction, else push_back -1000;
+ double EnergyDoppler = -1000;
+ if(MaxOuterId > -1){
+ EnergyDoppler= GetDoppler(E_AddBack, flange_nbr, crystal_nbr, MaxOuterId);
+ }
+ E_ABD.push_back(EnergyDoppler);
}
+}
-
- for(unsigned int j = 0 ; j < PreTreatedData -> GetGOCCEEMult(); j++) {
- GOCCE_E.push_back(PreTreatedData->GetGOCCEEEnergy(j));
- GOCCE_CloverNumber.push_back(PreTreatedData->GetGOCCEEClover(j));
- GOCCE_CristalNumber.push_back(PreTreatedData->GetGOCCEECristal(j));
- GOCCE_SegmentNumber.push_back(PreTreatedData->GetGOCCEESegment(j));
+int TExogamPhysics::GetMaxOuter(unsigned int EventId){
+ double OuterMax = 0;
+ int OuterId = -1;
+ if(m_EventData->GetExoOuter1(EventId) > OuterMax){
+ OuterMax = m_EventData->GetExoOuter1(EventId);
+ OuterId = 0;
}
-
-
- //int NumberOfHitClover = 0;
-
- int DetectorID = -1;
-
- for( unsigned short i = 0 ; i < PreTreatedData->GetECCEMult() ; i++ )
- {
- // cout << PreTreatedData->GetECCEClover(i) << endl;
- if( PreTreatedData->GetECCEClover(i) != DetectorID)
- {
- if(i==0)
- {
- NumberOfHitClover++;
- }
- else if(PreTreatedData->GetECCEClover(i)!= PreTreatedData->GetECCEClover(i-1) )
- {
- NumberOfHitClover++;
- }
+ if(m_EventData->GetExoOuter2(EventId) > OuterMax){
+ OuterMax = m_EventData->GetExoOuter2(EventId);
+ OuterId = 1;
}
- if(NumberOfHitClover == 4) break;
- //clover_mult -> Fill(NumberOfHitClover);
-
- }
-
- //cout << "NumberOfHitClover " << NumberOfHitClover << endl;
-
- map<int, vector<int> > MapCristal;
- map<int, vector<int> > MapSegment;
-
- map<int, vector<int> > :: iterator it; // iterator used with MapCristal
- map<int, vector<int> > :: iterator at; // iterator used with MapSegment
-
- vector<int> PositionOfCristal_Buffer_ECC;
- vector<int> PositionOfSegment_Buffer_GOCCE;
-
-
- //Fill map Cristal
- for(int clo = 0; clo < NumberOfClover; clo++)
- {
- for(unsigned int k = 0; k < ECC_CloverNumber.size(); k++)
- {
- if(ECC_CloverNumber.at(k) == clo) // && ECC_CristalNumber.at(k)== cri )
- PositionOfCristal_Buffer_ECC.push_back(k);
+ if(m_EventData->GetExoOuter3(EventId) > OuterMax){
+ OuterMax = m_EventData->GetExoOuter3(EventId);
+ OuterId = 2;
}
- if(PositionOfCristal_Buffer_ECC.size() != 0) MapCristal[clo] = PositionOfCristal_Buffer_ECC;
-
- PositionOfCristal_Buffer_ECC.clear();
-
- }
-
-
- //Fill map Segment
- for(int clo = 0; clo < NumberOfClover; clo++)
- {
- for(int cri = 0; cri < 4 ; cri++)
- {
- // for(int seg = 0; seg < 4 ; seg++)
- {
- for(unsigned int m = 0; m < GOCCE_CloverNumber.size(); m++)
- {
- if(GOCCE_CloverNumber.at(m) == clo && GOCCE_CristalNumber.at(m) == cri)// && GOCCE_SegmentNumber.at(m) == seg)
- {
- // PositionOfSegment_Buffer_GOCCE.push_back(4*clo+cri);
- PositionOfSegment_Buffer_GOCCE.push_back(m);
- }
- }
- }
- if(PositionOfSegment_Buffer_GOCCE.size() != 0) MapSegment[4*clo+cri] = PositionOfSegment_Buffer_GOCCE;
-
- PositionOfSegment_Buffer_GOCCE.clear();
+ if(m_EventData->GetExoOuter4(EventId) > OuterMax){
+ OuterMax = m_EventData->GetExoOuter4(EventId);
+ OuterId = 3;
}
- }
-
- // Treatment
- for(int clo = 0; clo < NumberOfClover ; clo++)
- {
- double E = 0; double T = 0;
- int mult_cristal = 0;
- int cristal = -1 , segment;
-
- int cristal_Emax = 0; int cristal_Emin = 0;
- int Emax = 0, Emin = 1000000;
- int Tmin = 0, Tmax = 0;
-
- //ADD-BACK
- it = MapCristal.find(clo);
-
- int cristal_cond = 0;
-
- if(it != MapCristal.end())
- {
- vector<int> PositionOfCristal = it -> second;
-
- mult_cristal = PositionOfCristal.size();
- //if(mult_cristal!=0) cristal_mult -> Fill(mult_cristal);
-
- // ADD-BACK
- //cout << "boucle" << endl;
-
- for(unsigned int k = 0; k < PositionOfCristal.size(); k++)
- {
- int indice = PositionOfCristal.at(k);
-
- cristal_cond += ECC_CristalNumber.at(indice);
- // cout << ECC_CristalNumber.at(k) << " " ECC_E.at(k) << endl;
-
- if(mult_cristal < 3)
- {
- E+= ECC_E.at(indice);
-
- if(ECC_E.at(indice) < Emin) {
- cristal_Emin = ECC_CristalNumber.at(indice);
- Emin = ECC_E.at(indice);
- Tmin = ECC_T.at(indice);
- }
-
- if(ECC_E.at(indice) > Emax) {
- cristal_Emax = ECC_CristalNumber.at(indice);
- Emax = ECC_E.at(indice);
- Tmax = ECC_T.at(indice);
- }
- }
-
- else // case of multiplicity = 3 or 4
- {
- E = -1; cristal_Emax = -1; cristal_Emin = -1; Tmax = -1; Tmin = -1;
- }
-
- // cout << ECC_E.at(indice) << " " << Emax << " " << cristal_Emax << " " << Emin << " " << cristal_Emin << endl;
-
- }
-
- if( (mult_cristal==1) || (mult_cristal ==2 && cristal_cond %2 == 1) )
- {
- // cout << cristal_cond << endl;
-
- //cristal = cristal_Emax; T = Tmax;
- //cout << Emax << " " << cristal_Emax << " " << Emin << " " << cristal_Emin << endl;
-
- if(E > 500) { cristal = cristal_Emax; T = Tmax; }
- else { cristal = cristal_Emin; T = Tmin; }
-
-
- // DOPPLER CORRECTION
-
- at = MapSegment.find(4*clo+cristal);
- segment = -1;
-
- if(at != MapSegment.end())
- {
- vector<int> PositionOfSegment = at -> second; // position of segment k in the vector
-
- int segment_max = -1, E_temp = -1;
-
- for(unsigned int m = 0; m < PositionOfSegment.size(); m++) // loop on hit segments of cristal cri of
- clover clo
- {
- int indice = PositionOfSegment.at(m);
-
- if(GOCCE_E.at(indice) > 0 && GOCCE_CristalNumber.at(indice) == cristal)
- {
- if( GOCCE_E.at(indice) > E_temp )
- {
- segment_max = GOCCE_SegmentNumber.at(indice) ;
- E_temp = GOCCE_E.at(indice);
- }
- }
- }
- segment = segment_max;
- }
-
- }
-
-
- if(E > 0 && cristal != -1 && segment != -1)
- {
- TotalEnergy_lab.push_back(E);
- Time.push_back(T);
- CloverNumber.push_back(clo);
- CristalNumber.push_back(cristal);
- SegmentNumber.push_back(segment);
-
- double theta = GetSegmentAngleTheta(clo, cristal, segment);
-
- Theta.push_back(theta);
-
- double doppler_E = DopplerCorrection(E, theta);
- DopplerCorrectedEnergy.push_back(doppler_E);
-
- // cout << E << " " << clo << " " << cristal << " " << segment << " " << theta << " " << doppler_E << endl;
-
- }
-
- } // end of condition over CristalMap
-
- } // loop over NumberOfClover
-
- CloverMult = GetClover_Mult();
+ return OuterId;
+}
- //cout << "Exogam fine" << endl;
- */
+double TExogamPhysics::GetDoppler(double Energy, unsigned int Flange, unsigned int Crystal, unsigned int Outer){
+ Exogam_struc = Ask_For_Angles(Flange, ComputeMeanFreePath(Energy));
+ double Theta_seg = Exogam_struc.Theta_Crystal_Seg[Crystal][Outer];
+ double Phi_seg = Exogam_struc.Phi_Crystal_Seg[Crystal][Outer];
+ return Doppler_Correction(Theta_seg,Phi_seg,0,0,Beta,Energy);
}
double TExogamPhysics::ComputeMeanFreePath(double Energy){
@@ -434,9 +277,21 @@ void TExogamPhysics::Clear() {
NumberOfHitCristal = 0;
E_Cal.clear();
- E_Doppler.clear();
- Flange_N.clear();
- Crystal_N.clear();
+ EHG_Cal.clear();
+ Outer1_Cal.clear();
+ Outer2_Cal.clear();
+ Outer3_Cal.clear();
+ Outer4_Cal.clear();
+ FlangeN.clear();
+ CrystalN.clear();
+ // E_Doppler.clear();
+ E_AB.clear();
+ FlangeN_AB.clear();
+ Size_AB.clear();
+ CrystalN_ABD.clear();
+ OuterN_ABD.clear();
+ E_ABD.clear();
+
//
// ECC_CloverNumber.clear();
// ECC_CristalNumber.clear();
@@ -887,7 +742,7 @@ void TExogamPhysics::DoCalibration() {
}
for (it = DoCalibrationE.begin(); it != DoCalibrationE.end(); it++) {
if (it->second) {
- DoCalibrationE_F(it->first,"E", calib_file, dispersion_file);
+ DoCalibrationE_F(it->first,"E", calib_file, dispersion_file, Threshold_E_Cal);
}
}
calib_file->close();
@@ -899,7 +754,7 @@ void TExogamPhysics::DoCalibration() {
}
for (it = DoCalibrationEHG.begin(); it != DoCalibrationEHG.end(); it++) {
if (it->second) {
- DoCalibrationE_F(it->first,"EHG", calib_file, dispersion_file);
+ DoCalibrationE_F(it->first,"EHG", calib_file, dispersion_file, Threshold_EHG_Cal);
}
}
calib_file->close();
@@ -912,7 +767,7 @@ void TExogamPhysics::DoCalibration() {
for (it2 = DoCalibrationOuter.begin(); it2 != DoCalibrationOuter.end(); it2++) {
for (it = (it2->second).begin(); it != (it2->second).end(); it++) {
if (it->second) {
- DoCalibrationE_F(it->first,Form("Outer%d_",it2->first), calib_file, dispersion_file);
+ DoCalibrationE_F(it->first,Form("Outer%d_",it2->first), calib_file, dispersion_file, Threshold_Outers_Cal);
}
}
}
@@ -935,7 +790,7 @@ void TExogamPhysics::MakeOuterCalibFolders(std::string make_folder) {
int sys =system((make_folder+"/Exogam_Outer").c_str());
}
-void TExogamPhysics::DoCalibrationE_F(unsigned int DetectorNumber,std::string CalibType, ofstream* calib_file, ofstream* dispersion_file) {
+void TExogamPhysics::DoCalibrationE_F(unsigned int DetectorNumber,std::string CalibType, ofstream* calib_file, ofstream* dispersion_file, unsigned int Threshold) {
auto TH1Map = RootHistogramsCalib::getInstance()->GetTH1Map();
auto TGraphMap = RootHistogramsCalib::getInstance()->GetTGraphMap();
@@ -961,7 +816,7 @@ void TExogamPhysics::DoCalibrationE_F(unsigned int DetectorNumber,std::string Ca
CubixEnergyCal->UseFirstDerivativeSearch();
- CubixEnergyCal->SetGlobalChannelLimits(hist->GetXaxis()->GetBinLowEdge(1),hist->GetXaxis()->GetBinLowEdge(hist->GetXaxis()->GetNbins())); // limit the search to this range in channels
+ CubixEnergyCal->SetGlobalChannelLimits(hist->GetXaxis()->GetBinLowEdge(1)+Threshold,hist->GetXaxis()->GetBinLowEdge(hist->GetXaxis()->GetNbins())); // limit the search to this range in channels
CubixEnergyCal->SetGlobalPeaksLimits(15,5); // default fwhm and minmum amplitude for the peaksearch [15 5]
CubixEnergyCal->StartCalib();
@@ -1160,7 +1015,7 @@ void TExogamPhysics::CreateCalibrationOuterFiles(ofstream* calib_file,
void TExogamPhysics::ReadDoCalibration(NPL::InputParser parser) {
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Exogam");
- vector<string> calibs = {"FirstCr","LastCr","FirstOuter","LastOuter","FitOrder","Source"};
+ vector<string> calibs = {"Threshold_E","Threshold_EHG","Threshold_Outers","FirstCr","LastCr","FirstOuter","LastOuter","FitOrder","Source"};
for (unsigned int i = 0; i < blocks.size(); i++) {
if (blocks[i]->HasTokenList(calibs)) {
@@ -1172,6 +1027,9 @@ void TExogamPhysics::ReadDoCalibration(NPL::InputParser parser) {
unsigned int LastOuter = blocks[i]->GetInt("LastOuter");
FitPolOrder = blocks[i]->GetInt("FitOrder");
Source_name = blocks[i]->GetString("Source");
+ Threshold_E_Cal = blocks[i]->GetInt("Threshold_E");
+ Threshold_EHG_Cal = blocks[i]->GetInt("Threshold_EHG");
+ Threshold_Outers_Cal = blocks[i]->GetInt("Threshold_Outers");
for(unsigned int k = FirstCr; k <= LastCr; k++){
DoCalibrationE[k] = true;
DoCalibrationEHG[k] = true;
@@ -1272,14 +1130,14 @@ namespace EXOGAM_LOCAL {
static string name;
name = "EXO/E";
name += NPL::itoa(m_EventData->GetExoCrystal(i));
- return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoE(i));
+ return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoE(i),1);
}
double fEXO_EHG(const TExogamData* m_EventData, const unsigned int& i) {
static string name;
name = "EXO/EHG";
name += NPL::itoa(m_EventData->GetExoCrystal(i));
- return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoEHG(i));
+ return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoEHG(i),1);
}
double fEXO_T(const TExogamData* m_EventData, const unsigned int& i) {
@@ -1287,7 +1145,7 @@ namespace EXOGAM_LOCAL {
name = "EXOGAM/Cr_";
name += NPL::itoa(m_EventData->GetExoCrystal(i));
name += "_TDC";
- return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoTDC(i));
+ return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoTDC(i),1);
}
double fEXO_Outer(const TExogamData* m_EventData, const unsigned int& i, const unsigned int OuterNumber) {
@@ -1297,13 +1155,13 @@ namespace EXOGAM_LOCAL {
name += "_";
name += NPL::itoa(OuterNumber);
if(OuterNumber == 0)
- return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoOuter1(i));
+ return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoOuter1(i),1);
else if(OuterNumber == 1)
- return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoOuter2(i));
+ return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoOuter2(i),1);
else if(OuterNumber == 2)
- return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoOuter3(i));
+ return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoOuter3(i),1);
else if(OuterNumber == 3)
- return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoOuter4(i));
+ return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoOuter4(i),1);
else{
std::cout << "WARNING: Outer number != 0-3, something is wrong\n";
return 0;
diff --git a/NPLib/Detectors/Exogam/TExogamPhysics.h b/NPLib/Detectors/Exogam/TExogamPhysics.h
index 149fcda3ddac803b7ee4a458b316dfdcfb517b03..a923e2cf4be32d735372236b5acdb588debe02dc 100644
--- a/NPLib/Detectors/Exogam/TExogamPhysics.h
+++ b/NPLib/Detectors/Exogam/TExogamPhysics.h
@@ -62,10 +62,24 @@ class TExogamPhysics : public TObject, public NPL::VDetector, public TExogamPhys
void Clear() ;
void Clear(const Option_t*) {};
+ // Only threshold and cal applied to exogam
std::vector<double> E_Cal;
- std::vector<double> E_Doppler;
- std::vector<unsigned int> Flange_N;
- std::vector<unsigned int> Crystal_N;
+ std::vector<double> EHG_Cal;
+ std::vector<double> Outer1_Cal;
+ std::vector<double> Outer2_Cal;
+ std::vector<double> Outer3_Cal;
+ std::vector<double> Outer4_Cal;
+ // std::vector<double> E_Doppler;
+ std::vector<unsigned int> FlangeN;
+ std::vector<unsigned int> CrystalN;
+
+ // Previous treatment + Add_Back (size of vectors are not the same because of AB !)
+ std::vector<double> E_AB;
+ std::vector<unsigned int> FlangeN_AB;
+ std::vector<unsigned int> Size_AB;
+ std::vector<unsigned int> CrystalN_ABD;
+ std::vector<int> OuterN_ABD;
+ std::vector<double> E_ABD;
@@ -180,7 +194,7 @@ class TExogamPhysics : public TObject, public NPL::VDetector, public TExogamPhys
void FillRootHistogramsCalib_F(); //!
- void DoCalibrationE_F(unsigned int Detector_Nbr,std::string CalibType, ofstream* calib_file, ofstream* dispersion_file); //!
+ void DoCalibrationE_F(unsigned int Detector_Nbr,std::string CalibType, ofstream* calib_file, ofstream* dispersion_file, unsigned int Threshold); //!
void DoCalibrationEHG_F(unsigned int Detector_Nbr, ofstream* calib_file, ofstream* dispersion_file){}; //!
@@ -244,9 +258,14 @@ class TExogamPhysics : public TObject, public NPL::VDetector, public TExogamPhys
unsigned int GetFlangeNbr(unsigned int crystal_nbr); //!
+ int GetMaxOuter(unsigned int EventId); //!
+
+ double GetDoppler(double Energy, unsigned int Flange, unsigned int Crystal, unsigned int Outer); //!
+
private: // Variables for analysis
unsigned int m_EXO_Mult;//!
+ double m_EXO_OuterUp_RAW_Threshold;//!
double m_EXO_E_RAW_Threshold;//!
double m_EXO_E_Threshold;//!
double m_EXO_EHG_RAW_Threshold;//!
@@ -260,21 +279,23 @@ class TExogamPhysics : public TObject, public NPL::VDetector, public TExogamPhys
double EXO_Outer4;//!
unsigned int flange_nbr;//!
unsigned int crystal_nbr;//!
- double EnergyDoppler; //!
double Beta = 0.257;//!
- double mean_free_path;//!
Clover_struc Exogam_struc;//!
std::map<unsigned int,std::pair<unsigned int,unsigned int>> MapCrystalFlangeCLover;//! Map key is raw crystal nbr, pair associated is flange nbr and crystal nbr in the flange
- const double GeDensity = 0.005323; //! g/mm3
+ double GeDensity = 0.005323; //! g/mm3
std::map<double, double> Map_PhotonCS;//!
map<int, bool> DoCalibrationE; //!
map<int, bool> DoCalibrationEHG; //!
map<int, bool> DoCalibrationT; //!
map<int, map<int,bool>> DoCalibrationOuter; //!
+
+ unsigned int Threshold_E_Cal;//!
+ unsigned int Threshold_EHG_Cal;//!
+ unsigned int Threshold_Outers_Cal;//!
std::vector<std::string> Source_isotope; //!
std::vector<double> Source_E; //!
diff --git a/NPLib/Detectors/Exogam/TExogamPhysicsReader.cxx b/NPLib/Detectors/Exogam/TExogamPhysicsReader.cxx
index 480a5292323f60b8696be46863aa966ca56e837c..5d7bec39e0491511a5f18c28ecdae205c358a697 100644
--- a/NPLib/Detectors/Exogam/TExogamPhysicsReader.cxx
+++ b/NPLib/Detectors/Exogam/TExogamPhysicsReader.cxx
@@ -29,5 +29,6 @@ TExogamPhysicsReader::TExogamPhysicsReader()
};
void TExogamPhysicsReader::r_SetTreeReader(TTreeReader* TreeReader){
-r_ReaderEventData = new TTreeReaderValue<TExogamData>(*TreeReader,"Exogam");
+ std::cout << "TEST " << TreeReader << std::endl;
+ r_ReaderEventData = new TTreeReaderValue<TExogamData>(*TreeReader,"Exogam");
};
\ No newline at end of file
diff --git a/NPLib/Detectors/MUST2/TMust2Physics.cxx b/NPLib/Detectors/MUST2/TMust2Physics.cxx
index b82b7d03f89db45b4267076896e023fe54da6fa8..e03492f3b77aadbd32a8b2e6646b802599b3c6b1 100644
--- a/NPLib/Detectors/MUST2/TMust2Physics.cxx
+++ b/NPLib/Detectors/MUST2/TMust2Physics.cxx
@@ -1232,11 +1232,12 @@ void TMust2Physics::AddParameterToCalibrationManager() {
vector<double> standardCsI;
for (int i = 0; i < m_NumberOfTelescope; ++i) {
- if (m_CsIOffset[i] == 1) {
+ if (m_CsIOffset[i+1] == 1) {
standardCsI = {0, 500. / 16384.};
}
- else
+ else{
standardCsI = {-250, 250. / 8192.};
+ }
for (int j = 0; j < 128; ++j) {
Cal->AddParameter("MUST2", "T" + NPL::itoa(i + 1) + "_Si_X" + NPL::itoa(j + 1) + "_E",
diff --git a/NPLib/Detectors/ZDD/TZDDPhysics.cxx b/NPLib/Detectors/ZDD/TZDDPhysics.cxx
index 8ac9064dd13e62921f694b008a7351395887048c..ff4b080edf1a754eb6bd6f3fb6d0319252d00e56 100644
--- a/NPLib/Detectors/ZDD/TZDDPhysics.cxx
+++ b/NPLib/Detectors/ZDD/TZDDPhysics.cxx
@@ -110,6 +110,7 @@ void TZDDPhysics::Match_IC(){
ICSum = 0;
for(auto it = SortIC.begin(); it != SortIC.end(); ++it){
ICSum+= (it->second).first;
+ // std::cout << "Test " << it->first << " " << it->second.first << std::endl;
IC_Nbr.push_back(it->first);
IC_E.push_back((it->second).first);
IC_TS.push_back((it->second).second);
@@ -145,6 +146,13 @@ void TZDDPhysics::Match_PL(){
PL_TS.push_back((it->second).second);
}
}
+
+// bool TZDDPhysics::CheckGoodEvent(){
+// if (NPOptionManager::getInstance()->IsReader() == true) {
+// m_EventData = &(**r_ReaderEventData);
+// }
+// return abs(ZDD->)
+// }
///////////////////////////////////////////////////////////////////////////
void TZDDPhysics::PreTreat() {
// This method typically applies thresholds and calibrations
diff --git a/NPLib/Detectors/ZDD/TZDDPhysics.h b/NPLib/Detectors/ZDD/TZDDPhysics.h
index 2847c289c42dfe321e948e6022df2f737520f67d..6ebfb408e37acbd50e37e7152470ad5c5e5dac03 100644
--- a/NPLib/Detectors/ZDD/TZDDPhysics.h
+++ b/NPLib/Detectors/ZDD/TZDDPhysics.h
@@ -150,6 +150,7 @@ class TZDDPhysics : public TObject, public NPL::VDetector, public TZDDPhysicsRea
void Match_IC1();
void Match_PL();
+ // bool CheckGoodEvent();
// PreTreating Energy for IC and Plastic
void PreTreatEnergy(std::string Detector, CalibrationManager* Cal);
diff --git a/NPLib/Detectors/ZDD/TZDDPhysicsReader.cxx b/NPLib/Detectors/ZDD/TZDDPhysicsReader.cxx
index ec5369f60e6b487357b6be571c723071abf1448b..49e5c22b3d977fd2d7211e302fdaa1a648055a1d 100644
--- a/NPLib/Detectors/ZDD/TZDDPhysicsReader.cxx
+++ b/NPLib/Detectors/ZDD/TZDDPhysicsReader.cxx
@@ -29,5 +29,6 @@ TZDDPhysicsReader::TZDDPhysicsReader()
};
void TZDDPhysicsReader::r_SetTreeReader(TTreeReader* TreeReader){
+ std::cout << "TEST " << TreeReader << std::endl;
r_ReaderEventData = new TTreeReaderValue<TZDDData>(*TreeReader,"ZDD");
};
\ No newline at end of file
diff --git a/NPLib/Utility/npreader.cxx b/NPLib/Utility/npreader.cxx
index ed1402777f2cf4226edd1cb78bffd2a2e8a6a4ba..bb17fc88bd2b7de2a0d446471bb3c63a1d5d6a3e 100644
--- a/NPLib/Utility/npreader.cxx
+++ b/NPLib/Utility/npreader.cxx
@@ -93,7 +93,7 @@ int main(int argc , char** argv){
UserAnalysis->SetDetectorManager(myDetector);
std::cout << "Salut ////////////////////" << std::endl;
UserAnalysis->Init();
- // UserAnalysis->InitTreeReader(inputTreeReader);
+ UserAnalysis->InitTreeReader(inputTreeReader);
}
else{
std::string str_error=error;
@@ -135,6 +135,7 @@ int main(int argc , char** argv){
unsigned long new_nentries = 0 ;
int current_tree = 0 ;
int total_tree = Chain->GetNtrees();
+ int entry_max = NPOptionManager::getInstance()->GetNumberOfEntryToAnalyse();
bool IsPhysics = myOptionManager->GetInputPhysicalTreeOption();
@@ -184,17 +185,18 @@ int main(int argc , char** argv){
else{
if(!IsPhysics){
- while(inputTreeReader->Next()){
+ while(inputTreeReader->Next()&& treated < entry_max){
// Build the current event
if(UserAnalysis->UnallocateBeforeBuild()){
myDetector->BuildPhysicalEvent();
// User Analysis;
if(UserAnalysis->UnallocateBeforeTreat()){
- UserAnalysis->TreatEvent();
+ UserAnalysis->TreatEvent();
+
//Fill the tree
if(UserAnalysis->FillOutputCondition())
- myDetector->FillOutputTree();
+ myDetector->FillOutputTree();
}
}
current_tree = Chain->GetTreeNumber()+1;
@@ -202,7 +204,7 @@ int main(int argc , char** argv){
}
}
- else{
+ /*else{
std::cout << "test npa branch 3" << std::endl;
while(inputTreeReader->Next()){
//for (unsigned long i = first_entry ; i < nentries + first_entry; i++) {
@@ -238,6 +240,7 @@ int main(int argc , char** argv){
//}
}
}
+*/
UserAnalysis->End();
}
diff --git a/NPSimulation/Detectors/MUST2/MUST2Array.cc b/NPSimulation/Detectors/MUST2/MUST2Array.cc
index 831edab0fd9535d954542613cfdf110a5642eb6a..ec47e723258b3c572ef3c83058cbdc7108e8f6a3 100644
--- a/NPSimulation/Detectors/MUST2/MUST2Array.cc
+++ b/NPSimulation/Detectors/MUST2/MUST2Array.cc
@@ -933,7 +933,7 @@ void MUST2Array::ReadSensitive(const G4Event*) {
vector<unsigned int> level = CsIScorer->GetLevel(i);
if (ECsI > ThresholdCsI) {
if (m_CsIOffset[level[0]] == 1)
- m_Event->SetCsIE(level[0], level[1], NPL::EnergyToADC(ECsI, 0, 500, 0, 16384));
+ m_Event->SetCsIE(level[0], level[1], NPL::EnergyToADC(ECsI, 0, 0, 0, 16384));
else
m_Event->SetCsIE(level[0], level[1], NPL::EnergyToADC(ECsI, 0, 250, 8192, 16384));
double timeCsI = RandGauss::shoot(CsIScorer->GetTime(i), ResoTimeMust);
diff --git a/Projects/E805/Analysis.cxx b/Projects/E805/Analysis.cxx
index 22906dc7054dea28e5e6356e9298f4c949dc3475..ec0431e654fe27bc0eddeb0c358028c8b5c5f148 100755
--- a/Projects/E805/Analysis.cxx
+++ b/Projects/E805/Analysis.cxx
@@ -35,204 +35,195 @@ Analysis::~Analysis(){
void Analysis::Init(){
InitInputBranch();
InitOutputBranch();
- // CATS = (TCATSPhysics*) m_DetectorManager -> GetDetector("CATSDetector");
+ CATS = (TCATSPhysics*) m_DetectorManager -> GetDetector("CATSDetector");
M2 = (TMust2Physics*) m_DetectorManager -> GetDetector("M2Telescope");
- //
- //reaction->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
- //OriginalBeamEnergy = reaction->GetBeamEnergy();
-
-
- //string Path = "../../Inputs/EnergyLoss/";
- //string TargetMaterial = m_DetectorManager->GetTargetMaterial();
- //TargetThickness = m_DetectorManager->GetTargetThickness();
- //string beam= NPL::ChangeNameToG4Standard(reaction->GetNucleus1()->GetName());
- //string heavy_ejectile= NPL::ChangeNameToG4Standard(reaction->GetNucleus4()->GetName());
- //string light=NPL::ChangeNameToG4Standard(reaction->GetNucleus3()->GetName());
-
- //string Reaction_pd_s = "48Cr(p,d)47Cr@1620";
- //string Reaction_pt_s = "48Cr(p,t)46Cr@1620";
- //string Reaction_p3He_s = "48Cr(p,3He)46V@1620";
- //Reaction_pd = new Reaction(Reaction_pd_s);
- //Reaction_pt = new Reaction(Reaction_pt_s);
- //Reaction_p3He = new Reaction(Reaction_p3He_s);
-
-////
- //ProtonSi = NPL::EnergyLoss(Path+ "proton_Si.G4table", "G4Table", 100);
- //
- //for(unsigned int i = 0; i < ParticleType.size(); i++){
- // LightAl[ParticleType[i]] = NPL::EnergyLoss(Path+ParticleType[i]+"_Al.G4table","G4Table",100);
- // LightTarget[ParticleType[i]] = NPL::EnergyLoss(Path+ParticleType[i]+"_CH2.G4table","G4Table",100);
- //}
- //BeamTarget["48Cr"] = NPL::EnergyLoss(Path+"Cr48_CH2.G4table","G4Table",100);
+ ZDD = (TZDDPhysics*) m_DetectorManager -> GetDetector("ZDD");
+ TAC = (TTACPhysics*) m_DetectorManager -> GetDetector("TAC");
+ Exogam = (TExogamPhysics*) m_DetectorManager -> GetDetector("Exogam");
+ reaction->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
+ OriginalBeamEnergy = reaction->GetBeamEnergy();
+
+
+ string Path = "../../Inputs/EnergyLoss/";
+ string TargetMaterial = m_DetectorManager->GetTargetMaterial();
+ TargetThickness = m_DetectorManager->GetTargetThickness();
+ string beam= NPL::ChangeNameToG4Standard(reaction->GetNucleus1()->GetName());
+ string heavy_ejectile= NPL::ChangeNameToG4Standard(reaction->GetNucleus4()->GetName());
+ string light=NPL::ChangeNameToG4Standard(reaction->GetNucleus3()->GetName());
+
+ string Reaction_pd_s = "48Cr(p,d)47Cr@1620";
+ string Reaction_pt_s = "48Cr(p,t)46Cr@1620";
+ string Reaction_p3He_s = "48Cr(p,3He)46V@1620";
+ Reaction_pd = new Reaction(Reaction_pd_s);
+ Reaction_pt = new Reaction(Reaction_pt_s);
+ Reaction_p3He = new Reaction(Reaction_p3He_s);
+
+//
+ ProtonSi = NPL::EnergyLoss(Path+ "proton_Si.G4table", "G4Table", 100);
+
+ for(unsigned int i = 0; i < ParticleType.size(); i++){
+ LightAl[ParticleType[i]] = NPL::EnergyLoss(Path+ParticleType[i]+"_Al.G4table","G4Table",100);
+ LightTarget[ParticleType[i]] = NPL::EnergyLoss(Path+ParticleType[i]+"_CH2.G4table","G4Table",100);
+ }
+ BeamTarget["48Cr"] = NPL::EnergyLoss(Path+"Cr48_CH2.G4table","G4Table",100);
- //Reaction_pd->SetBeamEnergy(BeamTarget["48Cr"].Slow(Reaction_pd->GetBeamEnergy(),TargetThickness*0.5,0));
- //Reaction_pt->SetBeamEnergy(BeamTarget["48Cr"].Slow(Reaction_pt->GetBeamEnergy(),TargetThickness*0.5,0));
- //Reaction_p3He->SetBeamEnergy(BeamTarget["48Cr"].Slow(Reaction_p3He->GetBeamEnergy(),TargetThickness*0.5,0));
- //Cal = CalibrationManager::getInstance();
+ Reaction_pd->SetBeamEnergy(BeamTarget["48Cr"].Slow(Reaction_pd->GetBeamEnergy(),TargetThickness*0.5,0));
+ Reaction_pt->SetBeamEnergy(BeamTarget["48Cr"].Slow(Reaction_pt->GetBeamEnergy(),TargetThickness*0.5,0));
+ Reaction_p3He->SetBeamEnergy(BeamTarget["48Cr"].Slow(Reaction_p3He->GetBeamEnergy(),TargetThickness*0.5,0));
+ Cal = CalibrationManager::getInstance();
}
///////////////////////////// Initialize some important parameters //////////////////////////////////
bool Analysis::UnallocateBeforeBuild(){
- // return true;
+ GATCONFMASTER.clear();
+ GATCONFMASTERTS.clear();
+
GATCONFMASTER = **GATCONFMASTER_;
return (GATCONFMASTER.size() == 1 && GATCONFMASTER[0] > 0);
- //return true;
+ // return true;
}
bool Analysis::UnallocateBeforeTreat(){
-/* for(int i = 0; i < 10; i++){
- PlasticRaw[i] = (*PlasticRaw_ )[i];
- PlasticRawTS[i] = (*PlasticRaw_TS_)[i];
- }
-
- for(int i = 0; i < 6; i++){
- IC_ZDDRaw[i] = (*IC_ZDDRaw_)[i];
- IC_ZDDRawTS[i] = (*IC_ZDDRaw_TS_)[i];
- }
-
- TAC_CATS_PL = **TAC_CATS_PL_;
- TAC_CATS_PLTS = **TAC_CATS_PL_TS_;
-
- TAC_CATS_HF = **TAC_CATS_HF_;
- TAC_CATS_HFTS = **TAC_CATS_HF_TS_;
-
- TAC_CATS_EXOGAM = **TAC_CATS_EXOGAM_;
- TAC_CATS_EXOGAMTS = **TAC_CATS_EXOGAM_TS_;
-
- TAC_MMG_CATS2 = **TAC_MMG_CATS2_;
- TAC_MMG_CATS2TS = **TAC_MMG_CATS2_TS_;
-
- TAC_MMG_CATS1 = **TAC_MMG_CATS1_;
- TAC_MMG_CATS1TS = **TAC_MMG_CATS1_TS_;
-
- TAC_MMG_EXOGAM = **TAC_MMG_EXOGAM_;
- TAC_MMG_EXOGAMTS = **TAC_MMG_EXOGAM_TS_;
-
- TAC_CATS1_CATS2 = **TAC_CATS1_CATS2_;
- TAC_CATS1_CATS2TS = **TAC_CATS1_CATS2_TS_;
-
- TAC_D4_CATS1 =**TAC_D4_CATS1_;
- TAC_D4_CATS1 =**TAC_D4_CATS1_TS_;
-
- TAC_PL_1 = **TAC_PL_1_;
- TAC_PL_1TS = **TAC_PL_1_TS_;
- TAC_PL_2 = **TAC_PL_2_;
- TAC_PL_2TS = **TAC_PL_2_TS_;
- TAC_PL_3 = **TAC_PL_3_;
- TAC_PL_3TS = **TAC_PL_3_TS_;
- TAC_PL_4 = **TAC_PL_4_;
- TAC_PL_4TS = **TAC_PL_4_TS_;
- TAC_PL_5 = **TAC_PL_5_;
- TAC_PL_5TS = **TAC_PL_5_TS_;
-*/
+ GATCONFMASTERTS = **GATCONFMASTERTS_;
return true;
}
bool Analysis::FillOutputCondition(){
- return true;
- //return (CATS->MapX.size() > 0 && CATS->MapY.size()>0);
+ return bCATS;
+ //return true;
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::TreatEvent(){
- // if(M2->CsI_E.size() > 0)
- // std::cout << "Analysis test " << M2->CsI_E[0] << " " << M2->CsI_E.size() << " " << "\n \n";
-/*
ReInit();
- // std::cout << CATS->PositionX.size() << std::endl;
+
//////////////////// MUST2 Part ////////////////////
- int M2_size = M2->Si_E.size();
- if(CATS->PositionX.size() == 2 && CATS->PositionY.size() == 2){
- for(unsigned int countMust2 = 0 ; countMust2 < M2_size ; countMust2++){
- M2_TelescopeM++;
- // MUST2
- int TelescopeNumber = M2->TelescopeNumber[countMust2];
-
- // Part 1 : Impact Angle
- ThetaM2Surface = 0;
- ThetaNormalTarget = 0;
-
- BeamImpact = TVector3(CATS->PositionOnTargetX,CATS->PositionOnTargetY,0);
- // std::cout << "Position On target : " << CATS->PositionOnTargetX << " " << CATS->PositionOnTargetY << std::endl;
+ //TreatCATS();
+ bCATS = true;
+ //if(bCATS){
+ // TreatZDD();
+ // TreatTAC();
+ // TreatMUST2();
+ // TreatEXO();
+ //}
+ /*//for(unsigned int countMust2 = 0 ; countMust2 < M2->Si_E.size() ; countMust2++){
+ //Si_E_M2 = M2->Si_E[countMust2];
+ //CsI_E_M2 = M2->CsI_E[countMust2];
+ //ThetaM2Surface = 0;
+ // if(Si_E_M2 > 0 && CsI_E_M2 > 8192){
+ // //double EfromdeltaE = ProtonSi.EvaluateEnergyFromDeltaE(
+ // // Si_E_M2, 300*um, ThetaM2Surface, 6.0 * MeV, 300.0 * MeV,
+ // // 0.001 * MeV, 10000);
+ // double EfromdeltaE = (CsI_E_M2-8192)*0.1;
+ // M2_ECsI_from_deltaE.push_back(EfromdeltaE);
+ // if(EfromdeltaE > 0){
+ // Beta_light = sqrt(1./(1.+1./(pow(EfromdeltaE/911. + 1,2)-1)));
+ // Beta_from_deltaE.push_back(Beta_light);
+ // if(Beta_light>0){
+ // double Beth = log(2*511.*Beta_light*Beta_light/(0.174*(1-Beta_light*Beta_light))) - Beta_light*Beta_light;
+ // Beth_from_deltaE.push_back(Beth);
+ // }
+ // }
+ // }
+ //}
+*/
+}
- //BeamImpact = TVector3(0,0,0);
-
- BeamDirection = TVector3(CATS->PositionX[0] - CATS->PositionX[1],CATS->PositionY[0] - CATS->PositionY[1],CATS->PositionZ[0] - CATS->PositionZ[1]);
- // std::cout << "Position XY " << CATS->PositionX[1] - CATS->PositionX[0] << " " << CATS->PositionY[1] - CATS->PositionY[0] << " " << CATS->PositionZ[1] - CATS->PositionZ[0] << std::endl;
- // BeamDirection = TVector3(0,0,1);
-
- TVector3 HitDirection = M2 -> GetPositionOfInteraction(countMust2) - BeamImpact ;
- M2_ThetaLab.push_back(HitDirection.Angle( BeamDirection ));
-
- //std::cout << M2 -> GetPositionOfInteraction(countMust2).X() << " " << M2 -> GetPositionOfInteraction(countMust2).Y() << " "<< M2 -> GetPositionOfInteraction(countMust2).Z() << std::endl;
- M2_X.push_back(M2 -> GetPositionOfInteraction(countMust2).X());
- M2_Y.push_back(M2 -> GetPositionOfInteraction(countMust2).Y());
- M2_Z.push_back(M2 -> GetPositionOfInteraction(countMust2).Z());
-
- ThetaM2Surface = HitDirection.Angle(- M2 -> GetTelescopeNormal(countMust2) );
- ThetaNormalTarget = HitDirection.Angle( TVector3(0,0,1) ) ;
-
- // Part 2 : Impact Energy
- int CristalNb = M2->CsI_N[countMust2];
-
- Si_E_M2 = M2->Si_E[countMust2];
- CsI_E_M2= M2->CsI_E[countMust2];
- // if(CsI_E_M2 > 0)
- // std::cout << "Analysis " << CsI_E_M2 << " " << M2->CsI_E.size() << " " << CristalNb << "\n \n";
- for(unsigned int i = 0; i < ParticleType.size(); i++){
- Energy[ParticleType[i]] = 0;
- CsI_Energy[ParticleType[i]] = 0;
-
- if(CsI_E_M2>8192 ){
- // The energy in CsI is calculate form dE/dx Table because
- CsI_Energy[ParticleType[i]] = Cal->ApplyCalibration("MUST2/"+ParticleType[i]+"_T"+NPL::itoa(TelescopeNumber)+"_CsI"+NPL::itoa(CristalNb)+"_E",CsI_E_M2);
- Energy[ParticleType[i]] = CsI_Energy[ParticleType[i]];
- // Energy[ParticleType[i]] = LightAl[ParticleType[i]].EvaluateInitialEnergy(Energy[ParticleType[i]], 0.4 * micrometer, ThetaM2Surface);
- //std::cout << ParticleType[i]+"MUST2/T"+NPL::itoa(TelescopeNumber)+"_CsI"+NPL::itoa(CristalNb)+"_E" << " " << Energy[ParticleType[i]] << "\n";
- //Energy = LightAl.EvaluateInitialEnergy( Energy ,0.4*micrometer , ThetaM2Surface);
- //Energy+=Si_E_M2;
- }
+void Analysis::TreatCATS(){
+ BeamImpact = TVector3(CATS->PositionOnTargetX,CATS->PositionOnTargetY,0);
+ // BeamImpact = TVector3(0,0,0);
+ // std::cout << "Position On target : " << CATS->PositionOnTargetX << " " << CATS->PositionOnTargetY << std::endl;
+ BeamDirection = TVector3(CATS->PositionX[0] - CATS->PositionX[1],CATS->PositionY[0] - CATS->PositionY[1],CATS->PositionZ[0] - CATS->PositionZ[1]);
+ // BeamDirection = TVector3(0,0,1);
+ // std::cout << "Position XY " << CATS->PositionX[1] - CATS->PositionX[0] << " " << CATS->PositionY[1] - CATS->PositionY[0] << " " << CATS->PositionZ[1] - CATS->PositionZ[0] << std::endl;
+ bCATS = true;
+}
-// else
- Energy[ParticleType[i]] += Si_E_M2;
- if(Energy[ParticleType[i]] > 0)
- Energy[ParticleType[i]] = LightTarget[ParticleType[i]].EvaluateInitialEnergy(Energy[ParticleType[i]] ,TargetThickness*0.5, ThetaNormalTarget);
- else
- Energy[ParticleType[i]] = -1000;
+void Analysis::TreatZDD(){
+
+}
+void Analysis::TreatTAC(){
+
+}
- }
- // Evaluate energy using the thickness
- //Energy = LightAl.EvaluateInitialEnergy( Energy ,0.4*micrometer , ThetaM2Surface);
- // Target Correction
- //M2_ELab.push_back(LightTarget.EvaluateInitialEnergy(Energy ,TargetThickness*0.5, ThetaNormalTarget));
- M2_ELab.push_back(Energy["proton"]);
- M2_dE.push_back(Si_E_M2);
-
- // Part 3 : Excitation Energy Calculation
- M2_Ex_p.push_back(reaction->ReconstructRelativistic( Energy["proton"] , M2_ThetaLab[countMust2] ));
- M2_Ex_d.push_back(Reaction_pd->ReconstructRelativistic( Energy["deuteron"] , M2_ThetaLab[countMust2] ));
- std::cout << "oui " << M2_Ex_d[countMust2] << std::endl;
- M2_Ex_t.push_back(Reaction_pt->ReconstructRelativistic( Energy["triton"] , M2_ThetaLab[countMust2] ));
- M2_Ex_a.push_back(Reaction_p3He->ReconstructRelativistic( Energy["alpha"] , M2_ThetaLab[countMust2] ));
-
- M2_CsI_E_p.push_back(CsI_Energy["proton"]);
- M2_CsI_E_d.push_back(CsI_Energy["deuteron"]);
- M2_CsI_E_t.push_back(CsI_Energy["triton"]);
- M2_CsI_E_a.push_back(CsI_Energy["alpha"]);
+void Analysis::TreatMUST2(){
+
+ int M2_size = M2->Si_E.size();
+ for(unsigned int countMust2 = 0 ; countMust2 < M2_size ; countMust2++){
+ M2_TelescopeM++;
+ int TelescopeNumber = M2->TelescopeNumber[countMust2];
+
+ // Part 1 : Impact Angle
+ ThetaM2Surface = 0;
+ ThetaNormalTarget = 0;
- M2_ThetaLab[countMust2]=M2_ThetaLab[countMust2]/deg;
+ TVector3 HitDirection = M2 -> GetPositionOfInteraction(countMust2) - BeamImpact ;
+ M2_ThetaLab.push_back(HitDirection.Angle( BeamDirection ));
+
+ M2_X.push_back(M2 -> GetPositionOfInteraction(countMust2).X());
+ M2_Y.push_back(M2 -> GetPositionOfInteraction(countMust2).Y());
+ M2_Z.push_back(M2 -> GetPositionOfInteraction(countMust2).Z());
+
+ ThetaM2Surface = HitDirection.Angle(- M2 -> GetTelescopeNormal(countMust2) );
+ ThetaNormalTarget = HitDirection.Angle( TVector3(0,0,1) ) ;
+
+ // Part 2 : Impact Energy
+ int CristalNb = M2->CsI_N[countMust2];
+
+ Si_E_M2 = M2->Si_E[countMust2];
+ CsI_E_M2= M2->CsI_E[countMust2];
+
+ for(unsigned int i = 0; i < ParticleType.size(); i++){
+ Energy[ParticleType[i]] = 0;
+ CsI_Energy[ParticleType[i]] = 0;
+
+ if(CsI_E_M2>8192 ){
+ // The energy in CsI is calculate form dE/dx Table because
+ std::string name = "MUST2/"+ParticleType[i]+"_T"+NPL::itoa(TelescopeNumber)+"_CsI"+NPL::itoa(CristalNb)+"_E";
+ CsI_Energy[ParticleType[i]] = Cal->ApplyCalibration(name,CsI_E_M2);
+ Energy[ParticleType[i]] = CsI_Energy[ParticleType[i]];
+ }
- // Part 4 : Theta CM Calculation
- M2_ThetaCM.push_back(reaction->EnergyLabToThetaCM( M2_ELab[countMust2] , M2_ThetaLab[countMust2])/deg);
+ Energy[ParticleType[i]] += Si_E_M2;
+
+ if(Energy[ParticleType[i]] > 0)
+ Energy[ParticleType[i]] = LightTarget[ParticleType[i]].EvaluateInitialEnergy(Energy[ParticleType[i]] ,TargetThickness*0.5, ThetaNormalTarget);
+ else
+ Energy[ParticleType[i]] = -1000;
+
+
+ }
+ // Target Correction
+ //M2_ELab.push_back(LightTarget.EvaluateInitialEnergy(Energy ,TargetThickness*0.5, ThetaNormalTarget));
+ M2_ELab.push_back(Energy["proton"]);
+ M2_dE.push_back(Si_E_M2);
+
+ // Part 3 : Excitation Energy Calculation
+ M2_Ex_p.push_back(reaction->ReconstructRelativistic( Energy["proton"] , M2_ThetaLab[countMust2] ));
+ M2_Ex_d.push_back(Reaction_pd->ReconstructRelativistic( Energy["deuteron"] , M2_ThetaLab[countMust2] ));
+ M2_Ex_t.push_back(Reaction_pt->ReconstructRelativistic( Energy["triton"] , M2_ThetaLab[countMust2] ));
+ M2_Ex_a.push_back(Reaction_p3He->ReconstructRelativistic( Energy["alpha"] , M2_ThetaLab[countMust2] ));
+
+ M2_CsI_E_p.push_back(CsI_Energy["proton"]);
+ M2_CsI_E_d.push_back(CsI_Energy["deuteron"]);
+ M2_CsI_E_t.push_back(CsI_Energy["triton"]);
+ M2_CsI_E_a.push_back(CsI_Energy["alpha"]);
+
+ M2_ThetaLab[countMust2]=M2_ThetaLab[countMust2]/deg;
+
+ // Part 4 : Theta CM Calculation
+ M2_ThetaCM.push_back(reaction->EnergyLabToThetaCM( M2_ELab[countMust2] , M2_ThetaLab[countMust2])/deg);
- //if(proton_cut[TelescopeNumber]){
- // std::cout << "Test :" << proton_cut[TelescopeNumber] << " \n";
- //
- }//end loop MUST2
}
+
+}
+
+void Analysis::TreatEXO(){
+ // }
+ /*
if(M2->Si_E.size() > 0){
if(Inner6MVM > 0){
@@ -272,36 +263,15 @@ void Analysis::TreatEvent(){
}
}
-
- //for(unsigned int countMust2 = 0 ; countMust2 < M2->Si_E.size() ; countMust2++){
- //Si_E_M2 = M2->Si_E[countMust2];
- //CsI_E_M2 = M2->CsI_E[countMust2];
- //ThetaM2Surface = 0;
- // if(Si_E_M2 > 0 && CsI_E_M2 > 8192){
- // //double EfromdeltaE = ProtonSi.EvaluateEnergyFromDeltaE(
- // // Si_E_M2, 300*um, ThetaM2Surface, 6.0 * MeV, 300.0 * MeV,
- // // 0.001 * MeV, 10000);
- // double EfromdeltaE = (CsI_E_M2-8192)*0.1;
- // M2_ECsI_from_deltaE.push_back(EfromdeltaE);
- // if(EfromdeltaE > 0){
- // Beta_light = sqrt(1./(1.+1./(pow(EfromdeltaE/911. + 1,2)-1)));
- // Beta_from_deltaE.push_back(Beta_light);
- // if(Beta_light>0){
- // double Beth = log(2*511.*Beta_light*Beta_light/(0.174*(1-Beta_light*Beta_light))) - Beta_light*Beta_light;
- // Beth_from_deltaE.push_back(Beth);
- // }
- // }
- // }
- //}
-*/
+*/
}
-
void Analysis::InitOutputBranch() {
RootOutput::getInstance()->GetTree()->Branch("GATCONF",&GATCONFMASTER);
- /*
+ RootOutput::getInstance()->GetTree()->Branch("GATCONFTS",&GATCONFMASTERTS);
+
RootOutput::getInstance()->GetTree()->Branch("M2_TelescopeM",&M2_TelescopeM,"M2_TelescopeM/s");
RootOutput::getInstance()->GetTree()->Branch("M2_CsI_E_p",&M2_CsI_E_p);
RootOutput::getInstance()->GetTree()->Branch("M2_CsI_E_d",&M2_CsI_E_d);
@@ -322,7 +292,7 @@ void Analysis::InitOutputBranch() {
RootOutput::getInstance()->GetTree()->Branch("M2_dE",&M2_dE);
RootOutput::getInstance()->GetTree()->Branch("CsI_E_M2",&CsI_E_M2);
// RootOutput::getInstance()->GetTree()->Branch("M2_ECsI_from_deltaE",&M2_ECsI_from_deltaE);
-
+ /*
RootOutput:: getInstance()->GetTree()->Branch("TAC_CATS_PL",&TAC_CATS_PL,"TAC_CATS_PL/s");
RootOutput:: getInstance()->GetTree()->Branch("TAC_CATS_PLTS",&TAC_CATS_PLTS,"TAC_CATS_PLTS/l");
@@ -388,6 +358,7 @@ void Analysis::InitInputBranch(){
TTreeReader* inputTreeReader = RootInput::getInstance()->GetTreeReader();
GATCONFMASTER_ = new TTreeReaderValue<vector<unsigned int>>(*inputTreeReader,"GATCONF");
+ GATCONFMASTERTS_ = new TTreeReaderValue<vector<unsigned long long>>(*inputTreeReader,"GATCONFTS");
//DATATRIG_CATS_ = new TTreeReaderValue<unsigned short>(*inputTreeReader,"DATATRIG_CATS");
/*PlasticRaw_ = new TTreeReaderArray<UShort_t>(*inputTreeReader,"PlasticRaw");
PlasticRaw_TS_ = new TTreeReaderArray<ULong64_t>(*inputTreeReader,"PlasticRawTS");
diff --git a/Projects/E805/Analysis.h b/Projects/E805/Analysis.h
index 6ef5d390e57243032fdde2184587b4d1bdc1d248..e891c248947033f2fd507acb90e0de6806141d12 100755
--- a/Projects/E805/Analysis.h
+++ b/Projects/E805/Analysis.h
@@ -21,10 +21,12 @@
* *
*****************************************************************************/
-#include"Geometry_Clover_Exogam.h"
+// #include"Geometry_Clover_Exogam.h"
#include "TMust2Physics.h"
#include "TCATSPhysics.h"
-#include "TMust2PhysicsReader.h"
+#include "TTACPhysics.h"
+#include "TExogamPhysics.h"
+// #include "TMust2PhysicsReader.h"
#include"NPVAnalysis.h"
#include"TZDDPhysics.h"
#include"NPEnergyLoss.h"
@@ -68,7 +70,18 @@ class Analysis: public NPL::VAnalysis{
bool CheckExoDeltaTV(float ExoTime);
static NPL::VAnalysis* Construct();
- TTreeReaderValue<TMust2Data>* r_ReaderEventData;
+ TTreeReaderValue<TMust2Data>* r_ReaderEventData;
+
+ private:
+ void TreatCATS();
+ void TreatZDD();
+ void TreatTAC();
+ void TreatMUST2();
+ void TreatEXO();
+
+ private:
+
+ bool bCATS;
private:
@@ -149,6 +162,8 @@ class Analysis: public NPL::VAnalysis{
vector<unsigned int> GATCONFMASTER;
TTreeReaderValue<vector<unsigned int>>* GATCONFMASTER_;
+ vector<unsigned long long> GATCONFMASTERTS;
+ TTreeReaderValue<vector<unsigned long long>>* GATCONFMASTERTS_;
unsigned short DATATRIG_CATS;
TTreeReaderValue<unsigned short>* DATATRIG_CATS_;
@@ -333,6 +348,9 @@ class Analysis: public NPL::VAnalysis{
TMust2Physics* M2;
TCATSPhysics* CATS;
TZDDPhysics* ZDD;
+ TTACPhysics* TAC;
+ TExogamPhysics* Exogam;
+
TMust2PhysicsReader* M2_Reader;
TInitialConditions* InitialConditions;
TReactionConditions* ReactionConditions;
@@ -374,7 +392,7 @@ class Analysis: public NPL::VAnalysis{
double Theta_seg;
double Phi_seg;
- Char_t Exogam[100];
+ // Char_t Exogam[100];
Int_t ExoNumb;
Int_t Flange_tmp;
int FlangeNumb[12];
@@ -382,7 +400,7 @@ class Analysis: public NPL::VAnalysis{
/////////////Exogam
- Clover_struc Exogam_Clovers_struc[12];
+ // Clover_struc Exogam_Clovers_struc[12];
TCutG* proton_cut[4];
diff --git a/Projects/E805/Calibration1.txt b/Projects/E805/Calibration1.txt
index 165ab033dbbb98f9cfb8964149b1341602947787..9f2c6a87027dda7f53a467c2e30f37c5b2e93f28 100644
--- a/Projects/E805/Calibration1.txt
+++ b/Projects/E805/Calibration1.txt
@@ -18,6 +18,11 @@ CalibrationFilePath
./data/NPRoot/Calibration/Time/Coeff/Cal_Str_X_T_MM4_r0187.cal
./data/NPRoot/Calibration/Time/Coeff/Cal_Str_Y_T_MM4_r0188_manual.cal
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM1.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM2.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM3.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM4.txt
+
./data/NPRoot/Calibration/Run_CSICalib/CsI_MM1_proton.txt
./data/NPRoot/Calibration/Run_CSICalib/CsI_MM2_proton.txt
./data/NPRoot/Calibration/Run_CSICalib/CsI_MM3_proton.txt
@@ -39,3 +44,8 @@ CalibrationFilePath
./calibration/CATS/CATS1Y/CATSLin1Y_r0164.txt
./calibration/CATS/CATS2X/CATSLin2X_r0166.txt
./calibration/CATS/CATS2Y/CATSLin2Y_r0166.txt
+
+ ./data/NPRoot/Calibration/NPrun_392_152Eu/Exogam_E/Cal_Exogam_E.cal
+ ./data/NPRoot/Calibration/NPrun_392_152Eu/Exogam_EHG/Cal_Exogam_EHG.cal
+ ./data/NPRoot/Calibration/NPrun_392_152Eu/Exogam_Outer/Cal_Exogam_Outer.cal
+
\ No newline at end of file
diff --git a/Projects/E805/DetectorConfiguration/MUST2_E805.detector b/Projects/E805/DetectorConfiguration/MUST2_E805.detector
index 98663a5d77eea74774214e348eeb49c2912a101f..0609faabff7d887a6fe1924e5aaa9c28732a8259 100644
--- a/Projects/E805/DetectorConfiguration/MUST2_E805.detector
+++ b/Projects/E805/DetectorConfiguration/MUST2_E805.detector
@@ -1,6 +1,7 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Target
- THICKNESS= 53.5 micrometer
+ %THICKNESS= 53.5 micrometer
+ THICKNESS= 0.000001 micrometer
ANGLE= 0 deg
RADIUS= 15 mm
MATERIAL= CH2
diff --git a/Projects/E805/project.config b/Projects/E805/project.config
deleted file mode 100644
index f4d49abf59994412aa11998b16b6959b5abd9a3e..0000000000000000000000000000000000000000
--- a/Projects/E805/project.config
+++ /dev/null
@@ -1,8 +0,0 @@
-Project E805
- AnalysisOutput= ./data/NPRoot/Analysis/
- SimulationOutput= ./data/NPRoot/Simulation/
- EnergyLoss= ../../Inputs/EnergyLoss/
- CalibrationOutput= ./data/NPRoot/Calibration/
- Cuts= ./data/NPRoot/Cuts/
-
-
diff --git a/Projects/E805_simu/48Cr_pd_test.reaction b/Projects/E805_simu/48Cr_pd_test.reaction
index ffeab1b9ffc69ddbec7d91f1cb4009eb792023f4..99ae3a287f0f802851376d23328908055a7e3743 100644
--- a/Projects/E805_simu/48Cr_pd_test.reaction
+++ b/Projects/E805_simu/48Cr_pd_test.reaction
@@ -4,12 +4,22 @@
Beam
Particle= 48Cr
Energy= 1511 MeV
- SigmaEnergy= 42 MeV
+ %SigmaEnergy= 42 MeV
+ %ExcitationEnergy= 0 MeV
+ %SigmaThetaX= 0.1 mrad
+ %SigmaPhiY= 0.1 mrad
+ %SigmaX= 5 mm
+ %SigmaY= 5 mm
+ %MeanThetaX= 0 mrad
+ %MeanPhiY= 0 mrad
+ %MeanX= 0 mm
+ %MeanY= 0 mm
+ SigmaEnergy= 0 MeV
ExcitationEnergy= 0 MeV
- SigmaThetaX= 0.1 mrad
- SigmaPhiY= 0.1 mrad
- SigmaX= 5 mm
- SigmaY= 5 mm
+ SigmaThetaX= 0 mrad
+ SigmaPhiY= 0 mrad
+ SigmaX= 0 mm
+ SigmaY= 0 mm
MeanThetaX= 0 mrad
MeanPhiY= 0 mrad
MeanX= 0 mm
@@ -24,7 +34,7 @@ TwoBodyReaction
Target= 1H
Light= 2H
Heavy= 47Cr
- ExcitationEnergyLight= 0.47 MeV
+ ExcitationEnergyLight= 0.0 MeV
ExcitationEnergyHeavy= 0.0 MeV
%CrossSectionPath= CrossSections/48Cr_tot_0.txt CS48Cr
CrossSectionPath= flat.txt CS
diff --git a/Projects/E805_simu/Analysis.cxx b/Projects/E805_simu/Analysis.cxx
index 95867f102122a4e792b2ca886b2706e90c229eb7..cd83a6b1d1de656813ec48f3d266fca08d4c0dbf 100755
--- a/Projects/E805_simu/Analysis.cxx
+++ b/Projects/E805_simu/Analysis.cxx
@@ -43,6 +43,7 @@ void Analysis::Init(){
string beam= NPL::ChangeNameToG4Standard(reaction->GetNucleus1()->GetName());
string heavy_ejectile= NPL::ChangeNameToG4Standard(reaction->GetNucleus4()->GetName());
string light=NPL::ChangeNameToG4Standard(reaction->GetNucleus3()->GetName());
+ std::cout << light << " " << heavy_ejectile << std::endl;
string Reaction_pd_s = "48Cr(p,d)47Cr@1620";
string Reaction_pt_s = "48Cr(p,t)46Cr@1620";
@@ -132,8 +133,11 @@ void Analysis::TreatEvent(){
}
else
Energy = -1000;
- if(Energy > 0)
- Energy = LightTarget["deuteron"].EvaluateInitialEnergy(Energy,TargetThickness*0.5, ThetaNormalTarget);
+ if(Energy > 0){
+ // Energy = LightTarget["deuteron"].EvaluateInitialEnergy(Energy,TargetThickness*0.5, ThetaNormalTarget);
+ Energy = LightTarget["alpha"].EvaluateInitialEnergy(Energy,TargetThickness*0.5, ThetaNormalTarget);
+ }
+
// Evaluate energy using the thickness
diff --git a/Projects/E805_simu/MUST2_E805.detector b/Projects/E805_simu/MUST2_E805.detector
index b15b96f6da402074f8006b7a710ec8a128a79437..1437b43cda0060c66a573e36b31faba4111527d4 100644
--- a/Projects/E805_simu/MUST2_E805.detector
+++ b/Projects/E805_simu/MUST2_E805.detector
@@ -1,6 +1,7 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Target
- THICKNESS= 53.5 micrometer
+ %THICKNESS= 53.5 micrometer
+ THICKNESS= 0.000001 micrometer
ANGLE= 0 deg
RADIUS= 15 mm
MATERIAL= CH2
diff --git a/Projects/E805_simu/project.config b/Projects/E805_simu/project.config
deleted file mode 100644
index 0edda1e97a786835c74354cfe800128eb48e8261..0000000000000000000000000000000000000000
--- a/Projects/E805_simu/project.config
+++ /dev/null
@@ -1,8 +0,0 @@
-Project E805_simu
- AnalysisOutput= ../E805/data/NPRoot/Analysis/
- SimulationOutput= ../E805/data/NPRoot/Simulation/
- EnergyLoss= ../../Inputs/EnergyLoss/
- CalibrationOutput= ../E805/data/NPRoot/Calibration/
- Cuts= ../E805/data/NPRoot/Cuts/
-
-
diff --git a/Projects/e870/.DS_Store b/Projects/e870/.DS_Store
index f5ddf34afd94504d47c0de268f5e147d9a740da6..c4bac702c96a1b1e9866a310b007dec4c64658b1 100644
Binary files a/Projects/e870/.DS_Store and b/Projects/e870/.DS_Store differ
diff --git a/Projects/e870/Analysis.cxx b/Projects/e870/Analysis.cxx
index db027a980e60e9255dd3a01b51a7321385f2d95f..af55fe7ba7ff505fa1267fe2f1d7178ab39cf748 100644
--- a/Projects/e870/Analysis.cxx
+++ b/Projects/e870/Analysis.cxx
@@ -50,17 +50,18 @@ void Analysis::Init() {
Initial = new TInitialConditions();
ReactionConditions = new TReactionConditions();
}
-
+ // ExogamData = new TExogamData();
InitOutputBranch();
InitInputBranch();
// get MUST2 and Gaspard objects
M2 = (TMust2Physics*)m_DetectorManager->GetDetector("M2Telescope");
- if (!simulation)
- CATS = (TCATSPhysics*)m_DetectorManager->GetDetector("CATSDetector");
+ // if (!simulation)
+ // CATS = (TCATSPhysics*)m_DetectorManager->GetDetector("CATSDetector");
// get reaction information
reaction.ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
OriginalBeamEnergy = reaction.GetBeamEnergy();
+ std::cout << OriginalBeamEnergy << std::endl;
// target thickness
TargetThickness = m_DetectorManager->GetTargetThickness();
string TargetMaterial = m_DetectorManager->GetTargetMaterial();
@@ -70,6 +71,7 @@ void Analysis::Init() {
string beam = NPL::ChangeNameToG4Standard(reaction.GetNucleus1()->GetName());
cout << light << " " << beam << " " << TargetMaterial << " " << TargetThickness << endl;
LightTarget = NPL::EnergyLoss(light + "_" + TargetMaterial + ".G4table", "G4Table", 100);
+ LightMyl = NPL::EnergyLoss(light + "_Mylar.G4table", "G4Table", 100);
LightAl = NPL::EnergyLoss(light + "_Al.G4table", "G4Table", 100);
LightSi = NPL::EnergyLoss(light + "_Si.G4table", "G4Table", 100);
BeamTarget = NPL::EnergyLoss(beam + "_" + TargetMaterial + ".G4table", "G4Table", 100);
@@ -80,17 +82,6 @@ void Analysis::Init() {
<< "MeV " << endl;
reaction.SetBeamEnergy(FinalBeamEnergy);
- if (WindowsThickness) {
- cout << "Cryogenic target with windows" << endl;
- // BeamWindow = new NPL::EnergyLoss(beam + "_" + WindowsMaterial + ".G4table", "G4Table", 100);
- // LightWindow = new NPL::EnergyLoss(light + "_" + WindowsMaterial + ".G4table", "G4Table", 100);
- }
-
- else {
- BeamWindow = NULL;
- LightWindow = NULL;
- }
-
// initialize various parameters
Rand = TRandom3();
DetectorNumber = 0;
@@ -114,19 +105,12 @@ void Analysis::TreatEvent() {
ReInitValue();
double XTarget, YTarget;
TVector3 BeamDirection;
- if (!simulation) {
- XTarget = CATS->GetPositionOnTarget().X();
- YTarget = CATS->GetPositionOnTarget().Y();
- BeamDirection = CATS->GetBeamDirection();
- }
- else {
- XTarget = 0;
- YTarget = 0;
- BeamDirection = TVector3(0, 0, 1);
- // OriginalELab = ReactionConditions->GetKineticEnergy(0);
- // OriginalThetaLab = ReactionConditions->GetTheta(0);
- // BeamEnergy = ReactionConditions->GetBeamEnergy();
- }
+
+ // For simulation
+ XTarget = 0;
+ YTarget = 0;
+ BeamDirection = TVector3(0, 0, 1);
+
BeamImpact = TVector3(XTarget, YTarget, 0);
// determine beam energy for a randomized interaction point in target
// 1% FWHM randominastion (E/100)/2.35
@@ -136,7 +120,8 @@ void Analysis::TreatEvent() {
////////////////////////////////////////////////////////////////////////////
//////////////////////////////// LOOP on MUST2 ////////////////////////////
////////////////////////////////////////////////////////////////////////////
- for (unsigned int countMust2 = 0; countMust2 < M2->Si_E.size(); countMust2++) {
+ int M2_size = M2->Si_E.size();
+ for (unsigned int countMust2 = 0; countMust2 < M2_size; countMust2++) {
/************************************************/
// Part 0 : Get the usefull Data
// MUST2
@@ -168,9 +153,8 @@ void Analysis::TreatEvent() {
if (CsI_E_M2 > 0) {
// The energy in CsI is calculate form dE/dx Table because
Energy = CsI_E_M2;
- if (simulation) {
- Energy = LightAl.EvaluateInitialEnergy(Energy, 0.4 * micrometer, ThetaM2Surface);
- }
+ Energy = LightMyl.EvaluateInitialEnergy(Energy, 3 * micrometer, ThetaM2Surface);
+ Energy = LightAl.EvaluateInitialEnergy(Energy, 0.4 * micrometer, ThetaM2Surface);
Energy += Si_E_M2;
}
else {
@@ -180,7 +164,7 @@ void Analysis::TreatEvent() {
Energy = LightAl.EvaluateInitialEnergy(Energy, 0.4 * micrometer, ThetaM2Surface);
// Evaluate energy using the thickness
// Target Correction
- Energy = LightTarget.EvaluateInitialEnergy(Energy, TargetThickness * 0.5, ThetaNormalTarget);
+ Energy = LightTarget.EvaluateInitialEnergy(Energy, TargetThickness * 0.5, Theta);
// What is written in the tree
ThetaLab.push_back(Theta / deg);
@@ -195,6 +179,7 @@ void Analysis::TreatEvent() {
void Analysis::End() {}
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitOutputBranch() {
+ // RootOutput::getInstance()->GetTree()->Branch("Exogam", &ExogamData);
RootOutput::getInstance()->GetTree()->Branch("Ex", &Ex);
RootOutput::getInstance()->GetTree()->Branch("ELab", &ELab);
RootOutput::getInstance()->GetTree()->Branch("ThetaLab", &ThetaLab);
@@ -226,6 +211,10 @@ void Analysis::InitInputBranch() {
RootInput::getInstance()->GetChain()->SetBranchStatus("fRC_*", true);
RootInput::getInstance()->GetChain()->SetBranchAddress("ReactionConditions", &ReactionConditions);
}
+ // RootInput::getInstance()->GetChain()->SetBranchStatus("Exogam", true);
+ RootInput::getInstance()->GetChain()->SetBranchStatus("fExo*", true);
+ // RootInput::getInstance()->GetChain()->SetBranchAddress("Exogam", &ExogamData);
+ // RootInput::getInstance()->GetChain()->SetBranchStatus("ZDD", true);
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::ReInitValue() {
diff --git a/Projects/e870/Analysis.h b/Projects/e870/Analysis.h
index 450e560cbbfbdb7f318a4c2194bad69af6204108..753d9fcdbfb81f25e6efff8ba3727f67561884b3 100644
--- a/Projects/e870/Analysis.h
+++ b/Projects/e870/Analysis.h
@@ -30,7 +30,8 @@
#include "TReactionConditions.h"
#include "TMust2Physics.h"
#include "TMugastPhysics.h"
-#include "TCATSPhysics.h"
+// #include "TCATSPhysics.h"
+#include "TExogamData.h"
#include <TRandom3.h>
#include <TVector3.h>
#include <TMath.h>
@@ -64,6 +65,7 @@ class Analysis: public NPL::VAnalysis{
NPL::Reaction reaction;
// Energy loss table: the G4Table are generated by the simulation
NPL::EnergyLoss LightTarget;
+ NPL::EnergyLoss LightMyl;
NPL::EnergyLoss LightAl;
NPL::EnergyLoss LightSi;
NPL::EnergyLoss BeamTarget;
@@ -118,9 +120,11 @@ class Analysis: public NPL::VAnalysis{
// Branches and detectors
TMust2Physics* M2;
TMugastPhysics* MG;
- TCATSPhysics* CATS;
+ // TCATSPhysics* CATS;
bool simulation;
TInitialConditions* Initial;
TReactionConditions* ReactionConditions;
+ TExogamData* ExogamData;
+
};
#endif
diff --git a/Projects/e870/DetectorConfiguration/MUGAST_LISE.detector b/Projects/e870/DetectorConfiguration/MUGAST_LISE.detector
index 670c6f6a727af047791c4db3f01708df155e76b7..a47e82527e81b493cf82ff6f858cd6e155a9c60e 100644
--- a/Projects/e870/DetectorConfiguration/MUGAST_LISE.detector
+++ b/Projects/e870/DetectorConfiguration/MUGAST_LISE.detector
@@ -2,6 +2,7 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1
Target
THICKNESS= 26 micrometer
+ %THICKNESS= 0.00001 micrometer
ANGLE= 0 deg
RADIUS= 13 mm
MATERIAL= CH2
diff --git a/Projects/e870/DetectorConfiguration/MUGAST_LISE_CD2.detector b/Projects/e870/DetectorConfiguration/MUGAST_LISE_CD2.detector
index 020abacbd17a7139fd87ac009d6e6fa893c11ee5..c8cb44e76f2c0332331696758d46f3e691f524ee 100644
--- a/Projects/e870/DetectorConfiguration/MUGAST_LISE_CD2.detector
+++ b/Projects/e870/DetectorConfiguration/MUGAST_LISE_CD2.detector
@@ -36,6 +36,8 @@ M2Telescope
SILI= 0
CSI= 1
VIS= all
+ CsIOffset= 1
+
%%%%%%% Telescope 2 %%%%%%%
M2Telescope
X1_Y1= -114.9 9.68 291.84 mm
@@ -46,6 +48,7 @@ M2Telescope
SILI= 0
CSI= 1
VIS= all
+ CsIOffset= 1
%%%%%%% Telescope 3 %%%%%%%
M2Telescope
@@ -57,18 +60,21 @@ M2Telescope
SILI= 0
CSI= 1
VIS= all
+ CsIOffset= 1
+
+%%%%%%% Telescope 4 %%%%%%%
+M2Telescope
+ X1_Y1= 113.56 -13.18 292.11 mm
+ X1_Y128= 23.23 -13.37 328.15 mm
+ X128_Y1= 102.39 -105.49 263.59 mm
+ X128_Y128= 12.04 -105.69 299.63 mm
+ SI= 1
+ SILI= 0
+ CSI= 1
+ VIS= all
+ CsIOffset= 1
-% %%%%%%% Telescope 4 %%%%%%%
-% M2Telescope
-% X1_Y1= 113.56 -13.18 292.11 mm
-% X1_Y128= 23.23 -13.37 328.15 mm
-% X128_Y1= 102.39 -105.49 263.59 mm
-% X128_Y128= 12.04 -105.69 299.63 mm
-% SI= 1
-% SILI= 0
-% CSI= 1
-% VIS= all
%
% %%%%%%% Telescope 1 %%%%%%%
% M2Telescope
diff --git a/Projects/e870/launch_simu.sh b/Projects/e870/launch_simu.sh
index 9005ba4a2f61168db03b0c9e79258597071a5b67..da74d29302d83ba6a884c2c0d32d0d501b8cd07d 100755
--- a/Projects/e870/launch_simu.sh
+++ b/Projects/e870/launch_simu.sh
@@ -47,9 +47,9 @@
# All of the above:
npcompilation
-npsimulation -D DetectorConfiguration/MUGAST_LISE.detector -E reaction/10Bepalpha.reaction -O test -B run.mac
+npsimulation -D DetectorConfiguration/MUGAST_LISE.detector -E reaction/10Bed6Li.reaction -O test -B run.mac
-npanalysis -D DetectorConfiguration/MUGAST_LISE.detector -E reaction/10Bepalpha.reaction -T $NPTOOL/Outputs/Simulation/test.root SimulatedTree -O test
+npanalysis -D DetectorConfiguration/MUGAST_LISE.detector -E reaction/10Bed6Li.reaction -T $NPTOOL/Outputs/Simulation/test.root SimulatedTree -O test
diff --git a/Projects/e870/macros/DrawExInvariant.cxx b/Projects/e870/macros/DrawExInvariant.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..0d191812285a13367b9dcfdda01464145e69b49d
--- /dev/null
+++ b/Projects/e870/macros/DrawExInvariant.cxx
@@ -0,0 +1,166 @@
+void DrawDalitz() {
+ ////////////////////////////////////////////////////////////////////////////////
+ TLorentzVector LV_alpha, LV_proton, LV_fragment;
+ double mu = 931.5; // MeV/c2
+ double m_proton = mu + 7.289; // Mev/c2
+ double m_alpha = 4 * mu + 2.425; // Mev/c2
+ double m_fragment = 44 * mu - 3.755; // Mev/c2
+ double m_5Li = 5 * mu + 11.68; // MeV/c2
+
+ TChain* t = new TChain("PhysicsTree"); // SimulatedTree
+ t->Add("NPRootA/NPr0371_*a.root"); // Name of simulated file
+
+ // YOUR ALPHA AND Li CUT
+ TFile* fcuta = new TFile("CUT_e805/CUT_alpha_must.root", "read");
+ TCutG* cuta = (TCutG*)gDirectory->FindObjectAny("CUT_alpha_must");
+ TFile* fcutp = new TFile("CUT_e805/CUT_proton_must.root", "read");
+ TCutG* cutp = (TCutG*)gDirectory->FindObjectAny("CUT_proton_must");
+
+ ////////////////////////////////////////////////////////////////////////////////
+ TMust2Physics* M2 = new TMust2Physics();
+ std::vector<double>* ELab = 0;
+ std::vector<double>* ThetaLab = 0;
+ std::vector<double>* M2_X = 0;
+ std::vector<double>* M2_Y = 0;
+ std::vector<double>* M2_Z = 0;
+ TBranch* b_ELab;
+ TBranch* b_ThetaLab;
+ TBranch* b_M2_X;
+ TBranch* b_M2_Y;
+ TBranch* b_M2_Z;
+
+ float CATS1_X, CATS1_Y, CATS2_X, CATS2_Y, Xf, Yf;
+
+ t->SetBranchStatus("*", false);
+ t->SetBranchStatus("MUST2", true);
+ t->SetBranchAddress("MUST2", &M2);
+ t->SetBranchStatus("M2_ELab", true);
+ t->SetBranchAddress("M2_ELab", &ELab, &b_ELab);
+ t->SetBranchStatus("M2_ThetaLab", true);
+ t->SetBranchAddress("M2_ThetaLab", &ThetaLab, &b_ThetaLab);
+ t->SetBranchStatus("CATS1_X", "true");
+ t->SetBranchAddress("CATS1_X", &CATS1_X);
+ t->SetBranchStatus("CATS1_Y", "true");
+ t->SetBranchAddress("CATS1_Y", &CATS1_Y);
+ t->SetBranchStatus("CATS2_X", "true");
+ t->SetBranchAddress("CATS2_X", &CATS2_X);
+ t->SetBranchStatus("CATS2_Y", "true");
+ t->SetBranchAddress("CATS2_Y", &CATS2_Y);
+ t->SetBranchStatus("M2_X", "true");
+ t->SetBranchAddress("M2_X", &M2_X, &b_M2_X);
+ t->SetBranchStatus("M2_Y", "true");
+ t->SetBranchAddress("M2_Y", &M2_Y, &b_M2_Y);
+ t->SetBranchStatus("M2_Z", "true");
+ t->SetBranchAddress("M2_Z", &M2_Z, &b_M2_Z);
+
+ TH1F* hExTot = new TH1F("hExTot", "hExTot", 1000, -100, 100);
+ TH1F* hEx5Li = new TH1F("hEx5Li", "hEx5Li", 1000, -100, 100);
+
+ TH2F* hE1E2 = new TH2F("hE1E2", "hE1E2", 1000, 0, 500, 1000, 0, 500);
+
+ ////////////////////////////////////////////////////////////////////////////////
+ // int n_entries = t->GetEntries();
+ int n_entries = 1e7;
+ for (unsigned int i = 0; i < n_entries; i++) {
+ t->GetEntry(i);
+
+ double e_alpha;
+ double theta_alpha;
+ double phi_alpha;
+
+ double e_proton;
+ double theta_proton;
+ double phi_proton;
+
+ double e_fragment;
+ double theta_fragment;
+ double phi_fragment;
+
+ double xtarget = 0;
+ double ytarget = 0;
+ double CATSX_diff = 0;
+ double CATSY_diff = 0;
+
+ if (CATS1_X > -1500 && CATS2_X > -1500 && Xf > -1500) {
+ CATSX_diff = -(CATS2_X - CATS1_X);
+ xtarget = -Xf;
+ }
+ if (CATS1_Y > -1500 && CATS2_X > -1500 && Yf > -1500) {
+ ytarget = Yf;
+ CATSY_diff = CATS2_X - CATS1_X;
+ }
+
+ TVector3 BeamImpact(xtarget, ytarget, 0);
+ TVector3 BeamDirection(CATSX_diff, CATSY_diff, 497.1);
+
+ if (M2->Si_E.size() == 2) {
+ if ((cuta->IsInside(M2->CsI_E[0], M2->Si_E[0]) && cutp->IsInside(M2->CsI_E[1], M2->Si_E[1]))) {
+ // Particle1 = alpha
+ TVector3 PositionOfInteraction_alpha(M2_X->at(0), M2_Y->at(0), M2_Z->at(0));
+ TVector3 HitDirection_alpha = PositionOfInteraction_alpha - BeamImpact;
+ e_alpha = ELab->at(0);
+ theta_alpha = HitDirection_alpha.Angle(BeamDirection);
+ phi_alpha = HitDirection_alpha.Phi();
+
+ // Particle2 = proton
+ TVector3 PositionOfInteraction_proton(M2_X->at(1), M2_Y->at(1), M2_Z->at(1));
+ TVector3 HitDirection_proton = PositionOfInteraction_proton - BeamImpact;
+ e_proton = ELab->at(1);
+ theta_proton = HitDirection_proton.Angle(BeamDirection);
+ phi_proton = HitDirection_proton.Phi();
+ }
+ else if (cuta->IsInside(M2->CsI_E[1], M2->Si_E[1]) && cutp->IsInside(M2->CsI_E[0], M2->Si_E[0])) {
+ // Particle2 = alpha
+ TVector3 PositionOfInteraction_alpha(M2_X->at(1), M2_Y->at(1), M2_Z->at(1));
+ TVector3 HitDirection_alpha = PositionOfInteraction_alpha - BeamImpact;
+ e_alpha = ELab->at(1);
+ theta_alpha = HitDirection_alpha.Angle(BeamDirection);
+ phi_alpha = HitDirection_alpha.Phi();
+
+ // Particle1 = proton
+ TVector3 PositionOfInteraction_proton(M2_X->at(0), M2_Y->at(0), M2_Z->at(0));
+ TVector3 HitDirection_proton = PositionOfInteraction_proton - BeamImpact;
+ e_proton = ELab->at(0);
+ theta_proton = HitDirection_proton.Angle(BeamDirection);
+ phi_proton = HitDirection_proton.Phi();
+ }
+
+ ////////////////////////////////////////////////////////////////////////////////
+ double gamma_alpha = e_alpha / m_alpha + 1;
+ double beta_alpha = sqrt(1 - 1 / (pow(gamma_alpha, 2.)));
+ double p_alpha_mag = gamma_alpha * m_alpha * beta_alpha;
+ double etot_alpha = sqrt(pow(p_alpha_mag, 2) + pow(m_alpha, 2));
+
+ TVector3 p_alpha(sin(theta_alpha) * cos(phi_alpha), // px
+ sin(theta_alpha) * sin(phi_alpha), // py
+ cos(theta_alpha)); // pz
+ p_alpha.SetMag(p_alpha_mag);
+ LV_alpha.SetPxPyPzE(p_alpha.x(), p_alpha.y(), p_alpha.z(), etot_alpha);
+
+ ////////////////////////////////////////////////////////////////////////////////
+ double gamma_proton = e_proton / m_proton + 1;
+ double beta_proton = sqrt(1 - 1 / (pow(gamma_proton, 2.)));
+ double p_proton_mag = gamma_proton * m_proton * beta_proton;
+ double etot_proton = sqrt(pow(p_proton_mag, 2) + pow(m_proton, 2));
+
+ TVector3 p_proton(sin(theta_proton) * cos(phi_proton), // px
+ sin(theta_proton) * sin(phi_proton), // py
+ cos(theta_proton)); // pz
+ p_proton.SetMag(p_proton_mag);
+ LV_proton.SetPxPyPzE(p_proton.x(), p_proton.y(), p_proton.z(), etot_proton);
+ TLorentzVector LV_total = LV_alpha + LV_proton;
+
+ double ExTot = LV_total.Mag() - m_alpha - m_proton;
+ hExTot->Fill(ExTot);
+ if (Ex.size() >= 2) {
+ hEx->Fill(Ex.at(0));
+ hEx->Fill(Ex.at(1));
+ }
+ hE1E2->Fill(e_alpha, e_proton);
+ }
+ }
+ TCanvas* c1 = new TCanvas();
+ hExTot->Draw();
+ hEx->Draw("same");
+}
+~
diff --git a/Projects/e870/reaction/10Bed6Li.reaction b/Projects/e870/reaction/10Bed6Li.reaction
index ecf3e3413291709c44c3db79303f3791111ce757..dcf2f7459d8d02300451c34a0ae6c792e9e8bf01 100644
--- a/Projects/e870/reaction/10Bed6Li.reaction
+++ b/Projects/e870/reaction/10Bed6Li.reaction
@@ -1,7 +1,7 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Beam
Particle= 10Be
- Energy= 300 MeV
+ Energy= 150 MeV
SigmaEnergy= 21.6 MeV
MeanThetaX= 0 deg
MeanPhiY= 0 deg
diff --git a/Projects/e870/reaction/10Bepalpha.reaction b/Projects/e870/reaction/10Bepalpha.reaction
index 5c5bf0441c3f73067dd0c3ff1777594c8c3bcf5d..a4b5f079592b21fa8e3deaba9be601711f87b4f3 100644
--- a/Projects/e870/reaction/10Bepalpha.reaction
+++ b/Projects/e870/reaction/10Bepalpha.reaction
@@ -5,13 +5,13 @@ Beam
% Stuff to try and test:
% Perfect case:
- SigmaEnergy= 0 MeV
- SigmaX= 0 mm
- SigmaY= 0 mm
+ % SigmaEnergy= 0 MeV
+ % SigmaX= 0 mm
+ % SigmaY= 0 mm
% Realistic case?
- %SigmaEnergy= 22.8 MeV
- %SigmaX= 6.216 mm
- %SigmaY= 6.109 mm
+ SigmaEnergy= 22.8 MeV
+ SigmaX= 6.216 mm
+ SigmaY= 6.109 mm
% Playground:
%SigmaEnergy= ?? MeV
%SigmaX= ?? mm
@@ -21,13 +21,13 @@ Beam
MeanPhiY= 0 deg
MeanX= 0 mm
MeanY= 0 mm
- %SigmaThetaX= 0.6138 deg
- %SigmaPhiY= 0.3812 deg
+ SigmaThetaX= 0.6138 deg
+ SigmaPhiY= 0.3812 deg
%SigmaEnergy= 0 MeV
%ExcitationEnergy= 0 MeV
- SigmaThetaX= 0 mrad
- SigmaPhiY= 0 mrad
+ %SigmaThetaX= 0 mrad
+ %SigmaPhiY= 0 mrad
%MeanThetaX= 0 mrad
%MeanPhiY= 0 mrad
%MeanX= 0 mm