diff --git a/.vscode/settings.json b/.vscode/settings.json index 4030e25..7cead8c 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -128,7 +128,8 @@ "test.C": "cpp", "SimTransfer.C": "cpp", "httpaccess.C": "cpp", - "httpcontrol.C": "cpp" + "httpcontrol.C": "cpp", + "SimTransfer2.C": "cpp" }, "better-comments.multilineComments": true, diff --git a/Armory/ClassDetGeo.h b/Armory/ClassDetGeo.h index 9f110d0..fbd968b 100644 --- a/Armory/ClassDetGeo.h +++ b/Armory/ClassDetGeo.h @@ -122,7 +122,7 @@ public: std::vector array; std::vector aux; - void Print( bool printAll = false) ; + void Print( bool printArray = false) ; private: @@ -234,7 +234,7 @@ inline bool DetGeo::LoadDetectorGeo(TMacro * macro, bool verbose){ } -inline void DetGeo::Print(bool printAll){ +inline void DetGeo::Print(bool printArray){ printf("#####################################################\n"); printf(" B-field : %8.2f T, %s\n", Bfield, Bfield > 0 ? "out of plan" : "into plan"); @@ -243,9 +243,9 @@ inline void DetGeo::Print(bool printAll){ printf(" z-Min : %8.2f mm\n", zMin); printf(" z-Max : %8.2f mm\n", zMax); - for( size_t i = 0; i < array.size() ; i++){ - printf("================================= %zu-th Detector Info (%s)\n", i, array[i].enable ? "enabled" : "disabled"); - if( printAll || array[i].enable ) { + if( printArray ) { + for( size_t i = 0; i < array.size() ; i++){ + printf("================================= %zu-th Detector Info (%s)\n", i, array[i].enable ? "enabled" : "disabled"); array[i].Print(); aux[i].Print(); } diff --git a/Armory/ClassReactionConfig.h b/Armory/ClassReactionConfig.h index 28a397c..e8cd448 100644 --- a/Armory/ClassReactionConfig.h +++ b/Armory/ClassReactionConfig.h @@ -108,8 +108,9 @@ public: bool isRedo; ///isReDo void SetReactionSimple(int beamA, int beamZ, - int targetA, int targetZ, - int recoilA, int recoilZ, float beamEnergy_AMeV, unsigned short ID); + int targetA, int targetZ, + int recoilA, int recoilZ, + float beamEnergy_AMeV); bool LoadReactionConfig(TString fileName); bool LoadReactionConfig(TMacro * macro); @@ -121,18 +122,21 @@ private: }; inline void ReactionConfig::SetReactionSimple(int beamA, int beamZ, - int targetA, int targetZ, - int recoilA, int recoilZ, float beamEnergy_AMeV, unsigned short ID){ + int targetA, int targetZ, + int recoilA, int recoilZ, + float beamEnergy_AMeV){ this->beamA = beamA; this->beamZ = beamZ; this->targetA = targetA; this->targetZ = targetZ; - this->recoil[ID].lightA = recoilA; - this->recoil[ID].lightZ = recoilZ; - recoil[ID].heavyA = this->beamA + this->targetA - recoil[ID].lightA; - recoil[ID].heavyZ = this->beamZ + this->targetZ - recoil[ID].lightZ; + this->recoil.push_back(Recoil()); + + this->recoil.back().lightA = recoilA; + this->recoil.back().lightZ = recoilZ; + recoil.back().heavyA = this->beamA + this->targetA - recoil.back().lightA; + recoil.back().heavyZ = this->beamZ + this->targetZ - recoil.back().lightZ; beamEnergy = beamEnergy_AMeV; beamEnergySigma = 0; @@ -260,8 +264,8 @@ inline void ReactionConfig::Print(int ID, bool withEx) const{ printf("================================= Number of recoil reactions : %zu\n", recoil.size()); for( size_t i = 0; i < recoil.size(); i ++ ){ - printf("------------------------------------------ Recoil-%zu\n", i); if( ID == i || ID < 0 ){ + printf("------------------------------------------ Recoil-%zu\n", i); recoil[i].Print(); if( withEx ) exList[i].Print(); } diff --git a/Cleopatra/Check_Simulation.C b/Cleopatra/Check_Simulation.C index d7c005f..7221cff 100644 --- a/Cleopatra/Check_Simulation.C +++ b/Cleopatra/Check_Simulation.C @@ -73,16 +73,16 @@ void Check_Simulation(TString filename = "transfer.root", TString gate = ExtractString(startLineNum+1, config); double elumRange = ExtractNumber(startLineNum+2, config); - vector thetaCMRange = doubleConvertor( StringToVector( ExtractString(startLineNum+3,config) )); + std::vector thetaCMRange = doubleConvertor( StringToVector( ExtractString(startLineNum+3,config) )); bool shownKELines = (ExtractString(startLineNum+4, config).Remove(4) == "true" ? true : false); bool isOverRideEx = (ExtractString(startLineNum+5, config).Remove(4) == "true" ? true : false); - vector oExRange = doubleConvertor( StringToVector ( ExtractString(startLineNum+6, config ))); + std::vector oExRange = doubleConvertor( StringToVector ( ExtractString(startLineNum+6, config ))); printf("%s \n", gate.Data()); ///==== config Canvas - vector plotConfig = StringToVector( ExtractString(startLineNum, config)); - vector canvas; + std::vector plotConfig = StringToVector( ExtractString(startLineNum, config)); + std::vector canvas; int colCount = 0; int colCount_new = 0; int rowCount = 1; @@ -100,7 +100,7 @@ void Check_Simulation(TString filename = "transfer.root", if( colCount == 0 ) colCount = colCount_new; //printf("plot row: %d, col: %d \n", rowCount, colCount); - vector Div = {colCount, rowCount}; + std::vector Div = {colCount, rowCount}; TFile * file = new TFile(filename, "read"); TTree * tree = (TTree*) file->Get("tree"); @@ -140,8 +140,8 @@ void Check_Simulation(TString filename = "transfer.root", DetGeo detGeo(detGeoTxt); Array array = detGeo.array[detGeoID]; - detGeo.PrintWithoutArray(); - array.PrintArray(); + detGeo.Print(); + array.Print(); printf("=================================\n"); @@ -177,13 +177,14 @@ void Check_Simulation(TString filename = "transfer.root", } }else{ - + numEx = exListTxt->GetListOfLines()->GetSize()-1; - for( int i = 0 ; i < numEx ; i++){ - double ex = atof(exListTxt->GetListOfLines()->At(i)->GetName()); + for( int i = 1 ; i <= numEx ; i++){ + std::vector tempStr = AnalysisLib::SplitStr(exListTxt->GetListOfLines()->At(i)->GetName(), " "); + double ex = atof(tempStr[0].c_str()); if( ex < ExRange[0] ) ExRange[0] = ex; if( ex > ExRange[1] ) ExRange[1] = ex; - exList.Add(ex, 0, 0, 0); + exList.Add(ex, atof(tempStr[1].c_str()), 1.0, 0.00); } } @@ -191,8 +192,8 @@ void Check_Simulation(TString filename = "transfer.root", exList.Print(); double dExRange = ExRange[1] - ExRange[0]; - ExRange[0] = ExRange[0] - dExRange * 0.1; - ExRange[1] = ExRange[1] + dExRange * 0.1; + ExRange[0] = ExRange[0] - 0.3 - dExRange * 0.1; + ExRange[1] = ExRange[1] + 0.3 + dExRange * 0.1; printf("Number of Ex states = %d \n", numEx); @@ -200,7 +201,13 @@ void Check_Simulation(TString filename = "transfer.root", //eRange by zRange and exList TransferReaction transfer; - transfer.SetReactionSimple( reactionConfig.beamA, reactionConfig.beamZ, reactionConfig.targetA, reactionConfig.targetZ, recoil.lightA, recoil.lightZ, reactionConfig.beamEnergy); + transfer.SetReactionSimple( reactionConfig.beamA, + reactionConfig.beamZ, + reactionConfig.targetA, + reactionConfig.targetZ, + recoil.lightA, + recoil.lightZ, + reactionConfig.beamEnergy); double QQ = transfer.GetCMTotalEnergy(); double gamm = transfer.GetReactionGamma(); @@ -208,7 +215,6 @@ void Check_Simulation(TString filename = "transfer.root", double slope = transfer.GetEZSlope( detGeo.Bfield); double eRange[2] = {0, 10}; - // double intercept = QQ/gamm - mass; eRange[1] = zRange[2] * slope; // printf("intercept of 0 MeV : %f MeV \n", intercept); @@ -272,14 +278,15 @@ void Check_Simulation(TString filename = "transfer.root", } if( pID == pRecoilXY ){ - TH2F * hRecoilXY = new TH2F("hRecoilXY", Form("RecoilXY [gated] @ %4.0f mm; X [mm]; Y [mm]", detGeo.recoilPos ), 400, -detGeo.recoilOuterRadius, detGeo.recoilOuterRadius, - 400, -detGeo.recoilOuterRadius, detGeo.recoilOuterRadius); + TH2F * hRecoilXY = new TH2F("hRecoilXY", Form("RecoilXY [gated] @ %4.0f mm; X [mm]; Y [mm]", detGeo.aux[detGeoID].detPos ), + 400, -detGeo.aux[detGeoID].outerRadius, detGeo.aux[detGeoID].outerRadius, + 400, -detGeo.aux[detGeoID].outerRadius, detGeo.aux[detGeoID].outerRadius); tree->Draw("yRecoil:xRecoil>>hRecoilXY", gate, "colz"); - TArc * detArc1 = new TArc(0,0, detGeo.recoilOuterRadius); + TArc * detArc1 = new TArc(0,0, detGeo.aux[detGeoID].outerRadius); detArc1->SetLineColor(kBlue-8); detArc1->SetFillStyle(0); detArc1->Draw("same"); - TArc * detArc2 = new TArc(0,0, detGeo.recoilInnerRadius); + TArc * detArc2 = new TArc(0,0, detGeo.aux[detGeoID].innerRadius); detArc2->SetLineColor(kBlue-8); detArc2->SetFillStyle(0); detArc2->Draw("same"); @@ -293,25 +300,27 @@ void Check_Simulation(TString filename = "transfer.root", } if( pID == pRecoilXY1 ){ - TH2F * hRecoilXY1 = new TH2F("hRecoilXY1", Form("RecoilXY-1 [gated] @ %4.0f mm; X [mm]; Y [mm]", detGeo.recoilPos1 ), 400, -detGeo.recoilOuterRadius, detGeo.recoilOuterRadius, - 400, -detGeo.recoilOuterRadius, detGeo.recoilOuterRadius); + TH2F * hRecoilXY1 = new TH2F("hRecoilXY1", Form("RecoilXY-1 [gated] @ %4.0f mm; X [mm]; Y [mm]", detGeo.aux[detGeoID].detPos1 ), + 400, -detGeo.aux[detGeoID].outerRadius, detGeo.aux[detGeoID].outerRadius, + 400, -detGeo.aux[detGeoID].outerRadius, detGeo.aux[detGeoID].outerRadius); tree->Draw("yRecoil1:xRecoil1>>hRecoilXY1", gate, "colz"); } if( pID == pRecoilXY2 ){ - TH2F * hRecoilXY2 = new TH2F("hRecoilXY2", Form("RecoilXY-2 [gated] @ %4.0f mm; X [mm]; Y [mm]", detGeo.recoilPos2 ), 400, -detGeo.recoilOuterRadius, detGeo.recoilOuterRadius, - 400, -detGeo.recoilOuterRadius, detGeo.recoilOuterRadius); + TH2F * hRecoilXY2 = new TH2F("hRecoilXY2", Form("RecoilXY-2 [gated] @ %4.0f mm; X [mm]; Y [mm]", detGeo.aux[detGeoID].detPos2 ), + 400, -detGeo.aux[detGeoID].outerRadius, detGeo.aux[detGeoID].outerRadius, + 400, -detGeo.aux[detGeoID].outerRadius, detGeo.aux[detGeoID].outerRadius); tree->Draw("yRecoil2:xRecoil2>>hRecoilXY2", gate, "colz"); } if( pID == pRecoilRZ ){ - TH2F * hRecoilRZ = new TH2F("hRecoilRZ", "RecoilR - Z [gated]; z [mm]; RecoilR [mm]", zRange[0], zRange[1], zRange[2], 400,0, detGeo.recoilOuterRadius); + TH2F * hRecoilRZ = new TH2F("hRecoilRZ", "RecoilR - Z [gated]; z [mm]; RecoilR [mm]", zRange[0], zRange[1], zRange[2], 400,0, detGeo.aux[detGeoID].outerRadius); tree->Draw("rhoRecoil:z>>hRecoilRZ", gate, "colz"); } if( pID == pRecoilRTR ){ FindRange("TB", gate, tree, recoilERange); - TH2F * hRecoilRTR = new TH2F("hRecoilRTR", "RecoilR - recoilE [gated]; recoil Energy [MeV]; RecoilR [mm]", 500, recoilERange[0], recoilERange[1], 500, 0, detGeo.recoilOuterRadius); + TH2F * hRecoilRTR = new TH2F("hRecoilRTR", "RecoilR - recoilE [gated]; recoil Energy [MeV]; RecoilR [mm]", 500, recoilERange[0], recoilERange[1], 500, 0, detGeo.aux[detGeoID].outerRadius); tree->Draw("rhoRecoil:TB>>hRecoilRTR", gate, "colz"); } @@ -391,7 +400,7 @@ void Check_Simulation(TString filename = "transfer.root", } if( pID == pRecoilRThetaCM ){ - TH2F * hRecoilRThetaCM = new TH2F("hRecoilRThetaCM", "RecoilR - thetaCM [gated]; thetaCM [deg]; RecoilR [mm]", 400, 0, 60, 400,0, detGeo.recoilOuterRadius); + TH2F * hRecoilRThetaCM = new TH2F("hRecoilRThetaCM", "RecoilR - thetaCM [gated]; thetaCM [deg]; RecoilR [mm]", 400, 0, 60, 400,0, detGeo.aux[detGeoID].outerRadius); tree->Draw("rhoRecoil:thetaCM>>hRecoilRThetaCM", gate, "colz"); } @@ -430,7 +439,7 @@ void Check_Simulation(TString filename = "transfer.root", } if( pID == pElum1XY ){ - TH2F * hElum1XY = new TH2F("hElum1XY", Form("Elum-1 XY [gated] @ %.0f mm ; X [mm]; Y [mm]", detGeo.elumPos1), 400, -elumRange, elumRange, 400, -elumRange, elumRange); + TH2F * hElum1XY = new TH2F("hElum1XY", Form("Elum-1 XY [gated] @ %.0f mm ; X [mm]; Y [mm]", detGeo.aux[detGeoID].elumPos1), 400, -elumRange, elumRange, 400, -elumRange, elumRange); tree->Draw("yElum1:xElum1>>hElum1XY", gate, "colz"); double count = hElum1XY->GetEntries(); diff --git a/Cleopatra/ClassTransfer.h b/Cleopatra/ClassTransfer.h index 29329da..33dc319 100644 --- a/Cleopatra/ClassTransfer.h +++ b/Cleopatra/ClassTransfer.h @@ -157,12 +157,14 @@ TransferReaction::TransferReaction(std::string configFile, unsigned short ID){ } TransferReaction::TransferReaction(int beamA, int beamZ, - int targetA, int targetZ, - int recoilA, int recoilZ, float beamEnergy_AMeV){ + int targetA, int targetZ, + int recoilA, int recoilZ, + float beamEnergy_AMeV){ Inititization(); SetReactionSimple(beamA, beamZ, targetA, targetZ, - recoilA, recoilZ, beamEnergy_AMeV); + recoilA, recoilZ, + beamEnergy_AMeV); } void TransferReaction::Inititization(){ @@ -243,12 +245,13 @@ void TransferReaction::SetIncidentEnergyAngle(double KEA, double theta, double p } void TransferReaction::SetReactionSimple(int beamA, int beamZ, - int targetA, int targetZ, - int recoilA, int recoilZ, float beamEnergy_AMeV){ + int targetA, int targetZ, + int recoilA, int recoilZ, + float beamEnergy_AMeV){ config.SetReactionSimple(beamA, beamZ, targetA, targetZ, - recoilA, recoilZ, beamEnergy_AMeV, 0); + recoilA, recoilZ, beamEnergy_AMeV); recoil = config.recoil[0]; @@ -281,6 +284,11 @@ void TransferReaction::SetReactionFromFile(std::string configFile, unsigned shor SetExA(config.beamEx); + if( ID > config.recoil.size() ){ + printf("Reaction Config only has %zu recoil settings. input is %u. Abort.\n", config.recoil.size(), ID); + return; + } + recoil = config.recoil[ID]; exList = config.exList[ID]; diff --git a/Cleopatra/SimTransfer.C b/Cleopatra/SimTransfer.C index 29adc52..faff24a 100644 --- a/Cleopatra/SimTransfer.C +++ b/Cleopatra/SimTransfer.C @@ -126,20 +126,29 @@ void Transfer( TFile * distFile = new TFile(ptolemyRoot, "read"); TObjArray * distList = nullptr; TMacro * dwbaExList = nullptr; + TMacro * dwbaReactList = nullptr; + + TMacro dwbaExList_Used; + if( distFile->IsOpen() ) { printf("--------- Found DWBA thetaCM distributions. Use the ExList from DWBA.\n"); distList = (TObjArray *) distFile->FindObjectAny("thetaCM_TF1"); // the function List - + dwbaExList = (TMacro *) distFile->FindObjectAny("ExList"); + dwbaExList_Used.AddLine(dwbaExList->GetListOfLines()->At(0)->GetName()); + dwbaReactList = (TMacro *) distFile->FindObjectAny("ReactionList"); exList.Clear(); - dwbaExList = (TMacro *) distFile->FindObjectAny("ExList"); int numEx = dwbaExList->GetListOfLines()->GetSize() - 1 ; for(int i = 1; i <= numEx ; i++){ - std::string temp = dwbaExList->GetListOfLines()->At(i)->GetName(); - if( temp[0] == '/' ) continue; - std::vector tempStr = AnalysisLib::SplitStr(temp, " "); - exList.Add( atof(tempStr[0].c_str()), atof(tempStr[1].c_str()), 1.0, 0.00); + std::string reactionName = dwbaReactList->GetListOfLines()->At(i-1)->GetName(); + if( reactionName.find( transfer.GetReactionName().Data() ) != std::string::npos) { + std::string temp = dwbaExList->GetListOfLines()->At(i)->GetName(); + dwbaExList_Used.AddLine(temp.c_str()); + if( temp[0] == '/' ) continue; + std::vector tempStr = AnalysisLib::SplitStr(temp, " "); + exList.Add( atof(tempStr[0].c_str()), atof(tempStr[1].c_str()), 1.0, 0.00); + } } }else{ @@ -180,7 +189,9 @@ void Transfer( detGeoTxt.Write("detGeo"); if( distList != NULL ) distList->Write("DWBA", 1); - if( dwbaExList != NULL ) dwbaExList->Write("DWBA_ExList", 1); + if( dwbaExList != NULL ) { + dwbaExList_Used.Write("DWBA_ExList", 1); + } TMacro idMacro; idMacro.AddLine(Form("%d", ID)); diff --git a/Cleopatra/makefile b/Cleopatra/makefile index 876a29f..d41d6fc 100644 --- a/Cleopatra/makefile +++ b/Cleopatra/makefile @@ -28,7 +28,7 @@ FindThetaCM: FindThetaCM.C FindThetaCM.h ClassTransfer.h ClassHelios.h ClassIsot SimTransfer: SimTransfer.C ClassTransfer.h ClassHelios.h ClassIsotope.h constant.h $(CC) SimTransfer.C -o SimTransfer `root-config --cflags --glibs` -SimTransfer2: SimTransfer2.C ClassTransfer.h ClassHelios.h ClassIsotope.h constant.h +SimTransfer2: SimTransfer2.C ClassTransfer.h ClassHelios.h ClassIsotope.h constant.h ../Armory/ClassReactionConfig.h ../Armory/ClassDetGeo.h $(CC) SimTransfer2.C -o SimTransfer2 `root-config --cflags --glibs` SimAlpha: SimAlpha.C ClassHelios.h