From 7a5107998a3c114811b6f79d6e009d29e9aa71e2 Mon Sep 17 00:00:00 2001 From: "Ryan@WorkStation" Date: Tue, 20 Feb 2024 18:15:16 -0500 Subject: [PATCH] update the Check_Simulation.C --- .gitignore | 2 +- Armory/ClassDetGeo.h | 39 +- Armory/ClassReactionConfig.h | 10 +- Cleopatra/Check_Simulation.C | 959 ++++++++++----------- Cleopatra/{alpha.C => SimAlpha.C} | 0 Cleopatra/{knockoutSim.C => SimKnockout.C} | 0 Cleopatra/{Transfer.h => SimTransfer.C} | 80 +- Cleopatra/Simulation_Helper.C | 3 +- Cleopatra/Transfer.C | 72 -- Cleopatra/makefile | 10 +- 10 files changed, 572 insertions(+), 603 deletions(-) rename Cleopatra/{alpha.C => SimAlpha.C} (100%) rename Cleopatra/{knockoutSim.C => SimKnockout.C} (100%) rename Cleopatra/{Transfer.h => SimTransfer.C} (89%) delete mode 100644 Cleopatra/Transfer.C diff --git a/.gitignore b/.gitignore index 6a646aa..abdd51d 100644 --- a/.gitignore +++ b/.gitignore @@ -25,4 +25,4 @@ Cleopatra/IsotopeShort Cleopatra/PlotSimulation Cleopatra/PlotTGraphTObjArray Cleopatra/SimAlpha -Cleopatra/Transfer \ No newline at end of file +Cleopatra/SimTransfer \ No newline at end of file diff --git a/Armory/ClassDetGeo.h b/Armory/ClassDetGeo.h index 0a239b4..847b94a 100644 --- a/Armory/ClassDetGeo.h +++ b/Armory/ClassDetGeo.h @@ -95,7 +95,9 @@ public: bool LoadDetectorGeo(TString fileName, bool verbose = true); bool LoadDetectorGeo(TMacro * macro, bool verbose = true); - void Print( bool printAll = true) const; + void PrintWithoutArray() ; + void Print( bool printAll = true) ; + private: @@ -203,24 +205,14 @@ inline bool DetGeo::LoadDetectorGeo(TMacro * macro, bool verbose){ } -inline void DetGeo::Print(bool printAll) const{ - +inline void DetGeo::PrintWithoutArray(){ printf("=====================================================\n"); - printf(" B-field: %8.2f T, Theta : %6.2f deg \n", Bfield, BfieldTheta); - if( BfieldTheta != 0.0 ) { - printf(" +---- field angle != 0 is not supported!!! \n"); - } + printf(" B-field: %8.2f T, %s\n", Bfield, Bfield > 0 ? "out of plan" : "into plan"); + printf(" B-field Theta : %6.2f deg \n", BfieldTheta); + + if( BfieldTheta != 0.0 ) printf(" +---- field angle != 0 is not supported!!! \n"); printf(" Recoil detector pos: %8.2f mm, radius: %6.2f - %6.2f mm \n", recoilPos, recoilInnerRadius, recoilOuterRadius); - for( int i = 0; i < 2 ; i++){ - - if( printAll || array[i].enable ) { - printf("-----------------------------------%d-th Detector Position \n", i); - array[i].PrintArray(); - } - - } - if( elumPos1 != 0 || elumPos2 != 0 || recoilPos1 != 0 || recoilPos2 != 0){ printf("=================================== Auxillary/Imaginary Detectors\n"); } @@ -229,6 +221,21 @@ inline void DetGeo::Print(bool printAll) const{ if( recoilPos1 != 0 ) printf(" Recoil 1 pos.: %f mm \n", recoilPos1); if( recoilPos2 != 0 ) printf(" Recoil 2 pos.: %f mm \n", recoilPos2); printf("=====================================================\n"); +} + +inline void DetGeo::Print(bool printAll){ + + PrintWithoutArray(); + + for( int i = 0; i < 2 ; i++){ + + if( printAll || array[i].enable ) { + printf("-----------------------------------%d-th Detector Position \n", i); + array[i].PrintArray(); + } + + } + printf("=====================================================\n"); } diff --git a/Armory/ClassReactionConfig.h b/Armory/ClassReactionConfig.h index b24c02d..e33a00d 100644 --- a/Armory/ClassReactionConfig.h +++ b/Armory/ClassReactionConfig.h @@ -114,7 +114,7 @@ public: bool LoadReactionConfig(TString fileName); bool LoadReactionConfig(TMacro * macro); - void Print() const; + void Print(int ID = -1, bool withEx = true) const; private: @@ -237,7 +237,7 @@ inline bool ReactionConfig::LoadReactionConfig(TMacro * macro){ return true; } -inline void ReactionConfig::Print() const{ +inline void ReactionConfig::Print(int ID, bool withEx) const{ printf("=====================================================\n"); @@ -261,8 +261,10 @@ inline void ReactionConfig::Print() const{ for( int i = 0; i < 2; i ++ ){ printf("------------------------------ Recoil-%d\n", i); - recoil[i].Print(); - exList[i].Print(); + if( ID == i || ID < 0 ){ + recoil[i].Print(); + if( withEx ) exList[i].Print(); + } } diff --git a/Cleopatra/Check_Simulation.C b/Cleopatra/Check_Simulation.C index 57eb9cb..4d78bd8 100644 --- a/Cleopatra/Check_Simulation.C +++ b/Cleopatra/Check_Simulation.C @@ -22,6 +22,7 @@ #include "../Armory/ClassDetGeo.h" #include "../Armory/ClassReactionConfig.h" #include "../Cleopatra/ClassIsotope.h" +#include "../Cleopatra/ClassTransfer.h" double * FindRange(TString branch, TString gate, TTree * tree, double output[2]); double ExtractNumber(int index, TMacro * macro); @@ -49,7 +50,9 @@ enum plotID { pEZ, /// 0 pElum1RThetaCM, /// 16 pEmpty }; /// 17 plotID StringToPlotID(TString str); -void Check_Simulation(TString filename = "transfer1.root", + +//*===================================================== +void Check_Simulation(TString filename = "transfer.root", TString configFile = "../working/Check_Simulation_Config.txt", Int_t padSize = 500, bool outputCanvas = false){ @@ -60,12 +63,12 @@ void Check_Simulation(TString filename = "transfer1.root", int numLine = config->GetListOfLines()->GetSize(); int startLineNum = 0; for( int i = 0; i < numLine ; i++){ - TString haha = config->GetListOfLines()->At(i)->GetName(); - haha.Remove(4); - if( haha != "////" ) { - startLineNum = i; - break; - } + TString haha = config->GetListOfLines()->At(i)->GetName(); + haha.Remove(4); + if( haha != "////" ) { + startLineNum = i; + break; + } } TString gate = ExtractString(startLineNum+1, config); @@ -75,8 +78,7 @@ void Check_Simulation(TString filename = "transfer1.root", bool isOverRideEx = (ExtractString(startLineNum+5, config).Remove(4) == "true" ? true : false); vector oExRange = doubleConvertor( StringToVector ( ExtractString(startLineNum+6, config ))); - printf("%s \n", gate.Data()); - + printf("%s \n", gate.Data()); ///==== config Canvas vector plotConfig = StringToVector( ExtractString(startLineNum, config)); @@ -85,481 +87,466 @@ void Check_Simulation(TString filename = "transfer1.root", int colCount_new = 0; int rowCount = 1; for( int i = 0; i < (int) plotConfig.size(); i++){ - if( plotConfig[i] == "break" ) { - rowCount ++; - if( colCount_new > colCount ) colCount = colCount_new; - colCount_new = 0; - continue; - } - canvas.push_back( StringToPlotID(plotConfig[i])); - colCount_new ++; + if( plotConfig[i] == "break" ) { + rowCount ++; + if( colCount_new > colCount ) colCount = colCount_new; + colCount_new = 0; + continue; + } + canvas.push_back( StringToPlotID(plotConfig[i])); + colCount_new ++; } if( colCount == 0 ) colCount = colCount_new; - ///printf("plot row: %d, col: %d \n", rowCount, colCount); + //printf("plot row: %d, col: %d \n", rowCount, colCount); vector Div = {colCount, rowCount}; - - + TFile * file = new TFile(filename, "read"); TTree * tree = (TTree*) file->Get("tree"); - TObjArray * fxList = (TObjArray *) file->FindObjectAny("fxList"); - TObjArray * txList = (TObjArray *) file->FindObjectAny("txList"); + TObjArray * fxList = (TObjArray *) file->FindObjectAny("EZCure"); + TObjArray * txList = (TObjArray *) file->FindObjectAny("EZ_thetaCM"); - //================== reactionConfig + 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); + + //*================== detGeoID + TMacro * detGeoIDTxt = (TMacro *) file->FindObjectAny("detGeoID"); + int detGeoID = atoi(detGeoIDTxt->GetListOfLines()->At(0)->GetName()); + + //*================== reactionConfig TMacro * reactionConfigTxt = (TMacro *) file->FindObjectAny("reactionConfig"); - TString Reaction=reactionConfigTxt->GetName(); - - ReactionConfig reactionConfig; - reactionConfig.LoadReactionConfig(reactionConfigTxt); + TString Reaction = reactionConfigTxt->GetName(); + ReactionConfig reactionConfig(reactionConfigTxt); + Recoil recoil = reactionConfig.recoil[detGeoID]; + int nEvent = reactionConfig.numEvents; - printf("number of events generated : %d \n", nEvent); + printf("number of events generated : %d \n", nEvent); - double xBeam = reactionConfig.beamX; - double yBeam = reactionConfig.beamY; - printf(" beam position : (%5.2f, %5.2f) mm \n", xBeam, yBeam); + reactionConfig.Print(detGeoID, false); - 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); - - double eRange[2] = {0, 10}; - double zRange[3] = {400, -1000, 1000}; /// zRange[0] = nBin - double recoilERange[2]; - vector exList; - double ExRange[2]; - - //================================== detetcor Geometry + //*================================== detetcor Geometry printf("=================================\n"); printf(" loading detector Geometry.\n"); TMacro * detGeoTxt = (TMacro *) file->FindObjectAny("detGeo"); - DetGeo detGeo; - detGeo.LoadDetectorGeo(detGeoTxt); + DetGeo detGeo(detGeoTxt); + Array array = detGeo.array[detGeoID]; - Array array; - if( detGeo.use2ndArray){ - array = detGeo.array2; - }else{ - array = detGeo.array1; - } - - double field = detGeo.Bfield; - TString fdmsg = field > 0 ? "out of plan" : "into plan"; - TString msg2; - msg2.Form("field = %.2f T, %s", field, fdmsg.Data()); - - double prepDist = array.detPerpDist; - double length = array.detLength; - - double posRecoil = detGeo.recoilPos; - double rhoRecoilIn = detGeo.recoilInnerRadius; - double rhoRecoilOut = detGeo.recoilOuterRadius; - - double posRecoil1 = detGeo.recoilPos1; - double posRecoil2 = detGeo.recoilPos2; - - vector pos = array.detPos; - - float firstPos = array.firstPos; - int rDet = array.nDet; - int cDet = array.mDet; - - double elum1 = detGeo.elumPos1; - - printf("number of row-Det : %d \n", rDet); - printf("number of col-Det : %d \n", cDet); - - for(int i = 0; i < rDet ; i++){ - if( firstPos > 0 ){ - printf("%d, %10.2f mm - %10.2f mm \n", i, pos[i], pos[i] + length); - }else{ - printf("%d, %10.2f mm - %10.2f mm \n", i, pos[i] - length , pos[i]); - } - } + detGeo.PrintWithoutArray(); + array.PrintArray(); printf("=================================\n"); - int numDet = rDet * cDet; + int numDet = array.nDet * array.mDet ; + double zRange[3] = {400, -1000, 1000}; /// zRange[0] = nBin zRange[1] = array.zMin - 50; zRange[2] = array.zMax + 50; printf(" zRange : %f - %f \n", zRange[1], zRange[2]); printf("=================================\n"); - //========================================= Ex List; - printf(" loading Ex list\n"); - - TMacro * exListMacro = (TMacro *) file->FindObjectAny("ExList"); - int numEx = exListMacro->GetListOfLines()->GetSize() - 1 ; - for(int i = 1; i <= numEx ; i++){ - string temp = exListMacro->GetListOfLines()->At(i)->GetName(); - if( temp[0] == '#' ) break; - if( temp[0] == '/' ) continue; - vector tempStr = AnalysisLib::SplitStr(temp, " "); - printf("%d | %s \n", i, tempStr[0].c_str()); - exList.push_back(atof(tempStr[0].c_str())); - } - - double exSpan = exList.back() - exList[0]; - - const int nExID = exList.size(); - printf("========= number of excited states : %d \n", nExID); + //*========================================= Ex List; + double ExRange[2]; + int numEx = 0; + ExcitedEnergies exList; - ExRange[0] = exList[0] - exSpan * 0.2; - ExRange[1] = exList.back() + exSpan * 0.2; - if( isOverRideEx ) { - ExRange[0] = oExRange[0]; - ExRange[1] = oExRange[1]; - } - - printf("=================================\n"); - - //========================================= reaction parameters - - printf(" loading reaction parameters \n"); - TMacro * reactionData = (TMacro *) file->FindObjectAny("reactionData"); - - double mass = ExtractNumber(0, reactionData); - double q = ExtractNumber(1, reactionData); - double beta = ExtractNumber(2, reactionData); - double Et = ExtractNumber(3, reactionData); - double massB = ExtractNumber(4, reactionData); - double alpha = ExtractNumber(5, reactionData); - - double gamm = 1./TMath::Sqrt(1-beta*beta); - double slope = alpha * beta; - printf("\tmass-b : %f MeV/c2 \n", mass); - printf("\tcharge-b : %f \n", q); - printf("\tE-total : %f MeV \n", Et); - printf("\tmass-B : %f MeV/c2 \n", massB); - printf("\tbeta : %f \n", beta); - printf("\tslope : %f MeV/mm \n", slope); - printf("=================================\n"); - - //=================================== calculate Ranges - - //eRange by zRange and exList - double QQ = (Et*Et + mass*mass - (massB-exList[0])*(massB-exList[0]))/2/Et; - double intercept = QQ/gamm - mass; - eRange[1] = intercept + zRange[2] * slope; - ///printf("intercept of 0 MeV : %f MeV \n", intercept); - ///printf("eRange 0 MeV : %f MeV \n", eRange[1]); - - //thetaCMRange - ///double momentum = sqrt(( Et*Et - pow(mass + massB - exList[0],2)) * ( Et*Et - pow(mass - massB + exList[0],2)))/2/Et; - ///double thetaMax = acos( (beta * QQ- alpha / gamm * zRange[2])/momentum) * TMath::RadToDeg(); - ///thetaCMRange[1] = (int) TMath::Ceil(thetaMax/10.)*10; - ///printf(" momentum : %f \n", momentum); - ///printf(" thetaCM Max : %f \n", thetaMax); - ///printf(" thetaCM Range : %d \n", thetaCMRange[1]); - - //=================================================== - printf("============================== Gate\n"); - printf("gate : %s\n", gate.Data()); - printf("====================================\n"); + // if DEBA_ExList exist, use this, else use the recoil ExList + TMacro * exListTxt = (TMacro *) file->FindObjectAny("DWBA_ExList"); - 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]); + ExRange[0] = 9999999; + ExRange[1] = -9999999; + + if( exListTxt == nullptr ){ + + exList = reactionConfig.exList[detGeoID]; + numEx = exList.ExList.size(); + + for( size_t i = 0; i < numEx; i++ ){ + double ex = exList.ExList[i].Ex; + if( ex < ExRange[0] ) ExRange[0] = ex; + if( ex > ExRange[1] ) ExRange[1] = ex; + } + + }else{ + + numEx = exListTxt->GetListOfLines()->GetSize(); + for( int i = 0 ; i < numEx ; i++){ + double ex = atof(exListTxt->GetListOfLines()->At(i)->GetName()); + if( ex < ExRange[0] ) ExRange[0] = ex; + if( ex > ExRange[1] ) ExRange[1] = ex; + exList.Add(ex, 0, 0, 0); + } + + } + + double dExRange = ExRange[1] - ExRange[0]; + ExRange[0] = ExRange[0] - dExRange * 0.1; + ExRange[1] = ExRange[1] + dExRange * 0.1; - for( int i = 1; i <= Div[0]*Div[1] ; i++){ - cCheck->cd(i); - cCheck->cd(i)->SetGrid(); - if( canvas[i-1] == pThetaCM ) { - cCheck->cd(i)->SetGrid(0,0); - cCheck->cd(i)->SetLogy(); - } - if( canvas[i-1] == pHitID ){ - cCheck->cd(i)->SetLogy(); - } - plotID pID = canvas[i-1]; - ///######################################## - if( pID == pEZ){ - TH2F * hez = new TH2F("hez", Form("e-z [gated] @ %5.0f mm; z [mm]; e [MeV]", firstPos), zRange[0], zRange[1], zRange[2], 400, eRange[0], eRange[1]); - tree->Draw("e:z>>hez", gate, "colz"); - if( shownKELines){ - for( int i = 0; i < nExID ; i++){ - fxList->At(i)->Draw("same"); - } - } - } - if( pID == pRecoilXY ){ - TH2F * hRecoilXY = new TH2F("hRecoilXY", Form("RecoilXY [gated] @ %4.0f mm; X [mm]; Y [mm]", posRecoil ), 400, -rhoRecoilOut, rhoRecoilOut, 400,-rhoRecoilOut, rhoRecoilOut); - tree->Draw("yRecoil:xRecoil>>hRecoilXY", gate, "colz"); - TArc * detArc1 = new TArc(0,0, rhoRecoilOut); - detArc1->SetLineColor(kBlue-8); - detArc1->SetFillStyle(0); - detArc1->Draw("same"); - TArc * detArc2 = new TArc(0,0, rhoRecoilIn); - detArc2->SetLineColor(kBlue-8); - detArc2->SetFillStyle(0); - detArc2->Draw("same"); - - if( xBeam != 0. || yBeam != 0. ){ - TArc * arc = new TArc(xBeam, yBeam, 1); - arc->SetLineColor(2); - detArc1->SetFillStyle(0); - arc->Draw("same"); - } - } - if( pID == pRecoilXY1 ){ - TH2F * hRecoilXY1 = new TH2F("hRecoilXY1", Form("RecoilXY-1 [gated] @ %4.0f mm; X [mm]; Y [mm]", posRecoil1 ), 400, -rhoRecoilOut, rhoRecoilOut, 400,-rhoRecoilOut, rhoRecoilOut); - 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]", posRecoil2 ), 400, -rhoRecoilOut, rhoRecoilOut, 400,-rhoRecoilOut, rhoRecoilOut); - 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, rhoRecoilOut); - 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, rhoRecoilOut); - tree->Draw("rhoRecoil:TB>>hRecoilRTR", gate, "colz"); - } - if( pID == pTDiffZ ){ - double tDiffRange [2]; - FindRange("t-tB", gate, tree, tDiffRange); - TH2F * hTDiffZ = new TH2F("hTDiffZ", "time(Array) - time(Recoil) vs Z [gated]; z [mm]; time diff [ns]", zRange[0], zRange[1], zRange[2], 500, tDiffRange[0], tDiffRange[1]); - tree->Draw("t - tB : z >> hTDiffZ", gate, "colz"); - } - if( pID == pThetaCM ){ - TH1F * hThetaCM[nExID]; - 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 < nExID; 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[i])); - double max = hThetaCM[i]->GetMaximum(); - if( max > maxCount ) maxCount = max; + printf("Number of Ex states = %d \n", numEx); + + //*=================================== calculate Ranges + //eRange by zRange and exList + + TransferReaction transfer; + transfer.SetReactionSimple( reactionConfig.beamA, reactionConfig.beamZ, reactionConfig.targetA, reactionConfig.targetZ, recoil.lightA, recoil.lightZ, reactionConfig.beamEnergy); + + double QQ = transfer.GetCMTotalEnergy(); + double gamm = transfer.GetReactionGamma(); + double mass = transfer.GetMass_b(); + 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); + printf("eRange 0 MeV : %f MeV \n", eRange[1]); + + double dERange = eRange[1] - eRange[0]; + + eRange[0] = eRange[0] - dERange * 0.1; + eRange[1] = eRange[1] + dERange * 0.1; + + + //thetaCMRange + double momentum = transfer.GetMomentumbCM(); + double beta = transfer.GetReactionBeta(); + double alpha = slope / beta; + double thetaMax = acos( (beta * QQ- alpha / gamm * zRange[2])/momentum) * TMath::RadToDeg(); + thetaCMRange[1] = (int) TMath::Ceil(thetaMax/10.)*10; + ///printf(" momentum : %f \n", momentum); + ///printf(" thetaCM Max : %f \n", thetaMax); + ///printf(" thetaCM Range : %d \n", thetaCMRange[1]); + + + double recoilERange[2] = {0, 100}; + + //=================================================== + printf("============================== Gate\n"); + printf("gate : %s\n", gate.Data()); + printf("====================================\n"); + + 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]); + + for( int i = 1; i <= Div[0]*Div[1] ; i++){ + cCheck->cd(i); + cCheck->cd(i)->SetGrid(); + + if( canvas[i-1] == pThetaCM ) { + cCheck->cd(i)->SetGrid(0,0); + cCheck->cd(i)->SetLogy(); + } + + if( canvas[i-1] == pHitID ){ + cCheck->cd(i)->SetLogy(); + } + + plotID pID = canvas[i-1]; + + ///######################################## + if( pID == pEZ){ + TH2F * hez = new TH2F("hez", Form("e-z [gated] @ %5.0f mm; z [mm]; e [MeV]", array.firstPos), zRange[0], zRange[1], zRange[2], + 400, eRange[0], eRange[1]); + tree->Draw("e:z>>hez", gate, "colz"); + if( shownKELines){ + for( int i = 0; i < numEx ; i++){ + fxList->At(i)->Draw("same"); } - - for( int i = startID; i < nExID; i++){ - hThetaCM[i]->GetYaxis()->SetRangeUser(1, maxCount * 1.2); - if( i == startID ) { - hThetaCM[i]->Draw(); - }else{ - hThetaCM[i]->Draw("same"); - } - } - legend->Draw(); } - if( pID == pThetaCM_Z ){ - TH2F *hThetaCM_Z = new TH2F("hThetaCM_Z","ThetaCM vs Z ; Z [mm]; thetaCM [deg]",zRange[0], zRange[1], zRange[2], 200, thetaCMRange[0], thetaCMRange[1]); - tree->Draw("thetaCM:z>>hThetaCM_Z",gate,"col"); - if( shownKELines){ - for( int i = 0; i < nExID ; i++){ - txList->At(i)->Draw("same"); - } - } + } + 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); + tree->Draw("yRecoil:xRecoil>>hRecoilXY", gate, "colz"); + TArc * detArc1 = new TArc(0,0, detGeo.recoilOuterRadius); + detArc1->SetLineColor(kBlue-8); + detArc1->SetFillStyle(0); + detArc1->Draw("same"); + TArc * detArc2 = new TArc(0,0, detGeo.recoilInnerRadius); + 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"); } - if( pID == pExCal ){ - TH1F * hExCal = new TH1F("hExCal", Form("calculated Ex [gated]; Ex [MeV]; count / %.2f keV", (ExRange[1]-ExRange[0])/400.*1000), 400, ExRange[0], ExRange[1]); - tree->Draw("ExCal>>hExCal", gate, ""); - Isotope hRecoil(reactionConfig.recoilHeavyA, reactionConfig.recoilHeavyZ); - 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); + 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); + tree->Draw("yRecoil1:xRecoil1>>hRecoilXY1", gate, "colz"); + } - 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"); + 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); + tree->Draw("yRecoil2:xRecoil2>>hRecoilXY2", gate, "colz"); + } - 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}"); - + 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); + 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); + tree->Draw("rhoRecoil:TB>>hRecoilRTR", gate, "colz"); + } + + if( pID == pTDiffZ ){ + double tDiffRange [2]; + FindRange("t-tB", gate, tree, tDiffRange); + TH2F * hTDiffZ = new TH2F("hTDiffZ", "time(Array) - time(Recoil) vs Z [gated]; z [mm]; time diff [ns]", zRange[0], zRange[1], zRange[2], 500, tDiffRange[0], tDiffRange[1]); + tree->Draw("t - tB : z >> hTDiffZ", gate, "colz"); + } + + if( pID == pThetaCM ){ + 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; } - if( pID == pRecoilRThetaCM ){ - TH2F * hRecoilRThetaCM = new TH2F("hRecoilRThetaCM", "RecoilR - thetaCM [gated]; thetaCM [deg]; RecoilR [mm]", 400, 0, 60, 400,0, rhoRecoilOut); - tree->Draw("rhoRecoil:thetaCM>>hRecoilRThetaCM", gate, "colz"); - } - if( pID == pArrayXY ){ - TH2F * hArrayXY = new TH2F("hArrayXY", "Array-XY [gated]; X [mm]; Y [mm]", 400, -prepDist*1.5, prepDist*1.5, 400, -prepDist*1.5, prepDist*1.5); - tree->Draw("yArray:xArray>>hArrayXY", gate, "colz"); - } - if( pID == pInfo ){ - TLatex text; - text.SetNDC(); - text.SetTextFont(82); - text.SetTextSize(0.06); - text.SetTextColor(2); - text.DrawLatex(0., 0.9, Reaction); - text.DrawLatex(0., 0.8, msg2); - 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 ) { - vector strList = AnalysisLib::SplitStr( (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); - } + 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(); + } - if( xBeam != 0.0 || yBeam != 0.0 ){ - text.DrawLatex(0.0, 0.1, Form("Bema pos: (%4.1f, %4.1f) mm", xBeam, yBeam)); - } + if( pID == pThetaCM_Z ){ + TH2F *hThetaCM_Z = new TH2F("hThetaCM_Z","ThetaCM vs Z ; Z [mm]; thetaCM [deg]",zRange[0], zRange[1], zRange[2], 200, thetaCMRange[0], thetaCMRange[1]); + tree->Draw("thetaCM:z>>hThetaCM_Z",gate,"col"); + if( shownKELines){ + for( int i = 0; i < numEx ; i++){ + txList->At(i)->Draw("same"); + } + } + } + + if( pID == pExCal ){ + TH1F * hExCal = new TH1F("hExCal", Form("calculated Ex [gated]; Ex [MeV]; count / %.2f keV", (ExRange[1]-ExRange[0])/400.*1000), 400, ExRange[0], ExRange[1]); + tree->Draw("ExCal>>hExCal", gate, ""); + Isotope hRecoil(recoil.heavyA, recoil.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}"); + + } + + if( pID == pRecoilRThetaCM ){ + TH2F * hRecoilRThetaCM = new TH2F("hRecoilRThetaCM", "RecoilR - thetaCM [gated]; thetaCM [deg]; RecoilR [mm]", 400, 0, 60, 400,0, detGeo.recoilOuterRadius); + tree->Draw("rhoRecoil:thetaCM>>hRecoilRThetaCM", gate, "colz"); + } + + if( pID == pArrayXY ){ + TH2F * hArrayXY = new TH2F("hArrayXY", "Array-XY [gated]; X [mm]; Y [mm]", 400, -array.detPerpDist*1.5, array.detPerpDist*1.5, 400, -array.detPerpDist*1.5, array.detPerpDist*1.5); + tree->Draw("yArray:xArray>>hArrayXY", gate, "colz"); + } + + if( pID == pInfo ){ + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.06); + text.SetTextColor(2); + + text.DrawLatex(0., 0.9, Reaction); + 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 ) { + vector strList = AnalysisLib::SplitStr( (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)); + } + } + + 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); + tree->Draw("yElum1:xElum1>>hElum1XY", gate, "colz"); + + double count = hElum1XY->GetEntries(); + + if( count < 2000. ) { + hElum1XY->SetMarkerStyle(7); + if( count < 500. ) hElum1XY->SetMarkerStyle(3); + hElum1XY->Draw("scat"); + } + } + + if( pID == pEElum1R ){ + TH2F * hEElum1Rho = new TH2F("hEElum1Rho", "Elum-1 E-R [gated]; R[mm]; Energy[MeV]", 400, 0, elumRange, 400, eRange[0], eRange[1]); + tree->Draw("Tb:rhoElum1>>hEElum1Rho", gate, "colz"); + } + + if( pID == pElum1RThetaCM){ + int angBin = 400; + + TH2F * hElum1RThetaCM = new TH2F("hElum1RThetaCM", "Elum-1 rho vs ThetaCM [gated]; thatCM [deg]; Elum- rho [mm]", angBin, thetaCMRange[0], thetaCMRange[1], 400, 0, elumRange); + tree->Draw("rhoElum1:thetaCM>>hElum1RThetaCM", gate, "colz"); + + TH1F * htemp = (TH1F *) hElum1RThetaCM->ProjectionX("htemp"); + + double rel = (thetaCMRange[1] - thetaCMRange[0])*1.0/angBin; + printf("angular resolution : %f deg \n", rel); + + vector xList; + double old_y = 0; + for( int i = 1; i <= angBin; i++){ + double y = htemp->GetBinContent(i); + if( old_y == 0 && y > 0) xList.push_back(htemp->GetBinCenter(i)); + if( old_y > 0 && y == 0 ) xList.push_back(htemp->GetBinCenter(i)); + old_y = y; + } + + printf("list of gaps :\n"); + for( int i = 0; i < (int) xList.size(); i+=2){ + printf("%d | %.3f - %.3f deg\n", i, xList[i], xList[i+1]); } - if( pID == pElum1XY ){ - - TH2F * hElum1XY = new TH2F("hElum1XY", Form("Elum-1 XY [gated] @ %.0f mm ; X [mm]; Y [mm]", elum1), 400, -elumRange, elumRange, 400, -elumRange, elumRange); - tree->Draw("yElum1:xElum1>>hElum1XY", gate, "colz"); - - double count = hElum1XY->GetEntries(); - - if( count < 2000. ) { - hElum1XY->SetMarkerStyle(7); - if( count < 500. ) hElum1XY->SetMarkerStyle(3); - hElum1XY->Draw("scat"); - } + TF1 f1("f1", "sin(x)"); + double acceptance = 0; + double err1 = 0; + double err2 = 0; + for( int i = 0; i < (int) xList.size(); i += 2 ){ + acceptance += f1.Integral(xList[i] * TMath::DegToRad(), xList[i+1] * TMath::DegToRad() ) * TMath::TwoPi(); + err1 += f1.Integral((xList[i]-rel) * TMath::DegToRad(), (xList[i+1] + rel) * TMath::DegToRad() ) * TMath::TwoPi(); + err2 += f1.Integral((xList[i]+rel) * TMath::DegToRad(), (xList[i+1] - rel) * TMath::DegToRad() ) * TMath::TwoPi(); } + printf("acceptance = %f sr +- %f \n", acceptance, (err1-err2)/2); + + TLatex text; + text.SetTextFont(82); + text.SetTextSize(0.06); + text.SetTextColor(2); + text.SetTextAngle(90); - if( pID == pEElum1R ){ - TH2F * hEElum1Rho = new TH2F("hEElum1Rho", "Elum-1 E-R [gated]; R[mm]; Energy[MeV]", 400, 0, elumRange, 400, eRange[0], eRange[1]); - tree->Draw("Tb:rhoElum1>>hEElum1Rho", gate, "colz"); + for( int i = 0; i < (int) xList.size(); i++){ + text.DrawLatex(xList[i], elumRange/2, Form("%.2f", xList[i])); } - if( pID == pElum1RThetaCM){ - int angBin = 400; - - TH2F * hElum1RThetaCM = new TH2F("hElum1RThetaCM", "Elum-1 rho vs ThetaCM [gated]; thatCM [deg]; Elum- rho [mm]", angBin, thetaCMRange[0], thetaCMRange[1], 400, 0, elumRange); - tree->Draw("rhoElum1:thetaCM>>hElum1RThetaCM", gate, "colz"); - - TH1F * htemp = (TH1F *) hElum1RThetaCM->ProjectionX("htemp"); - - double rel = (thetaCMRange[1] - thetaCMRange[0])*1.0/angBin; - printf("angular resolution : %f deg \n", rel); - - vector xList; - double old_y = 0; - for( int i = 1; i <= angBin; i++){ - double y = htemp->GetBinContent(i); - if( old_y == 0 && y > 0) xList.push_back(htemp->GetBinCenter(i)); - if( old_y > 0 && y == 0 ) xList.push_back(htemp->GetBinCenter(i)); - old_y = y; - } - - printf("list of gaps :\n"); - for( int i = 0; i < (int) xList.size(); i+=2){ - printf("%d | %.3f - %.3f deg\n", i, xList[i], xList[i+1]); - } - - TF1 f1("f1", "sin(x)"); - double acceptance = 0; - double err1 = 0; - double err2 = 0; - for( int i = 0; i < (int) xList.size(); i += 2 ){ - acceptance += f1.Integral(xList[i] * TMath::DegToRad(), xList[i+1] * TMath::DegToRad() ) * TMath::TwoPi(); - err1 += f1.Integral((xList[i]-rel) * TMath::DegToRad(), (xList[i+1] + rel) * TMath::DegToRad() ) * TMath::TwoPi(); - err2 += f1.Integral((xList[i]+rel) * TMath::DegToRad(), (xList[i+1] - rel) * TMath::DegToRad() ) * TMath::TwoPi(); - } - printf("acceptance = %f sr +- %f \n", acceptance, (err1-err2)/2); - - TLatex text; - text.SetTextFont(82); - text.SetTextSize(0.06); - text.SetTextColor(2); - text.SetTextAngle(90); + text.SetNDC(); + text.SetTextAngle(0); + text.DrawLatex(0.15, 0.15, Form("accp. = %.2f(%.2f) msr", acceptance * 1000., (err1-err2)*1000./2)); - for( int i = 0; i < (int) xList.size(); i++){ - text.DrawLatex(xList[i], elumRange/2, Form("%.2f", xList[i])); - } - - text.SetNDC(); - text.SetTextAngle(0); - text.DrawLatex(0.15, 0.15, Form("accp. = %.2f(%.2f) msr", acceptance * 1000., (err1-err2)*1000./2)); - - } + } - if( pID == pHitID ){ - printf("=======================meaning of Hit ID\n"); - printf(" 1 = light recoil hit array & heavy recoil hit recoil\n"); - printf(" 0 = no detector\n"); - printf(" -1 = light recoil go opposite side of array\n"); - printf(" -2 = light recoil hit > det width\n"); - printf(" -3 = light recoil hit > array \n"); - printf(" -4 = light recoil hit blocker \n"); - printf(" -10 = light recoil orbit radius too big \n"); - printf(" -11 = light recoil orbit radius too small\n"); - printf(" -12 = when reocol at the same side of array, light recoil blocked by recoil detector\n"); - printf(" -13 = more than 3 loops\n"); - printf(" -14 = heavy recoil did not hit recoil \n"); - printf(" -15 = cannot find hit on array\n"); - printf(" -20 = unknown\n"); - printf("===========================================\n"); - - TH1F * hHit = new TH1F("hHit", "hit; hit-ID; count", 13, -11, 2); - tree->Draw("hit>>hHit", "", ""); - } - ///####################################################### + if( pID == pHitID ){ + printf("=======================meaning of Hit ID\n"); + printf(" 1 = light recoil hit array & heavy recoil hit recoil\n"); + printf(" 0 = no detector\n"); + printf(" -1 = light recoil go opposite side of array\n"); + printf(" -2 = light recoil hit > det width\n"); + printf(" -3 = light recoil hit > array \n"); + printf(" -4 = light recoil hit blocker \n"); + printf(" -10 = light recoil orbit radius too big \n"); + printf(" -11 = light recoil orbit radius too small\n"); + printf(" -12 = when reocol at the same side of array, light recoil blocked by recoil detector\n"); + printf(" -13 = more than 3 loops\n"); + printf(" -14 = heavy recoil did not hit recoil \n"); + printf(" -15 = cannot find hit on array\n"); + printf(" -20 = unknown\n"); + printf("===========================================\n"); + + TH1F * hHit = new TH1F("hHit", "hit; hit-ID; count", 13, -11, 2); + tree->Draw("hit>>hHit", "", ""); + } - } + } - cCheck->Modified(); - cCheck->Update(); - - if( outputCanvas ){ - TDatime dateTime; - TString outPNGName = Form("Sim_%d%02d%02d_%06d.png", dateTime.GetYear(), dateTime.GetMonth(), dateTime.GetDay(), dateTime.GetTime()); - - cCheck->SaveAs(outPNGName); - printf("%s\n", outPNGName.Data()); - - gROOT->ProcessLine(".q"); - - } + cCheck->Modified(); + cCheck->Update(); + + if( outputCanvas ){ + TDatime dateTime; + TString outPNGName = Form("Sim_%d%02d%02d_%06d.png", dateTime.GetYear(), dateTime.GetMonth(), dateTime.GetDay(), dateTime.GetTime()); + + cCheck->SaveAs(outPNGName); + printf("%s\n", outPNGName.Data()); + + gROOT->ProcessLine(".q"); + + } } ///============================================================ ///============================================================ double * FindRange(TString branch, TString gate, TTree * tree, double output[2]){ - tree->Draw(Form("%s>>temp1", branch.Data()), gate); - TH1F * temp1 = (TH1F *) gROOT->FindObjectAny("temp1"); + tree->Draw(Form("%s>>temp1", branch.Data()), gate); + TH1F * temp1 = (TH1F *) gROOT->FindObjectAny("temp1"); - output[1] = temp1->GetXaxis()->GetXmax(); - output[0] = temp1->GetXaxis()->GetXmin(); + output[1] = temp1->GetXaxis()->GetXmax(); + output[0] = temp1->GetXaxis()->GetXmin(); - delete temp1; - return output; + delete temp1; + return output; } double ExtractNumber(int index, TMacro * macro){ @@ -584,75 +571,75 @@ TString ExtractString(int index, TMacro * macro){ vector StringToVector(TString str){ - vector temp; - - bool startFlag = false; - bool endFlag = false; - string jaja=""; - for(int i = 0; i < str.Length(); i++){ - - if( str[i] == '{' ) { - startFlag = true; - continue; - } - if( str[i] == ' '){ - continue; - } - if( startFlag && !endFlag){ - - if( str[i] == ',' ){ - temp.push_back(jaja); - jaja=""; - continue; - } - - if( str[i] == '}') { - temp.push_back(jaja); - endFlag = true; - continue; - } - jaja += str[i]; - - } - } - return temp; + vector temp; + + bool startFlag = false; + bool endFlag = false; + string jaja=""; + for(int i = 0; i < str.Length(); i++){ + + if( str[i] == '{' ) { + startFlag = true; + continue; + } + if( str[i] == ' '){ + continue; + } + if( startFlag && !endFlag){ + + if( str[i] == ',' ){ + temp.push_back(jaja); + jaja=""; + continue; + } + + if( str[i] == '}') { + temp.push_back(jaja); + endFlag = true; + continue; + } + jaja += str[i]; + + } + } + return temp; } vector intConvertor(vector arr){ - vector out ; - for( int i = 0 ; i < (int) arr.size(); i++){ - out.push_back( arr[i].Atoi()); - } - return out; + vector out ; + for( int i = 0 ; i < (int) arr.size(); i++){ + out.push_back( arr[i].Atoi()); + } + return out; } vector doubleConvertor(vector arr){ - vector out ; - for( int i = 0 ; i < (int) arr.size(); i++){ - out.push_back( arr[i].Atof()); - } - return out; + vector out ; + for( int i = 0 ; i < (int) arr.size(); i++){ + out.push_back( arr[i].Atof()); + } + return out; } plotID StringToPlotID(TString str){ - if( str == "pEZ") return plotID::pEZ; ///0 - if( str == "pRecoilXY") return plotID::pRecoilXY; /// 1 - if( str == "pRecoilXY1" ) return plotID::pRecoilXY1; /// 2 - if( str == "pRecoilXY2" ) return plotID::pRecoilXY2; /// 3 - if( str == "pRecoilRZ" ) return plotID::pRecoilRZ; /// 4 - if( str == "pRecoilRTR" ) return plotID::pRecoilRTR; /// 5 - if( str == "pTDiffZ" ) return plotID::pTDiffZ; /// 6 - if( str == "pThetaCM" ) return plotID::pThetaCM; /// 7 - if( str == "pThetaCM_Z" ) return plotID::pThetaCM_Z; /// 8 - if( str == "pExCal" ) return plotID::pExCal; /// 9 - if( str == "pRecoilRThetaCM" ) return plotID::pRecoilRThetaCM; /// 10 - if( str == "pArrayXY" ) return plotID::pArrayXY; /// 11 - if( str == "pInfo" ) return plotID::pInfo; /// 12 - if( str == "pHitID" ) return plotID::pHitID; /// 13 - if( str == "pElum1XY" ) return plotID::pElum1XY; /// 14 - if( str == "pEElum1R" ) return plotID::pEElum1R; /// 14 - if( str == "pElum1RThetaCM" ) return plotID::pElum1RThetaCM; /// 15 - if( str == "pEmpty" ) return plotID::pEmpty ; /// 16 - - return plotID::pEmpty; + if( str == "pEZ") return plotID::pEZ; ///0 + if( str == "pRecoilXY") return plotID::pRecoilXY; /// 1 + if( str == "pRecoilXY1" ) return plotID::pRecoilXY1; /// 2 + if( str == "pRecoilXY2" ) return plotID::pRecoilXY2; /// 3 + if( str == "pRecoilRZ" ) return plotID::pRecoilRZ; /// 4 + if( str == "pRecoilRTR" ) return plotID::pRecoilRTR; /// 5 + if( str == "pTDiffZ" ) return plotID::pTDiffZ; /// 6 + if( str == "pThetaCM" ) return plotID::pThetaCM; /// 7 + if( str == "pThetaCM_Z" ) return plotID::pThetaCM_Z; /// 8 + if( str == "pExCal" ) return plotID::pExCal; /// 9 + if( str == "pRecoilRThetaCM" ) return plotID::pRecoilRThetaCM; /// 10 + if( str == "pArrayXY" ) return plotID::pArrayXY; /// 11 + if( str == "pInfo" ) return plotID::pInfo; /// 12 + if( str == "pHitID" ) return plotID::pHitID; /// 13 + if( str == "pElum1XY" ) return plotID::pElum1XY; /// 14 + if( str == "pEElum1R" ) return plotID::pEElum1R; /// 14 + if( str == "pElum1RThetaCM" ) return plotID::pElum1RThetaCM; /// 15 + if( str == "pEmpty" ) return plotID::pEmpty ; /// 16 + + return plotID::pEmpty; } diff --git a/Cleopatra/alpha.C b/Cleopatra/SimAlpha.C similarity index 100% rename from Cleopatra/alpha.C rename to Cleopatra/SimAlpha.C diff --git a/Cleopatra/knockoutSim.C b/Cleopatra/SimKnockout.C similarity index 100% rename from Cleopatra/knockoutSim.C rename to Cleopatra/SimKnockout.C diff --git a/Cleopatra/Transfer.h b/Cleopatra/SimTransfer.C similarity index 89% rename from Cleopatra/Transfer.h rename to Cleopatra/SimTransfer.C index ef98d46..0274076 100644 --- a/Cleopatra/Transfer.h +++ b/Cleopatra/SimTransfer.C @@ -45,8 +45,7 @@ void Transfer( TString ptolemyRoot = "DWBA.root", TString saveFileName = "transfer.root"){ - //############################################# Set Reaction - + //*############################################# Set Reaction TransferReaction transfer; HELIOS helios; Decay decay; @@ -65,7 +64,7 @@ void Transfer( printf("\e[32m#################################### Reaction & HELIOS configuration\e[0m\n"); - transfer.PrintReaction(); + transfer.PrintReaction(false); if(transfer.GetRecoil().isDecay) { decay.SetMotherDaugther(transfer.GetRecoil()); @@ -118,19 +117,9 @@ void Transfer( // msB.LoadStoppingPower(reactConfig.recoilHeavyStoppingPowerFile); // } - //*############################################# Decay of particle-B - // Decay decay[2]; - // if(reactConfig.isDecay) { - // printf("\e[32m#################################### Decay\e[0m\n"); - // decay.SetMotherDaugther(reactConfig.recoilHeavyA, - // reactConfig.recoilHeavyZ, - // reactConfig.heavyDecayA, - // reactConfig.heavyDecayZ); - // } - ExcitedEnergies exList = transfer.GetRectionConfig().exList[ID]; - //############################################# Load DWBAroot for thetaCM distribution + //*############################################# Load DWBAroot for thetaCM distribution printf("\e[32m#################################### Load DWBA input : %s \e[0m\n", ptolemyRoot.Data()); TF1 * dist = NULL; TFile * distFile = new TFile(ptolemyRoot, "read"); @@ -178,7 +167,7 @@ void Transfer( } } - //############################################# build tree + //*############################################# build tree printf("\e[32m#################################### building Tree in %s\e[0m\n", saveFileName.Data()); TFile * saveFile = new TFile(saveFileName, "recreate"); TTree * tree = new TTree("tree", "tree"); @@ -192,6 +181,9 @@ void Transfer( if( distList != NULL ) distList->Write("DWBA", 1); if( dwbaExList != NULL ) dwbaExList->Write("DWBA_ExList", 1); + TMacro idMacro; + idMacro.AddLine(Form("%d", ID)); + idMacro.Write("detGeoID"); TMacro hitMeaning; hitMeaning.AddLine("======================= meaning of Hit\n"); @@ -421,7 +413,7 @@ void Transfer( double theta = reactConfig.beamTheta; double phi = 0.0; - //====================================================== calculate event + //*====================================================== calculate event int count = 0; for( int i = 0; i < numEvent; i++){ bool redoFlag = true; @@ -497,7 +489,6 @@ void Transfer( int decayID = 0; if( recoil.isDecay){ - //decayID = decay.CalDecay(PB, Ex, 0, phiCM + TMath::Pi()/2.); // decay to ground state decayID = decay.CalDecay(PB, Ex, 0, phiCM + TMath::Pi()/2); // decay to ground state if( decayID == 1 ){ PB = decay.GetDaugther_D(); @@ -661,3 +652,58 @@ void Transfer( return; } + + +int main (int argc, char *argv[]) { + + printf("=================================================================\n"); + printf("========== Simulate Transfer reaction in HELIOS ==========\n"); + printf("=================================================================\n"); + + if(argc == 2 || argc > 7) { + printf("Usage: ./Transfer [1] [2] [3] [4] [5] [6]\n"); + printf(" default file name \n"); + printf(" [1] reactionConfig.txt (input) reaction Setting \n"); + printf(" [2] detectorGeo.txt (input) detector Setting \n"); + printf(" [3] ID (input) detector & reaction ID (default = 0 ) \n"); + printf(" [4] DWBA.root (input) thetaCM distribution from DWBA \n"); + printf(" [5] transfer.root (output) rootFile name for output \n"); + printf(" [6] plot (input) will it plot stuffs [1/0] \n"); + + printf("------------------------------------------------------\n"); + return 0 ; + } + + std::string basicConfig = "reactionConfig.txt"; + std::string heliosDetGeoFile = "detectorGeo.txt"; + int ID = 0; + TString ptolemyRoot = "DWBA.root"; // when no file, use isotropic distribution of thetaCM + TString saveFileName = "transfer.root"; + bool isPlot = false; + + if( argc >= 2) basicConfig = argv[1]; + if( argc >= 3) heliosDetGeoFile = argv[2]; + if( argc >= 4) ID = atoi(argv[3]); + if( argc >= 5) ptolemyRoot = argv[4]; + if( argc >= 6) saveFileName = argv[5]; + if( argc >= 7) isPlot = atoi(argv[7]); + + Transfer( basicConfig, heliosDetGeoFile, ID, ptolemyRoot, saveFileName); + + //run Armory/Check_Simulation + if( isPlot ){ + ifstream file_in; + file_in.open("../Cleopatra/Check_Simulation.C", ios::in); + if( file_in){ + printf("---- running ../Cleopatra/Check_Simulation.C on %s \n", saveFileName.Data()); + TString cmd; + cmd.Form("root -l '../Cleopatra/Check_Simulation.C(\"%s\")'", saveFileName.Data()); + + system(cmd.Data()); + }else{ + printf("cannot find ../Cleopatra/Check_Simulation.C \n"); + } + } + +} + diff --git a/Cleopatra/Simulation_Helper.C b/Cleopatra/Simulation_Helper.C index 9f306f2..c12b75c 100644 --- a/Cleopatra/Simulation_Helper.C +++ b/Cleopatra/Simulation_Helper.C @@ -12,7 +12,7 @@ #include -#include "../Cleopatra/Transfer.h" +#include "../Cleopatra/SimTransfer.C" #include "../Cleopatra/InFileCreator.h" #include "../Cleopatra/ExtractXSec.h" #include "../Cleopatra/PlotTGraphTObjArray.h" @@ -58,7 +58,6 @@ private: TGComboBox * extractFlag; - TGTextEntry * txtName ; TGTextEntry * txtEx ; diff --git a/Cleopatra/Transfer.C b/Cleopatra/Transfer.C deleted file mode 100644 index 3f851e3..0000000 --- a/Cleopatra/Transfer.C +++ /dev/null @@ -1,72 +0,0 @@ -/*********************************************************************** - * - * This is Transfer.C for simulation of transfer reaction. - * - * ----------------------------------------------------- - * This program will call the root library and compile in g++ - * compilation: - * g++ Transfer.C -o Transfer `root-config --cflags --glibs` - * - * ------------------------------------------------------ - * created by Ryan (Tsz Leung) Tang, Feb-4, 2019 - * email: goluckyryan@gmail.com - * ********************************************************************/ - -#include -#include -#include "Transfer.h" - -using namespace std; - -int main (int argc, char *argv[]) { - - printf("=================================================================\n"); - printf("========== Simulate Transfer reaction in HELIOS ==========\n"); - printf("=================================================================\n"); - - if(argc == 2 || argc > 7) { - printf("Usage: ./Transfer [1] [2] [3] [4] [5] [6]\n"); - printf(" default file name \n"); - printf(" [1] reactionConfig.txt (input) reaction Setting \n"); - printf(" [2] detectorGeo.txt (input) detector Setting \n"); - printf(" [3] ID (input) detector & reaction ID (default = 0 ) \n"); - printf(" [4] DWBA.root (input) thetaCM distribution from DWBA \n"); - printf(" [5] transfer.root (output) rootFile name for output \n"); - printf(" [6] plot (input) will it plot stuffs [1/0] \n"); - - printf("------------------------------------------------------\n"); - return 0 ; - } - - string basicConfig = "reactionConfig.txt"; - string heliosDetGeoFile = "detectorGeo.txt"; - int ID = 0; - TString ptolemyRoot = "DWBA.root"; // when no file, use isotropic distribution of thetaCM - TString saveFileName = "transfer.root"; - bool isPlot = false; - - if( argc >= 2) basicConfig = argv[1]; - if( argc >= 3) heliosDetGeoFile = argv[2]; - if( argc >= 4) ID = atoi(argv[3]); - if( argc >= 5) ptolemyRoot = argv[4]; - if( argc >= 6) saveFileName = argv[5]; - if( argc >= 7) isPlot = atoi(argv[7]); - - Transfer( basicConfig, heliosDetGeoFile, ID, ptolemyRoot, saveFileName); - - //run Armory/Check_Simulation - if( isPlot ){ - ifstream file_in; - file_in.open("../Cleopatra/Check_Simulation.C", ios::in); - if( file_in){ - printf("---- running ../Cleopatra/Check_Simulation.C on %s \n", saveFileName.Data()); - TString cmd; - cmd.Form("root -l '../Cleopatra/Check_Simulation.C(\"%s\")'", saveFileName.Data()); - - system(cmd.Data()); - }else{ - printf("cannot find ../Cleopatra/Check_Simulation.C \n"); - } - } - -} diff --git a/Cleopatra/makefile b/Cleopatra/makefile index 6e0bc3b..f74a9cb 100644 --- a/Cleopatra/makefile +++ b/Cleopatra/makefile @@ -1,6 +1,6 @@ CC=g++ -ALL = Isotope InFileCreator ExtractXSec ExtractXSecFromText PlotTGraphTObjArray Cleopatra FindThetaCM Transfer SimAlpha +ALL = Isotope InFileCreator ExtractXSec ExtractXSecFromText PlotTGraphTObjArray Cleopatra FindThetaCM SimTransfer SimAlpha all: $(ALL) @@ -25,11 +25,11 @@ Cleopatra: Cleopatra.C FindThetaCM: FindThetaCM.C FindThetaCM.h ../Cleopatra/ClassTransfer.h ../Cleopatra/ClassHelios.h ../Cleopatra/ClassIsotope.h ../Cleopatra/constant.h $(CC) FindThetaCM.C -o FindThetaCM `root-config --cflags --glibs` -Transfer: Transfer.C Transfer.h ../Cleopatra/ClassTransfer.h ../Cleopatra/ClassHelios.h ../Cleopatra/ClassIsotope.h ../Cleopatra/constant.h - $(CC) Transfer.C -o Transfer `root-config --cflags --glibs` +SimTransfer: SimTransfer.C ../Cleopatra/ClassTransfer.h ../Cleopatra/ClassHelios.h ../Cleopatra/ClassIsotope.h ../Cleopatra/constant.h + $(CC) SimTransfer.C -o SimTransfer `root-config --cflags --glibs` -SimAlpha: alpha.C ../Cleopatra/ClassHelios.h - $(CC) alpha.C -o SimAlpha `root-config --cflags --glibs` +SimAlpha: SimAlpha.C ../Cleopatra/ClassHelios.h + $(CC) SimAlpha.C -o SimAlpha `root-config --cflags --glibs` clean: /bin/rm -f $(ALL) \ No newline at end of file