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/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/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