diff --git a/Cleopatra/SimTransfer.C b/Cleopatra/SimTransfer.C index faff24a..2b40557 100644 --- a/Cleopatra/SimTransfer.C +++ b/Cleopatra/SimTransfer.C @@ -189,9 +189,8 @@ void Transfer( detGeoTxt.Write("detGeo"); if( distList != NULL ) distList->Write("DWBA", 1); - if( dwbaExList != NULL ) { - dwbaExList_Used.Write("DWBA_ExList", 1); - } + if( dwbaExList != NULL ) dwbaExList_Used.Write("DWBA_ExList", 1); + TMacro idMacro; idMacro.AddLine(Form("%d", ID)); diff --git a/Cleopatra/SimTransfer2.C b/Cleopatra/SimTransfer2.C index 7f00720..25891d6 100644 --- a/Cleopatra/SimTransfer2.C +++ b/Cleopatra/SimTransfer2.C @@ -95,6 +95,9 @@ void Transfer( TObjArray * distList = nullptr; TMacro * dwbaExList = nullptr; TMacro * dwbaReactList = nullptr; + + TMacro dwbaExList_Used; + TMacro dwbaReactList_Used; if( distFile->IsOpen() ) { printf("\e[32m#################################### Load DWBA input : %s \e[0m\n", ptolemyRoot.Data()); @@ -102,6 +105,7 @@ void Transfer( 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"); int numEx = dwbaExList->GetListOfLines()->GetSize() - 1 ; @@ -115,6 +119,8 @@ void Transfer( for( int j = 0; j < numTransfer; j++){ if( reactionName.find( transfer[j].GetReactionName().Data() ) != std::string::npos) { std::string temp = dwbaExList->GetListOfLines()->At(i)->GetName(); + dwbaReactList_Used.AddLine((reactionName + " | " + std::to_string(j)).c_str()); + dwbaExList_Used.AddLine(temp.c_str()); if( temp[0] == '/' ) continue; std::vector tempStr = AnalysisLib::SplitStr(temp, " "); transfer[j].GetExList()->Add( atof(tempStr[0].c_str()), atof(tempStr[1].c_str()), 1.0, 0.00); @@ -155,7 +161,10 @@ 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); + dwbaReactList_Used.Write("DWBA_ReactionList", 1); + } TMacro hitMeaning; @@ -299,6 +308,15 @@ void Transfer( tree->Branch("yRecoil2", &yRecoil2, "yRecoil2/D"); tree->Branch("rhoRecoil2", &rhoRecoil2, "rhoRecoil2/D"); } + + //======= list of reaction used. + TMacro listOfReaction; + for( int i = 0; i < numTransfer ; i++){ + if( !listOfTransfer[i] ) continue; + listOfReaction.AddLine(Form("%2d | %s", i, transfer[i].GetReactionName_Latex().Data())); + } + + listOfReaction.Write("ListOfReactions"); //======= function for e-z plot for ideal case printf("++++ generate functions\n"); @@ -393,8 +411,6 @@ void Transfer( delete [] fx; delete fxList; - printf("============================== dasjdlasdj\n"); - // //--- cal modified thetaCM vs z // TObjArray * txList = new TObjArray(); // TGraph ** tx = new TGraph*[numEx]; diff --git a/Cleopatra/SimulationChecker.C b/Cleopatra/SimulationChecker.C new file mode 100644 index 0000000..4358047 --- /dev/null +++ b/Cleopatra/SimulationChecker.C @@ -0,0 +1,240 @@ +#include "TString.h" +#include "TFile.h" +#include "TTree.h" +#include "TH2.h" +#include "TStyle.h" +#include "TCanvas.h" +#include "TArc.h" +#include "TLine.h" +#include "TLatex.h" +#include "TMacro.h" + +#include "../Armory/AnalysisLib.h" +#include "../Armory/ClassDetGeo.h" +#include "../Armory/ClassReactionConfig.h" +#include "../Cleopatra/ClassIsotope.h" +#include "../Cleopatra/ClassTransfer.h" + +void SimulationChecker(TString transferRoot = "transfer.root"){ + + gStyle->SetOptStat(""); + gStyle->SetStatY(0.9); + gStyle->SetStatX(0.9); + gStyle->SetStatW(0.4); + gStyle->SetStatH(0.2); + gStyle->SetLabelSize(0.05, "XY"); + gStyle->SetTitleFontSize(0.1); + + TFile * file = new TFile(transferRoot, "read"); + TTree * tree = (TTree*) file->Get("tree"); + + TObjArray * fxList = (TObjArray *) file->FindObjectAny("EZCurve"); + + //*================================== List Of Reactions + TMacro * listOfReaction = (TMacro *) file->FindObjectAny("ListOfReactions"); + + int numReaction = listOfReaction->GetListOfLines()->GetSize(); + + std::vector reactionIDList(numReaction, -1); + std::vector reactionNameList(numReaction, ""); + for( int i = 0; i < numReaction; i++){ + std::vector tempStr = AnalysisLib::SplitStr(listOfReaction->GetListOfLines()->At(i)->GetName(), "|"); + reactionIDList[i] = atoi(tempStr[0].c_str()); + reactionNameList[i] = tempStr[1]; + } + + // for( int i = 0; i < numReaction; i++) printf("%2d | %s\n", reactionIDList[i], reactionNameList[i].Data()); + + //*================================== DWBA, if any + TMacro * dwbaReaction = (TMacro *) file->FindObjectAny("DWBA_ReactionList"); + TMacro * dwbaExList = (TMacro *) file->FindObjectAny("DWBA_ExList"); + // dwbaReaction->Print(); + // dwbaExList->Print(); + + + //*================================== reactionConfig + TMacro * reactionConfigTxt = (TMacro *) file->FindObjectAny("reactionConfig"); + + ReactionConfig reactionConfig(reactionConfigTxt); + + int nEvent = reactionConfig.numEvents; + printf("number of events generated : %d \n", nEvent); + + //=== clear the Ex in reactionCondig if DWBA_ReactionList + if( dwbaReaction || listOfReaction ){ + for( int i = 0; i < numReaction; i ++){ + reactionConfig.exList[reactionIDList[i]].Clear(); + + for( int j = 0; j < dwbaReaction->GetListOfLines()->GetSize(); j++){ + std::vector haha = AnalysisLib::SplitStr(dwbaReaction->GetListOfLines()->At(j)->GetName(), "|"); + + if( atoi(haha[1].c_str()) == reactionIDList[i] ) { + std::vector dudu = AnalysisLib::SplitStr(dwbaExList->GetListOfLines()->At(j+1)->GetName(), " "); + + reactionConfig.exList[reactionIDList[i]].Add(atof(dudu[0].c_str()), atof(dudu[1].c_str()), 1.0, 0.00); + + } + + } + + reactionConfig.Print(reactionIDList[i]); + } + } + + + + //*================================== detector Geometry + printf("=================================\n"); + printf(" loading detector Geometry.\n"); + TMacro * detGeoTxt = (TMacro *) file->FindObjectAny("detGeo"); + + DetGeo detGeo(detGeoTxt); + + detGeo.Print(true); + + printf("%f %f \n", detGeo.array[0].zMin, detGeo.array[0].zMax); + + printf("=================================\n"); + + //*================================== + + int padSize = 500; + int Div[2] = {5, numReaction}; + Int_t size[2] = {padSize,padSize}; ///x,y, single Canvas size + TCanvas * cCheck = new TCanvas("cCheck", "Check For Simulation", 0, 0, size[0]*Div[0], size[1]*Div[1]); + if(cCheck->GetShowEditor() )cCheck->ToggleEditor(); + if(cCheck->GetShowToolBar() )cCheck->ToggleToolBar(); + cCheck->Divide(Div[0],Div[1]); + + bool shownKELines = true; + TString gate = "hit == 1 && loop <= 1 && thetaCM > 10"; + + //^------------------------------- + cCheck->cd(1); + { + TH2F * hez = new TH2F("hez", Form("e-z [gated] @ %5.0f mm; z [mm]; e [MeV]", detGeo.array[0].firstPos), + 400, detGeo.array[0].zMin, detGeo.array[0].zMax, + 400, 0, 10); + + tree->Draw("e:z>>hez", gate + " && rID == 0", "colz"); + } + + //^------------------------------- + cCheck->cd(2); + { + TH2F * hRecoilXY = new TH2F("hRecoilXY", Form("RecoilXY [gated] @ %4.0f mm; X [mm]; Y [mm]", detGeo.aux[0].detPos ), + 400, -detGeo.aux[0].outerRadius, detGeo.aux[0].outerRadius, + 400, -detGeo.aux[0].outerRadius, detGeo.aux[0].outerRadius); + tree->Draw("yRecoil:xRecoil>>hRecoilXY", gate + " && rID == 0", "colz"); + TArc * detArc1 = new TArc(0,0, detGeo.aux[0].outerRadius); + detArc1->SetLineColor(kBlue-8); + detArc1->SetFillStyle(0); + detArc1->Draw("same"); + TArc * detArc2 = new TArc(0,0, detGeo.aux[0].innerRadius); + detArc2->SetLineColor(kBlue-8); + detArc2->SetFillStyle(0); + detArc2->Draw("same"); + + if( reactionConfig.beamX != 0. || reactionConfig.beamY != 0. ){ + TArc * arc = new TArc(reactionConfig.beamX, reactionConfig.beamY, 1); + arc->SetLineColor(2); + detArc1->SetFillStyle(0); + arc->Draw("same"); + } + } + //^------------------------------- + cCheck->cd(3); + { + // TH1F * hThetaCM[numEx]; + // TLegend * legend = new TLegend(0.8,0.2,0.99,0.8); + // double maxCount = 0; + // int startID = 0; // set the start ExID + // for( int i = startID; i < numEx; i++){ + // hThetaCM[i] = new TH1F(Form("hThetaCM%d", i), Form("thetaCM [gated] (ExID=%d); thetaCM [deg]; count", i), 200, thetaCMRange[0], thetaCMRange[1]); + // hThetaCM[i]->SetLineColor(i+1-startID); + // hThetaCM[i]->SetFillColor(i+1-startID); + // hThetaCM[i]->SetFillStyle(3000+i-startID); + // tree->Draw(Form("thetaCM>>hThetaCM%d", i), gate + Form("&& ExID==%d", i), ""); + // legend->AddEntry(hThetaCM[i], Form("Ex=%5.1f MeV", exList.ExList[i].Ex)); + // double max = hThetaCM[i]->GetMaximum(); + // if( max > maxCount ) maxCount = max; + // } + + // for( int i = startID; i < numEx; i++){ + // hThetaCM[i]->GetYaxis()->SetRangeUser(1, maxCount * 1.2); + // if( i == startID ) { + // hThetaCM[i]->Draw(); + // }else{ + // hThetaCM[i]->Draw("same"); + // } + // } + // legend->Draw(); + } + + //^------------------------------- + cCheck->cd(4); + { + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.06); + text.SetTextColor(2); + + text.DrawLatex(0., 0.9, reactionNameList[0]); + text.DrawLatex(0., 0.8, detGeo.Bfield > 0 ? "out of plan" : "into plan"); + text.SetTextColor(1); + text.DrawLatex(0., 0.7, "gate:"); + + text.SetTextColor(2); + //check gate text length, if > 30, break by "&&" + int ll = gate.Length(); + if( ll > 30 ) { + std::vector strList = AnalysisLib::SplitStr( (std::string) gate.Data(), "&&"); + for( int i = 0; i < strList.size(); i++){ + text.DrawLatex(0., 0.6 - 0.05*i, (TString) strList[i]); + } + }else{ + text.DrawLatex(0., 0.6, gate); + } + + if( reactionConfig.beamX != 0.0 || reactionConfig.beamY != 0.0 ){ + text.DrawLatex(0.0, 0.1, Form("Bema pos: (%4.1f, %4.1f) mm", reactionConfig.beamX, reactionConfig.beamY)); + } + } + //^------------------------------- + cCheck->cd(5); + { + TH1F * hExCal = new TH1F("hExCal", Form("calculated Ex [gated]; Ex [MeV]; count / %.2f keV", 10.), 400, -1, 3); + tree->Draw("ExCal>>hExCal", gate + " && rID == 0", ""); + Isotope hRecoil(reactionConfig.recoil[0].heavyA, reactionConfig.recoil[0].heavyZ); + double Sn = hRecoil.CalSp(0,1); + double Sp = hRecoil.CalSp(1,0); + double Sa = hRecoil.CalSp2(4,2); + double S2n = hRecoil.CalSp(0, 2); + + printf("Heavy recoil: %s \n", hRecoil.Name.c_str()); + printf("Sn : %f MeV/u \n", Sn); + printf("Sp : %f MeV/u \n", Sp); + printf("Sa : %f MeV/u \n", Sa); + printf("S2n : %f MeV/u \n", S2n); + + double yMax = hExCal->GetMaximum(); + TLine * lineSn = new TLine(Sn, 0, Sn, yMax); lineSn->SetLineColor(2); lineSn->Draw(""); + TLine * lineSp = new TLine(Sp, 0, Sp, yMax); lineSp->SetLineColor(4); lineSp->Draw("same"); + TLine * lineSa = new TLine(Sa, 0, Sa, yMax); lineSa->SetLineColor(6); lineSa->Draw("same"); + TLine * lineS2n = new TLine(S2n, 0, S2n, yMax); lineS2n->SetLineColor(8); lineS2n->Draw("same"); + + TLatex * text = new TLatex(); + text->SetTextFont(82); + text->SetTextSize(0.06); + text->SetTextColor(2); text->DrawLatex(Sn, yMax*0.9, "S_{n}"); + text->SetTextColor(4); text->DrawLatex(Sp, yMax*0.9, "S_{p}"); + text->SetTextColor(6); text->DrawLatex(Sa, yMax*0.9, "S_{a}"); + text->SetTextColor(8); text->DrawLatex(S2n, yMax*0.9, "S_{2n}"); + + } + + cCheck->Modified(); + cCheck->Update(); + +} \ No newline at end of file