diff --git a/.vscode/settings.json b/.vscode/settings.json index 9d8c18b..4fa9342 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -3,7 +3,7 @@ "detectorGeo.txt": "shellscript", "reactionConfig.txt": "shellscript", "DWBA": "shellscript", - "SimChecker_Config.txt": "cpp", + "SimChecker_Config.txt": "fsharp", "*.C": "cpp" }, diff --git a/Cleopatra/Check_Simulation.C b/Cleopatra/Check_Simulation_obsolete.C similarity index 99% rename from Cleopatra/Check_Simulation.C rename to Cleopatra/Check_Simulation_obsolete.C index 7221cff..65b2ad2 100644 --- a/Cleopatra/Check_Simulation.C +++ b/Cleopatra/Check_Simulation_obsolete.C @@ -252,7 +252,6 @@ void Check_Simulation(TString filename = "transfer.root", 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); @@ -633,23 +632,24 @@ vector doubleConvertor(vector arr){ 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 == "pEZ") return plotID::pEZ; /// 0 + if( str == "pRecoilXY") return plotID::pRecoilXY; /// 1 + if( str == "pThetaCM" ) return plotID::pThetaCM; /// 2 + if( str == "pExCal" ) return plotID::pExCal; /// 2 + if( str == "pArrayXY" ) return plotID::pArrayXY; /// 3 + if( str == "pInfo" ) return plotID::pInfo; /// 4 + if( str == "pElum1XY" ) return plotID::pElum1XY; /// 5 + if( str == "pRecoilXY1" ) return plotID::pRecoilXY1; /// 6 + if( str == "pRecoilXY2" ) return plotID::pRecoilXY2; /// 7 + if( str == "pTDiffZ" ) return plotID::pTDiffZ; /// 8 + if( str == "pRecoilRThetaCM" ) return plotID::pRecoilRThetaCM; /// 9 + if( str == "pRecoilRZ" ) return plotID::pRecoilRZ; /// 10 + if( str == "pEElum1R" ) return plotID::pEElum1R; /// 11 + if( str == "pRecoilRTR" ) return plotID::pRecoilRTR; /// 12 + if( str == "pThetaCM_Z" ) return plotID::pThetaCM_Z; /// 13 + if( str == "pElum1RThetaCM" ) return plotID::pElum1RThetaCM; /// 14 + 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/ClassSimPlotter.h b/Cleopatra/ClassSimPlotter.h new file mode 100644 index 0000000..453e009 --- /dev/null +++ b/Cleopatra/ClassSimPlotter.h @@ -0,0 +1,501 @@ +#ifndef SimPlotter_h +#define SimPlotter_h + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "../Armory/ClassDetGeo.h" +#include "../Armory/ClassReactionConfig.h" +#include "../Cleopatra/ClassIsotope.h" + +enum plotID { pEZ, /// 0 + pRecoilXY, /// 1 + pRecoilXY1, /// 2 + pRecoilXY2, /// 3 + pRecoilRZ, /// 4 + pRecoilRTR, /// 5 + pTDiffZ, /// 6 + pThetaCM, /// 7 + pThetaCM_Z, /// 8 + pExCal, /// 9 + pRecoilRThetaCM, /// 10 + pArrayXY, /// 11 + pInfo, /// 12 + pElum1XY, /// 13 + pEElum1R, /// 14 + pElum1RThetaCM, /// 15 + pEmpty }; /// 16 + +class Plotter{ +public: + Plotter(TTree * tree, short reactionID, std::string reactionName, DetGeo detGeo, ExcitedEnergies exList, TString gate, std::vector padPlotID); + ~Plotter(); + + void SetRanges(double zMin, double zMax, double eMax, double elumMax, double thetaCM_Max){ + this->zRange[0] = zMin; + this->zRange[1] = zMax; + this->eMax = eMax; + this->elumMax = elumMax; + this->thetaCM_Max = thetaCM_Max; + } + void SetCanvas( int colSize, int rowSize , int padSize, std::vector padPlotID ){ + cCheck = new TCanvas(Form("SimChecker-%d", rID), + "Simulation Checker", + 200 * rID, + 200 * rID, + padSize*colSize, + padSize*rowSize); + + if(cCheck->GetShowEditor() ) cCheck->ToggleEditor(); + if(cCheck->GetShowToolBar() )cCheck->ToggleToolBar(); + cCheck->Divide(colSize,rowSize); + + numPad = colSize * rowSize; + } + + void Plot(); + +private: + + TTree * tree; + DetGeo detGeo; + std::string reactionName; + short rID; + ExcitedEnergies exList; + int numEx; + TString gate; + TString gateTrue; + std::vector padPlotID; + + double zRange[2]; + double eMax; + double recoilOut; + double recoilIn; + double elumMax; + double thetaCM_Max; + + int numPad; + TCanvas * cCheck; + TH2F * hez; + TH2F * hArrayXY; + TH2F * heavyXY; + TH2F * hElum1XY; + TH1F ** hThetaCM; + TH1F * hExCal; + TH1F * hHit; + TH2F * hRecoilXY1; + TH2F * hRecoilXY2; + TH2F * hTDiffZ; + TH2F * hRecoilRThetaCM; + TH2F * hRecoilRZ; + TH2F * hEElum1Rho; + TH2F * hRecoilRTR; + TH2F * hElum1RThetaCM; + TH2F * hThetaCM_Z; + + Isotope heavy; + + std::vector FindRange(TString branch); + +}; + +Plotter::Plotter(TTree * tree, short reactionID, std::string reactionName, DetGeo detGeo, ExcitedEnergies exList, TString gate, std::vector padPlotID){ + this->tree = tree; + this->rID = reactionID; + this->reactionName = reactionName; + this->exList = exList; + numEx = this->exList.ExList.size(); + + this->detGeo = detGeo; + + this->recoilOut = detGeo.aux[rID].outerRadius; + this->recoilIn = detGeo.aux[rID].innerRadius; + + this->gate = gate; + this->gateTrue = gate + Form(" && rID == %d", rID); + + this->padPlotID = padPlotID; + + cCheck = nullptr; + hez = nullptr; + hArrayXY = nullptr; + heavyXY = nullptr; + hElum1XY = nullptr; + hThetaCM = nullptr; + hExCal = nullptr; + hHit = nullptr; + hRecoilXY1 = nullptr; + hRecoilXY2 = nullptr; + hTDiffZ = nullptr; + hRecoilRThetaCM = nullptr; + hRecoilRZ = nullptr; + hEElum1Rho = nullptr; + hRecoilRTR = nullptr; + hElum1RThetaCM = nullptr; + hThetaCM_Z = nullptr; + + std::size_t pos1 = reactionName.find(")") + 1; + std::size_t pos2 = reactionName.find(" ", pos1); + std::string result = reactionName.substr(pos1, pos2 - pos1); + + pos1 = result.find("{") + 1; + pos2 = result.find("}", pos1); + + std::string massHeavy = result.substr(pos1, pos2-pos1); + std::string symHeavy = result.substr(pos2+1); + + heavy.SetIsoByName(massHeavy + symHeavy); + +} + +Plotter::~Plotter(){ + + if( cCheck ) delete cCheck; + if( hez ) delete hez; + if( hArrayXY ) delete hArrayXY; + if( heavyXY ) delete heavyXY; + if( hElum1XY ) delete hElum1XY; + if( hExCal ) delete hExCal; + if( hRecoilXY1 ) delete hRecoilXY1; + if( hRecoilXY2 ) delete hRecoilXY2; + if( hTDiffZ ) delete hTDiffZ; + if( hRecoilRThetaCM ) delete hRecoilRThetaCM; + if( hRecoilRZ ) delete hRecoilRZ; + if( hEElum1Rho ) delete hEElum1Rho; + if( hRecoilRTR ) delete hRecoilRTR; + if( hElum1RThetaCM ) delete hElum1RThetaCM; + if( hThetaCM_Z ) delete hThetaCM_Z; + + if(hThetaCM ){ + for( int i = 0; i < numEx; i++) delete [] hThetaCM[i]; + delete [] hThetaCM; + } + +} + +inline void Plotter::Plot(){ + + for( int i = 1; i <= numPad; i++ ){ + cCheck->cd(i); + plotID pID = padPlotID[i-1]; + + //&======================================= 0 + if( pID == pEZ){ + hez = new TH2F(Form("hez%d",rID), //need unqie name in root. + Form("e-z [gated] @ %5.0f mm; z [mm]; e [MeV]", detGeo.array[rID].firstPos), + 400, zRange[0], zRange[1], + 400, 0, eMax); + + tree->Draw(Form("e:z>>hez%d", rID), gateTrue , "colz"); + } + + //&======================================= 1 + if( pID == pRecoilXY ){ + heavyXY = new TH2F(Form("heavyXY%d",rID), + Form("RecoilXY [gated] @ %4.0f mm; X [mm]; Y [mm]", detGeo.aux[rID].detPos ), + 400, -recoilOut, recoilOut, + 400, -recoilOut, recoilOut); + tree->Draw(Form("yRecoil:xRecoil>>heavyXY%d", rID), gateTrue, "colz"); + + // printf("%f, %f\n", recoilOut, recoilIn); + TArc * detArc1 = new TArc(0, 0, recoilOut); + detArc1->SetLineColor(kBlue-8); + detArc1->SetFillStyle(0); + detArc1->Draw("same"); + TArc * detArc2 = new TArc(0, 0, recoilIn); + 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"); + // } + } + + //&======================================= 2 + if( pID == pThetaCM ){ + hThetaCM = new TH1F *[numEx]; + TLegend * legend = new TLegend(0.8,0.2,0.99,0.8); + double maxCount = 0; + for( int h = 0; h < numEx; h++){ + hThetaCM[h] = new TH1F(Form("hThetaCM%d-%d", rID, h), Form("thetaCM [gated] (ExID=%d); thetaCM [deg]; count", h), 200, 0, thetaCM_Max); + hThetaCM[h]->SetLineColor(h+1); + hThetaCM[h]->SetFillColor(h+1); + hThetaCM[h]->SetFillStyle(3000+h); + tree->Draw(Form("thetaCM>>hThetaCM%d-%d", rID, h), gateTrue + Form("&& ExID==%d", h), ""); + legend->AddEntry(hThetaCM[h], Form("Ex=%5.1f MeV", exList.ExList[h].Ex)); + double max = hThetaCM[h]->GetMaximum(); + if( max > maxCount ) maxCount = max; + } + + for( int h = 0; h < numEx; h++){ + hThetaCM[h]->GetYaxis()->SetRangeUser(1, maxCount * 1.2); + if( h == 0 ) { + hThetaCM[h]->Draw(); + }else{ + hThetaCM[h]->Draw("same"); + } + } + legend->Draw(); + } + + //&======================================= 3 + if( pID == pExCal ){ + double exMin = 9999; + double exMax = -1; + for( size_t k = 0 ; k < exList.ExList.size(); k++){ + double kuku = exList.ExList[k].Ex; + if( kuku > exMax ) exMax = kuku; + if( kuku < exMin ) exMin = kuku; + } + + double exPadding = (exMax - exMin) < 1 ? 1 : (exMax - exMin) * 0.3; + + exMin = exMin - exPadding ; + exMax = exMax + exPadding ; + + hExCal = new TH1F(Form("hExCal%d",rID), + Form("calculated Ex [gated]; Ex [MeV]; count / %.2f keV", (exMax-exMin)/400.*1000), + 400, exMin, exMax); + tree->Draw(Form("ExCal>>hExCal%d", rID), gateTrue, ""); + + + double Sn = heavy.CalSp(0,1); + double Sp = heavy.CalSp(1,0); + double Sa = heavy.CalSp2(4,2); + double S2n = heavy.CalSp(0, 2); + + printf("======= Heavy recoil: %s \n", heavy.Name.c_str()); + printf("Sn : %6.2f MeV/u \n", Sn); + printf("Sp : %6.2f MeV/u \n", Sp); + printf("Sa : %6.2f MeV/u \n", Sa); + printf("S2n : %6.2f MeV/u \n", S2n); + + double yMax = hExCal->GetMaximum(); + TLatex * text = new TLatex(); + text->SetTextFont(82); + text->SetTextSize(0.06); + + if( Sn < exMax ) {TLine * lineSn = new TLine(Sn, 0, Sn, yMax); lineSn->SetLineColor(2); lineSn->Draw(""); text->SetTextColor(2); text->DrawLatex(Sn, yMax*0.9, "S_{n}");} + if( Sp < exMax ) {TLine * lineSp = new TLine(Sp, 0, Sp, yMax); lineSp->SetLineColor(4); lineSp->Draw("same"); text->SetTextColor(4); text->DrawLatex(Sp, yMax*0.9, "S_{p}");} + if( Sa < exMax ) {TLine * lineSa = new TLine(Sa, 0, Sa, yMax); lineSa->SetLineColor(6); lineSa->Draw("same"); text->SetTextColor(6); text->DrawLatex(Sa, yMax*0.9, "S_{a}");} + if( S2n < exMax ) {TLine * lineS2n = new TLine(S2n, 0, S2n, yMax); lineS2n->SetLineColor(8); lineS2n->Draw("same"); text->SetTextColor(8); text->DrawLatex(S2n, yMax*0.9, "S_{2n}");} + + } + + //&======================================= 4 + if( pID == pArrayXY ){ + + hArrayXY = new TH2F(Form("hArrayXY%d", rID), + "Array-XY [gated]; X [mm]; Y [mm]", + 400, -detGeo.array[rID].detPerpDist*1.5, detGeo.array[rID].detPerpDist*1.5, + 400, -detGeo.array[rID].detPerpDist*1.5, detGeo.array[rID].detPerpDist*1.5); + tree->Draw(Form("yArray:xArray>>hArrayXY%d", rID), gateTrue, "colz"); + } + + //&======================================= 5 + if( pID == pInfo ){ + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.06); + text.SetTextColor(2); + + text.DrawLatex(0., 0.9, reactionName.c_str()); + 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)); + // } + } + + //&======================================= 6 + if( pID == pElum1XY ){ + hElum1XY = new TH2F(Form("hElum1XY%d", rID), + Form("Elum-1 XY [gated] @ %.0f mm ; X [mm]; Y [mm]", detGeo.aux[rID].elumPos1), + 400, -elumMax, elumMax, + 400, -elumMax, elumMax); + tree->Draw(Form("yElum1:xElum1>>hElum1XY%d", rID), gateTrue, "colz"); + + double count = hElum1XY->GetEntries(); + if( count < 2000. ) { + hElum1XY->SetMarkerStyle(7); + if( count < 500. ) hElum1XY->SetMarkerStyle(3); + hElum1XY->Draw("scat"); + } + } + + //&======================================= 7 + if( pID == pRecoilXY1 ){ + TH2F * hRecoilXY1 = new TH2F(Form("hRecoilXY1-%d", rID), + Form("RecoilXY-1 [gated] @ %4.0f mm; X [mm]; Y [mm]", detGeo.aux[rID].detPos1 ), + 400, -detGeo.aux[rID].outerRadius, detGeo.aux[rID].outerRadius, + 400, -detGeo.aux[rID].outerRadius, detGeo.aux[rID].outerRadius); + tree->Draw(Form("yRecoil1:xRecoil1>>hRecoilXY1-%d", rID), gateTrue, "colz"); + } + + //&======================================= 8 + if( pID == pRecoilXY2 ){ + hRecoilXY2 = new TH2F(Form("hRecoilXY2=%d", rID), + Form("RecoilXY-2 [gated] @ %4.0f mm; X [mm]; Y [mm]", detGeo.aux[rID].detPos2 ), + 400, -detGeo.aux[rID].outerRadius, detGeo.aux[rID].outerRadius, + 400, -detGeo.aux[rID].outerRadius, detGeo.aux[rID].outerRadius); + tree->Draw(Form("yRecoil2:xRecoil2>>hRecoilXY2-%d", rID), gate, "colz"); + } + + //&======================================= 9 + if( pID == pTDiffZ ){ + std::vector tDiffRange = FindRange("t-tB"); + hTDiffZ = new TH2F(Form("hTDiffZ%d", rID), + "time(Array) - time(Recoil) vs Z [gated]; z [mm]; time diff [ns]", + 400, zRange[0], zRange[1], + 400, tDiffRange[0], tDiffRange[1]); + tree->Draw(Form("t - tB : z >> hTDiffZ%d", rID), gateTrue, "colz"); + } + + //&======================================= 10 + if( pID == pRecoilRThetaCM ){ + hRecoilRThetaCM = new TH2F(Form("hRecoilRThetaCM%d", rID), + "RecoilR - thetaCM [gated]; thetaCM [deg]; RecoilR [mm]", + 400, 0, 60, + 400, 0, detGeo.aux[rID].outerRadius); + tree->Draw(Form("rhoRecoil:thetaCM>>hRecoilRThetaCM%d", rID), gateTrue, "colz"); + } + + //&======================================= 11 + if( pID == pRecoilRZ ){ + hRecoilRZ = new TH2F(Form("hRecoilRZ%d", rID), + "RecoilR - Z [gated]; z [mm]; RecoilR [mm]", + 400, zRange[0], zRange[1], + 400,0, detGeo.aux[rID].outerRadius); + tree->Draw(Form("rhoRecoil:z>>hRecoilRZ%d", rID), gateTrue, "colz"); + } + + //&======================================= 12 + if( pID == pEElum1R ){ + hEElum1Rho = new TH2F(Form("hEElum1Rho%d", rID), + "Elum-1 E-R [gated]; R[mm]; Energy[MeV]", + 400, 0, elumMax, + 400, 0, eMax); + tree->Draw(Form("Tb:rhoElum1>>hEElum1Rho%d", rID), gateTrue, "colz"); + } + + //&======================================= 13 + if( pID == pRecoilRTR ){ + std::vector recoilERange = FindRange("TB"); + hRecoilRTR = new TH2F(Form("hRecoilRTR%d", rID), + "RecoilR - recoilE [gated]; recoil Energy [MeV]; RecoilR [mm]", + 500, recoilERange[0], recoilERange[1], + 500, 0, detGeo.aux[rID].outerRadius); + tree->Draw(Form("rhoRecoil:TB>>hRecoilRTR%d", rID), gateTrue, "colz"); + } + + //&======================================= 14 + if( pID == pThetaCM_Z ){ + hThetaCM_Z = new TH2F(Form("hThetaCM_Z%d", rID), + "ThetaCM vs Z ; Z [mm]; thetaCM [deg]", + 400, zRange[0], zRange[1], + 400, 0, thetaCM_Max); + tree->Draw(Form("thetaCM:z>>hThetaCM_Z%d", rID), gateTrue,"col"); + } + + //&======================================= 15 + if( pID == pElum1RThetaCM){ + int angBin = 400; + + hElum1RThetaCM = new TH2F(Form("hElum1RThetaCM%d", rID), + "Elum-1 rho vs ThetaCM [gated]; thatCM [deg]; Elum- rho [mm]", + angBin, 0, thetaCM_Max, + 400, 0, elumMax); + tree->Draw(Form("rhoElum1:thetaCM>>hElum1RThetaCM%d", rID), gateTrue, "colz"); + + TH1F * htemp = (TH1F *) hElum1RThetaCM->ProjectionX("htemp"); + + double rel = (thetaCM_Max)*1.0/angBin; + printf("angular resolution : %f deg \n", rel); + + std::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); + + for( int i = 0; i < (int) xList.size(); i++){ + text.DrawLatex(xList[i], elumMax/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)); + + delete htemp; + + } + + }///===== end of pad loop + +} + +std::vector Plotter::FindRange(TString branch){ + + tree->Draw(Form("%s>>temp1", branch.Data()), gateTrue); + TH1F * temp1 = (TH1F *) gROOT->FindObjectAny("temp1"); + + std::vector output; + output.push_back( temp1->GetXaxis()->GetXmax() ); + output.push_back( temp1->GetXaxis()->GetXmin() ); + + delete temp1; + return output; +} + + +#endif \ No newline at end of file diff --git a/Cleopatra/SimChecker.C b/Cleopatra/SimChecker.C index 43e3ad1..43c9914 100644 --- a/Cleopatra/SimChecker.C +++ b/Cleopatra/SimChecker.C @@ -23,25 +23,8 @@ #include "../Armory/ClassReactionConfig.h" #include "../Cleopatra/ClassIsotope.h" #include "../Cleopatra/ClassTransfer.h" +#include "../Cleopatra/ClassSimPlotter.h" -enum plotID { pEZ, /// 0 - pRecoilXY, /// 1 - pRecoilXY1, /// 2 - pRecoilXY2, /// 3 - pRecoilRZ, /// 4 - pRecoilRTR, /// 5 - pTDiffZ, /// 6 - pThetaCM, /// 7 - pThetaCM_Z, /// 8 - pExCal, /// 9 - pRecoilRThetaCM, /// 10 - pArrayXY, /// 11 - pInfo, /// 12 - pHitID, /// 13 - pElum1XY, /// 14 - pEElum1R, /// 15 - pElum1RThetaCM, /// 16 - pEmpty }; /// 17 plotID StringToPlotID(TString str); void SimChecker(TString filename = "transfer.root", @@ -132,7 +115,6 @@ void SimChecker(TString filename = "transfer.root", // } // //^################################################ - TMacro * config = new TMacro(configFile); int numLine = config->GetListOfLines()->GetSize(); @@ -190,8 +172,6 @@ void SimChecker(TString filename = "transfer.root", } } - - int Div[2] = {colCount, rowCount}; gStyle->SetOptStat(""); gStyle->SetStatY(0.9); @@ -217,357 +197,39 @@ void SimChecker(TString filename = "transfer.root", printf("#####################################################\n"); - // //^################################################ - Int_t size[2] = {padSize,padSize}; ///x,y, single Canvas size - - TCanvas * cCheck[numReact]; - - TH2F * hez[numReact]; - TH2F * hArrayXY[numReact]; - TH2F * hRecoilXY[numReact]; - TH2F * hElum1XY[numReact]; - TH1F ** hThetaCM[numReact]; - TH1F * hExCal[numReact]; - TH1F * hHit[numReact]; - - for( int i = 0; i < numReact; i++ ){ - - TString gateThis = gate + Form(" && rID == %d", i); + Plotter * plotter[numReact]; - cCheck[i] = new TCanvas(Form("cCheck-%d", i), Form("Simulation Checker (%d)", i), 200 * i , 200 * i, size[0]*Div[0], size[1]*Div[1]); - if(cCheck[i]->GetShowEditor() ) cCheck[i]->ToggleEditor(); - if(cCheck[i]->GetShowToolBar() )cCheck[i]->ToggleToolBar(); - cCheck[i]->Divide(Div[0],Div[1]); - - for( int j = 1; j <= Div[0]*Div[1]; j ++){ - cCheck[i]->cd(j); - - if( padPlotID[i-1] == pThetaCM ) { - cCheck[i]->cd(j)->SetGrid(0,0); - cCheck[i]->cd(j)->SetLogy(); - } - - plotID pID = padPlotID[j-1]; - - ///######################################## - if( pID == pEZ){ - hez[i] = new TH2F(Form("hez%d", i), - Form("e-z [gated] @ %5.0f mm; z [mm]; e [MeV]", detGeo.array[i].firstPos), - 200, zRange[i][0], zRange[i][1], - 200, 0, eMax[i]); - - tree->Draw(Form("e:z>>hez%d", i), gateThis , "colz"); - } - - - if( pID == pRecoilXY ){ - hRecoilXY[i] = new TH2F(Form("hRecoilXY%d",i), - Form("RecoilXY [gated] @ %4.0f mm; X [mm]; Y [mm]", detGeo.aux[i].detPos ), - 200, -detGeo.aux[i].outerRadius, detGeo.aux[i].outerRadius, - 200, -detGeo.aux[i].outerRadius, detGeo.aux[i].outerRadius); - tree->Draw(Form("yRecoil:xRecoil>>hRecoilXY%d", i), gateThis, "colz"); - TArc * detArc1 = new TArc(0,0, detGeo.aux[i].outerRadius); - detArc1->SetLineColor(kBlue-8); - detArc1->SetFillStyle(0); - detArc1->Draw("same"); - TArc * detArc2 = new TArc(0,0, detGeo.aux[i].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"); - // } - } - - /* - if( pID == pRecoilXY1 ){ - 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.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.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.aux[detGeoID].outerRadius); - 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 ){ - int numExTemp = reactEx[i].ExList.size(); - hThetaCM[i] = new TH1F *[numExTemp]; - TLegend * legend = new TLegend(0.8,0.2,0.99,0.8); - double maxCount = 0; - for( int h = 0; h < numExTemp; h++){ - hThetaCM[i][h] = new TH1F(Form("hThetaCM%d-%d", i, h), Form("thetaCM [gated] (ExID=%d); thetaCM [deg]; count", h), 200, 0, thetaCMMax); - hThetaCM[i][h]->SetLineColor(h+1); - hThetaCM[i][h]->SetFillColor(h+1); - hThetaCM[i][h]->SetFillStyle(3000+h); - tree->Draw(Form("thetaCM>>hThetaCM%d-%d", i, h), gateThis + Form("&& ExID==%d", h), ""); - legend->AddEntry(hThetaCM[i][h], Form("Ex=%5.1f MeV", reactEx[i].ExList[h].Ex)); - double max = hThetaCM[i][h]->GetMaximum(); - if( max > maxCount ) maxCount = max; - } - - for( int h = 0; h < numExTemp; h++){ - hThetaCM[i][h]->GetYaxis()->SetRangeUser(1, maxCount * 1.2); - if( h == 0 ) { - hThetaCM[i][h]->Draw(); - }else{ - hThetaCM[i][h]->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 < numEx ; i++){ - txList->At(i)->Draw("same"); - } - } - } - */ - - if( pID == pExCal ){ - - double exMin = 9999; - double exMax = -1; - for( size_t k = 0 ; k < reactEx[i].ExList.size(); k++){ - double kuku = reactEx[i].ExList[k].Ex; - if( kuku > exMax ) exMax = kuku; - if( kuku < exMin ) exMin = kuku; - } - - double exPadding = (exMax - exMin) < 1 ? 1 : (exMax - exMin) * 0.3; - - exMin = exMin - exPadding ; - exMax = exMax + exPadding ; - - hExCal[i] = new TH1F(Form("hExCal%d",i), - Form("calculated Ex [gated]; Ex [MeV]; count / %.2f keV", (exMax-exMin)/400.*1000), - 400, exMin, exMax); - tree->Draw(Form("ExCal>>hExCal%d", i), gateThis, ""); - // 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.aux[detGeoID].outerRadius); - tree->Draw("rhoRecoil:thetaCM>>hRecoilRThetaCM", gate, "colz"); - } - */ - if( pID == pArrayXY ){ - - hArrayXY[i] = new TH2F(Form("hArrayXY%d", i), - "Array-XY [gated]; X [mm]; Y [mm]", - 400, -detGeo.array[i].detPerpDist*1.5, detGeo.array[i].detPerpDist*1.5, - 400, -detGeo.array[i].detPerpDist*1.5, detGeo.array[i].detPerpDist*1.5); - tree->Draw(Form("yArray:xArray>>hArrayXY%d", i), gateThis, "colz"); - } - - if( pID == pInfo ){ - TLatex text; - text.SetNDC(); - text.SetTextFont(82); - text.SetTextSize(0.06); - text.SetTextColor(2); - - text.DrawLatex(0., 0.9, reactionList[i].c_str()); - 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)); - // } - } - - - if( pID == pElum1XY ){ - hElum1XY[i] = new TH2F(Form("hElum1XY%d", i), - Form("Elum-1 XY [gated] @ %.0f mm ; X [mm]; Y [mm]", detGeo.aux[i].elumPos1), - 400, -elumMax, elumMax, - 400, -elumMax, elumMax); - tree->Draw(Form("yElum1:xElum1>>hElum1XY%d", i), gateThis, "colz"); - - double count = hElum1XY[i]->GetEntries(); - if( count < 2000. ) { - hElum1XY[i]->SetMarkerStyle(7); - if( count < 500. ) hElum1XY[i]->SetMarkerStyle(3); - hElum1XY[i]->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]); - } - - 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); - - 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"); - - hHit[i] = new TH1F(Form("hHit%d", i), "hit; hit-ID; count", 13, -11, 2); - tree->Draw(Form("hit>>hHit%d", i), "", ""); - } - - }///===== end of pad loop - - }///===== end of reaction loop + for( int i = 0; i < numReact; i++){ + plotter[i] = new Plotter(tree, i, reactionList[i], detGeo, reactEx[i], gate, padPlotID); + plotter[i]->SetRanges(zRange[i][0], zRange[i][1], eMax[i], elumMax, thetaCMMax); + plotter[i]->SetCanvas(colCount, rowCount, 500, padPlotID); + plotter[i]->Plot(); + } + return; } 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 + if( str == "pEZ") return plotID::pEZ; /// 0 + if( str == "pRecoilXY") return plotID::pRecoilXY; /// 1 + if( str == "pThetaCM" ) return plotID::pThetaCM; /// 2 + if( str == "pExCal" ) return plotID::pExCal; /// 3 + if( str == "pArrayXY" ) return plotID::pArrayXY; /// 4 + if( str == "pInfo" ) return plotID::pInfo; /// 5 + if( str == "pElum1XY" ) return plotID::pElum1XY; /// 6 + if( str == "pRecoilXY1" ) return plotID::pRecoilXY1; /// 7 + if( str == "pRecoilXY2" ) return plotID::pRecoilXY2; /// 8 + if( str == "pTDiffZ" ) return plotID::pTDiffZ; /// 9 + if( str == "pRecoilRThetaCM" ) return plotID::pRecoilRThetaCM; /// 10 + if( str == "pRecoilRZ" ) return plotID::pRecoilRZ; /// 11 + if( str == "pEElum1R" ) return plotID::pEElum1R; /// 12 + if( str == "pRecoilRTR" ) return plotID::pRecoilRTR; /// 13 + if( str == "pThetaCM_Z" ) return plotID::pThetaCM_Z; /// 14 + if( str == "pElum1RThetaCM" ) return plotID::pElum1RThetaCM; /// 15 + if( str == "pEmpty" ) return plotID::pEmpty ; /// 16 return plotID::pEmpty; } \ No newline at end of file diff --git a/Cleopatra/SimTransfer.C b/Cleopatra/SimTransfer.C index 6c0cabe..b5f882f 100644 --- a/Cleopatra/SimTransfer.C +++ b/Cleopatra/SimTransfer.C @@ -758,13 +758,14 @@ int main (int argc, char *argv[]) { printf("========== Simulate Transfer reaction in HELIOS ==========\n"); printf("=================================================================\n"); - if(argc == 2 || argc > 5 ) { - printf("Usage: ./Transfer [1] [2] [3] [4]\n"); + if(argc == 2 || argc > 6 ) { + printf("Usage: ./Transfer [1] [2] [3] [4] [5]\n"); printf(" default file name \n"); printf(" [1] reactionConfig.txt (input) reaction Setting \n"); printf(" [2] detectorGeo.txt (input) detector Setting \n"); printf(" [3] DWBA.root (input) thetaCM distribution from DWBA \n"); printf(" [4] transfer.root (output) rootFile name for output \n"); + printf(" [5] 1 (input) 0 = no plot, 1 = plot \n"); printf("-----------------------------------------------------------------\n"); printf(" When DWBA.root provided.\n"); @@ -784,23 +785,26 @@ int main (int argc, char *argv[]) { if( argc >= 3) detGeoFile = argv[2]; if( argc >= 4) ptolemyRoot = argv[3]; if( argc >= 5) saveFileName = argv[4]; + if( argc >= 6) isPlot = atoi(argv[5]); Transfer( basicConfig, detGeoFile, ptolemyRoot, saveFileName); - //run Armory/Check_Simulation - // if( isPlot ){ - // std::ifstream file_in; - // file_in.open("../Cleopatra/Check_Simulation.C", std::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()); + //run Cleopatra/SimChecker.C + if( isPlot ){ + std::ifstream file_in; + file_in.open("../Cleopatra/SimChecker.C", std::ios::in); + if( file_in){ + printf("---- running ../Cleopatra/SimChecker.C on %s \n", saveFileName.Data()); + TString cmd; + cmd.Form("root -l '../Cleopatra/SimChecker.C(\"%s\")'", saveFileName.Data()); - // system(cmd.Data()); - // }else{ - // printf("cannot find ../Cleopatra/Check_Simulation.C \n"); - // } - // } + system(cmd.Data()); + }else{ + printf("cannot find ../Cleopatra/SimChecker.C \n"); + } + } + + return 0; }