diff --git a/compareTH2D.C b/compareTH2D.C index 5af0ab1..367210f 100644 --- a/compareTH2D.C +++ b/compareTH2D.C @@ -1,10 +1,10 @@ #include "armory/AnalysisLibrary.h" double xMin1, xMax1; -TGraph * fitGraph1 = new TGraph(); +TGraph fitGraph1; Double_t fitFunc1(Double_t *x, Double_t *par) { double pos = par[1]*x[0] + par[0]; - return (xMax1 > pos && pos > xMin1 ? par[2] * fitGraph1->Eval(pos): 0); + return (xMax1 > pos && pos > xMin1 ? par[2] * fitGraph1.Eval(pos): 0); } TF1 * FitHist1(TH1F * hist1, TH1F * hist2, double a0, double a1, double scale){ @@ -15,36 +15,43 @@ TF1 * FitHist1(TH1F * hist1, TH1F * hist2, double a0, double a1, double scale){ xMax1 = hist2->GetXaxis()->GetXmax(); int nBin = hist1->GetNbinsX(); + fitGraph1.Clear(); for( int i = 0; i < nBin; i++){ - fitGraph1->AddPoint(hist1->GetBinCenter(i), hist1->GetBinContent(i)); + fitGraph1.AddPoint(hist1->GetBinCenter(i), hist1->GetBinContent(i)); } double par[3] = {a0, a1, scale}; - TF1 * fit = new TF1("fit", fitFunc1, xMin1, xMax1, 3); - fit->SetParameters(par); - fit->SetLineWidth(2); - fit->SetLineColor(4); - fit->SetNpx(1000); + TF1 * fit1 = new TF1("fit1", fitFunc1, xMin1, xMax1, 3); + fit1->SetParameters(par); + fit1->SetParLimits(0, -20, 15); + fit1->SetParLimits(1, 0.6, 1.3); + //fit1->SetParLimits(2, 0.5*scale, 1.4*scale); + fit1->SetLineWidth(2); + fit1->SetLineColor(4); + fit1->SetNpx(1000); - hist2->Fit("fit", "Rq"); + hist2->Fit("fit1", "Rq"); - return fit; + return fit1; } double xMin2, xMax2; -TGraph * fitGraph2 = new TGraph(); +TGraph fitGraph2; Double_t fitFunc2(Double_t *x, Double_t *par) { - double pos = par[1]*x[0] + par[0]; + + + double pos = par[3]*x[0]*x[0] + par[1]*x[0] + par[0]; //return (xMax2 > pos && pos > xMin2 ? par[2] * fitGraph2->Eval(pos) * exp(pos/1000./par[3]): 0); - return (xMax2 > pos && pos > xMin2 ? par[2] * fitGraph2->Eval(pos) : 0); + + return (xMax2 > pos && pos > xMin2 ? par[2] * fitGraph2.Eval(pos) : 0); - //if( x[0] < 2200 ){ - // return (xMax2 > pos && pos > xMin2 ? par[2] * fitGraph2->Eval(pos): 0); + //if( x[0] < 3900 ){ + // return (xMax2 > pos && pos > xMin2 ? par[2] * fitGraph2.Eval(pos): 0); //}else{ - // return (xMax2 > pos && pos > xMin2 ? par[3] * fitGraph2->Eval(pos): 0); + // return (xMax2 > pos && pos > xMin2 ? par[3] * fitGraph2.Eval(pos): 0); //} @@ -58,25 +65,28 @@ TF1 * FitHist2(TH1F * hist1, TH1F * hist2, double a0, double a1, double scale, d xMax2 = hist2->GetXaxis()->GetXmax(); int nBin = hist1->GetNbinsX(); + fitGraph2.Clear(); for( int i = 0; i < nBin; i++){ - fitGraph2->AddPoint(hist1->GetBinCenter(i), hist1->GetBinContent(i)); + fitGraph2.AddPoint(hist1->GetBinCenter(i), hist1->GetBinContent(i)); } double par[4] = {a0, a1, scale, scale2}; - TF1 * fit = new TF1("fit", fitFunc2, xMin2, xMax2, 4); - fit->SetParameters(par); - fit->SetLineWidth(2); - fit->SetLineColor(4); - fit->SetNpx(1000); + TF1 * fit2 = new TF1("fit2", fitFunc2, xMin2, xMax2, 4); + fit2->SetParameters(par); + fit2->SetParLimits(0, 0, 15); + fit2->SetParLimits(1, 0.6, 1.2); + fit2->SetParLimits(3, -0.00001, 0.00001); + //fit2->FixParameter(3, 0); + fit2->SetLineWidth(2); + fit2->SetLineColor(4); + fit2->SetNpx(1000); - hist2->Fit("fit", "Rq"); + hist2->Fit("fit2", "Rq"); - return fit; + return fit2; } - - void compareTH2D(int run1, int run2){ TFile * f1 = new TFile(Form("PID_%d.root", run1)); @@ -89,16 +99,28 @@ void compareTH2D(int run1, int run2){ TH1F * hdE2 = (TH1F*) f2->Get("hdE"); TH2F * hPID2 = (TH2F*) f2->Get("hPID0"); - TCanvas * cc = new TCanvas("cc", Form("cc %d, %d", run1, run2), 3000, 1000); + TCanvas * cc = (TCanvas *) gROOT->FindObjectAny("cc"); + if( cc == NULL )cc = new TCanvas("cc", Form("cc %d, %d", run1, run2), 3000, 1000); + cc->Clear(); cc->Divide(3,1); //================= TOF cc->cd(1); - hTOF1->Draw(""); - hTOF2->SetLineColor(2); - hTOF2->Draw("same"); - TF1 * fitTOF = FitHist1(hTOF1, hTOF2, 6, 1.0, 0.01); + double yMax1 = hTOF1->GetMaximum(); + double yMax2 = hTOF2->GetMaximum(); + hTOF2->SetLineColor(2); + + double scaleGuess = yMax2/yMax1; + + if( yMax1 > yMax2 ) { + hTOF1->Draw(""); + hTOF2->Draw("same"); + }else{ + hTOF2->Draw(""); + hTOF1->Draw("same"); + } + TF1 * fitTOF = FitHist1(hTOF1, hTOF2, -6, 1.0, scaleGuess); fitTOF->Draw("same"); //TPad * chaha = (TPad*) (cc->cd(1))->Clone(); @@ -120,7 +142,7 @@ void compareTH2D(int run1, int run2){ hdE2->SetLineColor(2); hdE2->Draw("same"); - TF1 * fitdE = FitHist2(hdE1, hdE2, -1, 1, 1, -2); + TF1 * fitdE = FitHist2(hdE1, hdE2, -1, 1, 1, 0.0000001); fitdE->Draw("same"); @@ -134,7 +156,7 @@ void compareTH2D(int run1, int run2){ printf("B2 : %7.4f +- %7.4f \n", paraA_dE[3], paraE_dE[3]); printf("################################\n"); - printf("%3d %9.6f %9.6f %9.6f %9.6f\n", run2, paraA_TOF[0], paraA_TOF[1], paraA_dE[0], paraA_dE[1]); + printf("%3d %9.6f %9.6f %9.6f %9.6f %14.10f\n", run2, paraA_TOF[0], paraA_TOF[1], paraA_dE[0], paraA_dE[1], paraA_dE[3]); printf("\n\n\n"); cc->cd(3); diff --git a/correction_PID.dat b/correction_PID.dat index 8906f79..0844d46 100644 --- a/correction_PID.dat +++ b/correction_PID.dat @@ -1,20 +1,56 @@ -246 5.202374 0.972581 1.823645 0.982412 -247 5.474717 0.973615 2.228469 0.989046 -248 -2.016251 0.974793 3.005200 0.993577 -249 -2.670544 0.972111 2.883073 0.996861 -250 0.000000 1.000000 0.000000 1.000000 -251 -0.657032 0.997220 -0.209729 1.001230 -252 -0.440989 0.998095 -0.235229 1.002683 -253 -0.736023 0.996942 -0.291993 1.003505 -254 8.379539 1.005064 -14.304124 1.011162 -256 12.176462 1.001000 10.480354 0.988767 -257 6.429354 1.010119 10.926955 0.995218 -258 8.853158 1.019825 11.033525 0.998783 -259 15.584682 1.014310 10.904478 1.000129 -261 14.648562 1.010594 9.730146 1.004049 -262 6.638140 1.010922 9.764069 1.005226 -263 7.778790 0.987777 9.657852 1.005952 -264 6.860067 1.011819 9.671465 1.005891 -269 4.648145 1.006698 9.352800 0.921318 -270 12.961422 1.008092 9.634760 0.892950 -271 5.361827 1.010096 9.085634 0.884655 +246 5.202374 0.972581 1.823645 0.982412 0.0 +247 5.474717 0.973615 2.228469 0.989046 0.0 +248 -2.016251 0.974793 3.005200 0.993577 0.0 +249 -2.670544 0.972111 2.883073 0.996861 0.0 +250 0.000000 1.000000 0.000000 1.000000 0.0 +251 -0.657032 0.997220 -0.209729 1.001230 0.0 +252 -0.440989 0.998095 -0.235229 1.002683 0.0 +253 -0.736023 0.996942 -0.291993 1.003505 0.0 +254 8.379539 1.005064 -14.304124 1.011162 0.0 +256 12.176462 1.001000 10.480354 0.988767 0.0 +257 6.429354 1.010119 10.926955 0.995218 0.0 +258 8.853158 1.019825 11.033525 0.998783 0.0 +259 15.584682 1.014310 10.904478 1.000129 0.0 +261 14.410873 1.009665 8.460009 1.008042 -0.0000015982 +262 6.617968 1.010836 8.726306 1.007930 -0.0000009269 +263 13.650833 1.006689 9.110117 1.008211 -0.0000009318 +264 6.860014 1.011819 9.188193 1.008153 -0.0000009184 +269 4.695810 1.006896 4.634574 0.936230 -0.0000056271 +270 12.957748 1.008077 3.095747 0.912196 -0.0000066846 +271 5.362512 1.010098 2.531388 0.905879 -0.0000074197 +283 4.897151 1.008010 10.666408 0.863863 -0.0000090642 +284 5.022649 1.008598 0.689671 0.890106 -0.0000090894 +286 8.558948 1.006674 0.471578 0.889442 -0.0000094556 +287 1.018873 1.008725 -0.007570 0.889010 -0.0000095513 +288 1.068998 1.008322 0.202725 0.889160 -0.0000095366 +292 1.728610 1.010990 0.613623 0.887610 -0.0000090418 +294 1.257935 1.009023 0.682529 0.886398 -0.0000086273 +295 9.308503 1.009003 0.987144 0.887736 -0.0000092990 +296 7.543449 1.001872 0.740777 0.887669 -0.0000093976 +298 1.429204 1.005651 0.365736 0.887218 -0.0000096348 +299 8.355584 1.001215 0.906154 0.886579 -0.0000095789 +306 1.732491 1.007028 0.540370 0.886010 -0.0000096222 +309 1.455301 1.005791 0.511240 0.886757 -0.0000094339 +310 8.305157 1.001002 1.057534 0.886542 -0.0000092710 +311 1.626763 1.006480 0.503249 0.886968 -0.0000094761 +312 8.370479 1.001189 1.224966 0.885002 -0.0000093943 +313 0.868675 1.003133 0.437214 0.886037 -0.0000095128 +314 -7.391993 0.967753 3.602238 0.883946 -0.0000088520 +315 1.453498 1.005403 0.288068 0.887976 -0.0000094549 +316 0.738072 1.002615 0.307353 0.887159 -0.0000091529 +317 1.823295 1.006856 0.481761 0.888012 -0.0000094753 +318 8.209070 1.000661 1.207103 0.887063 -0.0000093263 +319 8.011329 0.999823 1.600204 0.886653 -0.0000091960 +320 0.326489 1.001083 0.552779 0.888471 -0.0000094284 +322 0.241687 1.000738 0.340305 0.889576 -0.0000096254 +323 0.167235 1.000447 0.501191 0.888241 -0.0000095490 +323 -9.981457 0.980308 0.501191 0.888241 -0.0000095490 +324 0.054908 0.999928 0.364381 0.889151 -0.0000095068 +325 8.050909 0.999893 1.242025 0.888227 -0.0000091104 +326 0.250943 1.000691 0.338226 0.889374 -0.0000093779 + + + + + + diff --git a/peachCake.C b/peachCake.C index 3f322be..d3f3f11 100644 --- a/peachCake.C +++ b/peachCake.C @@ -42,8 +42,10 @@ TH1F * hDecay_veto; TH1F * hFlag; TH1F * hvetoFlag; -TH2F * heVrun; -TH2F * htVrun; +TH1F * hZ3; +TH1F * hZ6; +TH1F * hZ10; +TH1F * hZ14; //============ parameters //TString cutFileName="PIDCuts.root"; @@ -99,9 +101,9 @@ void peachCake::Begin(TTree * /*tree*/){ hCloverID = createTH2F("hCloverID", "Clover; ID; energy [keV]", 52, 0, 52, 400, 0, 2000); - hPID = createTH2F("hPID", "PID; ns; msx100", tofRange[0], tofRange[1], tofRange[2], dERange[0], dERange[1], dERange[2]); hPID0 = createTH2F("hPID0", "PID raw; ns; msx100", tofRange[0], tofRange[1], tofRange[2], dERange[0], dERange[1], dERange[2]); - hPID2 = createTH2F("hPID2", "PID2; A/Q; Z", tofRange[0], 2.4, 3.6, dERange[0], 0, 16.5); + hPID = createTH2F("hPID", "PID corrected; ns; msx100", tofRange[0], tofRange[1], tofRange[2], dERange[0], dERange[1], dERange[2]); + hPID2 = createTH2F("hPID2", "PID; A/Q; Z", tofRange[0], 2.0, 4.2, dERange[0], 0, 16.5); hTOF = createTH1F("hTOF", "TOF", tofRange[0], tofRange[1], tofRange[2]); hdE = createTH1F("hdE", "dE", dERange[0], dERange[1], dERange[2]); @@ -123,9 +125,11 @@ void peachCake::Begin(TTree * /*tree*/){ hFlag = createTH1F("hFlag", "Flag ( 1 = beam, 2 = Ions, 4 = Beta)", 8, 0, 8); hvetoFlag = createTH1F("hvetoFlag", "vetoFlag ( 1 = front, 4 = rear)", 8, 0, 8); - - heVrun = createTH2F("heVrun", "energy vs run; run; energy", 200, 237, 257, 400, 0, 2000); - htVrun = createTH2F("hrVrun", "TOF vs run; run; energy", 200, 237, 257, 400, -265, -220); + + hZ3 = createTH1F("hZ3", "Z=3; A ; count ", 200, 6, 12); + hZ6 = createTH1F("hZ6", "Z=6; A ; count ", 200, 14, 21); + hZ10 = createTH1F("hZ10", "Z=10; A ; count ", 200, 25, 33); + hZ14 = createTH1F("hZ14", "Z=14; A ; count ", 200, 38, 44); eCorr = LoadCorrectionParameters("correction_e.dat"); @@ -301,7 +305,7 @@ Bool_t peachCake::Process(Long64_t entry){ if( int(pidCorr[k][0]) != runID/100 ) continue; TOF = pidCorr[k][1] + pidCorr[k][2] * TOF; - energy = pidCorr[k][3] + pidCorr[k][4] * energy; + energy = pidCorr[k][3] + pidCorr[k][4] * energy + pidCorr[k][5] * energy * energy; } } @@ -315,37 +319,45 @@ Bool_t peachCake::Process(Long64_t entry){ hdE->Fill(energy); } - double newTOF = TOF - TOFCorrection(energy) - 240.; + double newTOF = TOF - TOFCorrection(energy) - 234; hPID->Fill(newTOF, energy); - heVrun->Fill(runID/100., energy); - htVrun->Fill(runID/100., TOF); - flag += 1; ///======== Z vs A/Q double c = 299.792458; /// mm/ns - double FL = 37656 ; /// mm - double Brho = 4.90370; /// T.mm + /// exp value + ///double FL = 37656 ; /// mm + ///double Brho = 4.90370; /// T.mm + double FL = 56423 ; /// mm + double Brho = 9.04; /// T.mm double ma = 931.5; double beta0 = Brho * c / TMath::Sqrt( 9*ma*ma + Brho*Brho*c*c); // for AoQ = 3 - double tofOffset = 509.79; /// ns + double tofOffset = 504.69; /// ns double beta = FL/c/(newTOF+tofOffset); double gamma = 1./TMath::Sqrt(1.-beta*beta); - double Z = sqrt((energy + 11.9473) / 23.097) * TMath::Power(beta / beta0,1.3) ; + double Z = sqrt((energy + 11.9473) / 23.097) * TMath::Power(beta / beta0, 1.3) ; double AoQ = c * Brho / gamma/ beta / ma; //printf("tof : %f, beta: %f (%f), Z : %f, A/Q : %f \n", TOF + tofOffset, beta, beta0, Z, AoQ); //Z = -0.15938 + 1.01736 *Z + 0.000203316 * Z*Z; //Z = -0.254632 + 1.06285 * Z - 0.00539634 * Z * Z + 0.000169443 * Z * Z * Z; ///2022-06-09 + + Z = 0.0153343 + 1.00339 * Z; + hPID2->Fill(AoQ, Z); hZ->Fill(Z); + + if( 3.5 > Z && Z > 2.5 ) hZ3->Fill( AoQ * 3); + if( 6.5 > Z && Z > 5.5 ) hZ6->Fill( AoQ * 6); + if( 10.5 > Z && Z > 9.5 ) hZ10->Fill( AoQ * 10); + if( 14.5 > Z && Z > 13.5 ) hZ14->Fill( AoQ * 14); } @@ -471,8 +483,8 @@ void peachCake::Terminate(){ if( plotHists ){ gStyle->SetOptStat(111111); - int div[2] = {3, 2} ; ///x, y - TCanvas * c1 = new TCanvas("c1", "c1", 700 * div[0], 700 * div[1]); + int div[2] = {4, 2} ; ///x, y + TCanvas * c1 = new TCanvas("c1", canvasTitle, 700 * div[0], 700 * div[1]); c1->Divide(div[0], div[1]); @@ -494,7 +506,7 @@ void peachCake::Terminate(){ hPID->Draw("colz"); padID ++; - c1->cd(padID); //c1->cd(padID)->SetLogz(); + c1->cd(padID); c1->cd(padID)->SetLogz(); hPID2->Draw("colz"); padID ++; @@ -503,11 +515,19 @@ void peachCake::Terminate(){ padID ++; c1->cd(padID); //c1->cd(padID)->SetLogz(); - heVrun->Draw("colz"); + hZ3->Draw("colz"); padID ++; c1->cd(padID); //c1->cd(padID)->SetLogz(); - htVrun->Draw("colz"); + hZ6->Draw("colz"); + + padID ++; + c1->cd(padID); //c1->cd(padID)->SetLogz(); + hZ10->Draw("colz"); + + padID ++; + c1->cd(padID); //c1->cd(padID)->SetLogz(); + hZ14->Draw("colz"); //padID ++; //c1->cd(padID); c1->cd(padID)->SetGrid(); c1->cd(padID)->SetLogz(); diff --git a/peachCake.h b/peachCake.h index 03fc0ca..02a5f6c 100644 --- a/peachCake.h +++ b/peachCake.h @@ -52,6 +52,7 @@ public : pidCorrFileName = ""; fHistRootName = ""; plotHists = true; + canvasTitle = ""; } virtual ~peachCake() { } virtual Int_t Version() const { return 2; } @@ -122,6 +123,8 @@ public : bool plotHists; + TString canvasTitle; + }; #endif @@ -149,6 +152,9 @@ void peachCake::Init(TTree *tree){ fChain->SetBranchAddress("multiBeam", &multiBeam, &b_multiplicity_Beam); fChain->SetBranchAddress("runID", &runID, &b_runID); + //============ Canvas Title + + //============ new tree printf(" is save tree? :%d \n", saveNewTree); diff --git a/plotPID.C b/plotPID.C index 89a196f..180258c 100644 --- a/plotPID.C +++ b/plotPID.C @@ -10,7 +10,7 @@ void plotPID() { selector->SetPIDCorrectionFile(""); selector->SetPlotHist(false); - for( int run = 269 ; run < 272 ; run ++){ + for( int run = 272 ; run < 330 ; run ++){ printf("\033[31m######################################### run : %d \033[0m\n", run); diff --git a/script.C b/script.C index 47ccf92..d4178ed 100644 --- a/script.C +++ b/script.C @@ -17,32 +17,23 @@ void script() { //chain->Add("root_data/run-0241-*.root"); // new beam - //chain->Add("root_data/run-0246-*.root"); - //chain->Add("root_data/run-0247-*.root"); - //chain->Add("root_data/run-0248-*.root"); - //chain->Add("root_data/run-0249-*.root"); - //chain->Add("root_data/run-0250-*.root"); - //chain->Add("root_data/run-0251-*.root"); - //chain->Add("root_data/run-0252-*.root"); - //chain->Add("root_data/run-0253-*.root"); - //chain->Add("root_data/run-0254-*.root"); - //chain->Add("root_data/run-0257-*.root"); - //chain->Add("root_data/run-0258-*.root"); - //chain->Add("root_data/run-0259-*.root"); - //chain->Add("root_data/run-0261-*.root"); - //chain->Add("root_data/run-0262-*.root"); - //chain->Add("root_data/run-0263-*.root"); - //chain->Add("root_data/run-0264-*.root"); - //chain->Add("root_data/run-0269-*.root"); - //chain->Add("root_data/run-0270-*.root"); - chain->Add("root_data/run-0271-*.root"); - - int nFile = chain->GetListOfFiles()->GetEntries(); - - printf("================================ num of Files : %d \n", nFile); - - if( nFile == 0 ) return; + chain->Add("root_data/run-024[6-8]*.root"); + chain->Add("root_data/run-0250*.root"); + chain->Add("root_data/run-025[1-4].root"); + chain->Add("root_data/run-025[6-9].root"); + chain->Add("root_data/run-026[1-4]*.root"); + chain->Add("root_data/run-0269*.root"); + chain->Add("root_data/run-027[0-1]*.root"); + chain->Add("root_data/run-028[3-4]*.root"); + chain->Add("root_data/run-028[6-8]*.root"); + chain->Add("root_data/run-0292*.root"); + chain->Add("root_data/run-029[4-6]*.root"); + chain->Add("root_data/run-029[8-9]*.root"); + chain->Add("root_data/run-03*.root"); + + //chain->Add("root_data/run-0250*.root"); + bool isSaveNewTree = false; TString histRootFileName = "";// Form("PID_%03d.root", runNum); @@ -50,6 +41,9 @@ void script() { TString pidCorrFileName = "correction_PID.dat"; chain->GetListOfFiles()->Print(); + int nFile = chain->GetListOfFiles()->GetEntries(); + if( nFile == 0 ) return; + printf("================================ num of Files : %d \n", nFile); printf("================================\n");