rtang/FRIB_e21062
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

snapshot from workstation

Browse Source
This commit is contained in:
Ryan Tang 2023-05-18 16:33:51 -04:00
parent 44087a0df2
commit a1d83534ba
11 changed files with 597 additions and 281 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitignore vendored
View File

@ -3,6 +3,7 @@
*.so *.so
*.root *.root
*.evt *.evt
*.pdf
root_data root_data
EventBuilder EventBuilder

139
DecayFinder.C
View File

@ -13,11 +13,11 @@
//===================== Settings //===================== Settings
Bool_t debug = false; Bool_t debug = false;
TString pidCutFileName = "PIDCuts.root"; TString pidCutFileName = ""; //"PIDCuts.root";
double distThreshold = 0.2; double distThreshold = 0.1;
double forwardTime = 250; //ms double forwardTime = 100; //ms
double backwardTime = 100; //ms double backwardTime = 50; //ms
//===================== histograms //===================== histograms
@ -31,7 +31,7 @@ TH2F * hBetaPos;
//===================== Parameters //===================== Parameters
int tick2ns = 8; int tick2ns = 4;
double tick2ms = tick2ns / 1e6; double tick2ms = tick2ns / 1e6;
TFile * cutFile; TFile * cutFile;
@ -41,12 +41,44 @@ int numCut = 0;
int numMatch = 0; int numMatch = 0;
vector<vector<int>> isotopes = {
{19, 7}, {20, 7}, {21, 7},
{22, 8}, {23, 8}, {24, 8},
{24, 9}, {25, 9}, {26, 9}, {27, 9}, {29, 9},
{27, 10}, {28, 10}, {29, 10}, {30, 10}, {31, 10}, {32, 10},
{29, 11}, {30, 11}, {31, 11}, {32, 11}, {33, 11}, {34, 11}, {35, 11},
{32, 12}, {33, 12}, {34, 12}, {35, 12}, {36, 12}, {37, 12}, {38, 12},
{35, 13}, {36, 13}, {37, 13}, {38, 13}, {39, 13}, {40, 13}, {41, 13},
{39, 14}, {40, 14}, {41, 14}, {42, 14},
{43, 15}, {44, 15}
};
string isoSym(int Z){
switch (Z) {
case 3 : return "Li"; break;
case 4 : return "Be"; break;
case 5 : return "B"; break;
case 6 : return "C"; break;
case 7 : return "N"; break;
case 8 : return "O"; break;
case 9 : return "F"; break;
case 10 : return "Ne"; break;
case 11 : return "Na"; break;
case 12 : return "Mg"; break;
case 13 : return "Al"; break;
case 14 : return "Si"; break;
case 15 : return "P"; break;
case 16 : return "S"; break;
}
return "XX";
}
//===================== Begins //===================== Begins
void DecayFinder::Begin(TTree * /*tree*/){ void DecayFinder::Begin(TTree * /*tree*/){
TString option = GetOption(); TString option = GetOption();
hPID = new TH2F("hPID", "PID; TOF; dE", 400, -260, -220, 400, 0, 6500); hPID = new TH2F("hPID", "PID (beam & Implant & back veto); A/Q; Z", 400, 2.2, 4.0, 400, 3, 16);
hDist = new TH2F("hDist", "dist; dx; dy", 100, -distThreshold, distThreshold, 100, -distThreshold, distThreshold); hDist = new TH2F("hDist", "dist; dx; dy", 100, -distThreshold, distThreshold, 100, -distThreshold, distThreshold);
@ -54,10 +86,35 @@ void DecayFinder::Begin(TTree * /*tree*/){
hBetaPos = new TH2F("hBetaPos", "Beta Pos; x; y", 200, 0, 1, 200, 0, 1); hBetaPos = new TH2F("hBetaPos", "Beta Pos; x; y", 200, 0, 1, 200, 0, 1);
//-------------- GetCut; //-------------- GetCut;
cutFile = new TFile(pidCutFileName, "UPDATE"); cutFile = new TFile(pidCutFileName, "READ");
bool listExist = cutFile->GetListOfKeys()->Contains("cutList"); bool listExist = cutFile->GetListOfKeys()->Contains("cutList");
if( !listExist ) { if( !listExist ) {
cutList = new TObjArray(); cutList = new TObjArray();
numCut = (int) isotopes.size();
printf("============= no Cut found, creating new default Cuts\n");
for( int i = 0 ; i < numCut ; i++){
cut = new TCutG();
printf("%2d | A: %2d, Z: %d \n", i, isotopes[i][0], isotopes[i][1]);
TString name; name.Form("cut%02d%s", isotopes[i][0], isoSym(isotopes[i][1]).c_str());
cut->SetName(name);
cut->SetVarX("AoQ");
cut->SetVarY("Z");
cut->SetTitle(Form("cut%2d", i));
cut->SetLineColor(i+1);
cut->SetPoint(0, (isotopes[i][0] - 0.45)/isotopes[i][1], isotopes[i][1] - 0.45);
cut->SetPoint(1, (isotopes[i][0] - 0.45)/isotopes[i][1], isotopes[i][1] + 0.45);
cut->SetPoint(2, (isotopes[i][0] + 0.45)/isotopes[i][1], isotopes[i][1] + 0.45);
cut->SetPoint(3, (isotopes[i][0] + 0.45)/isotopes[i][1], isotopes[i][1] - 0.45);
cut->SetPoint(4, (isotopes[i][0] - 0.45)/isotopes[i][1], isotopes[i][1] - 0.45);
cutList->AddAtAndExpand(cut, i);
}
}else{ }else{
cutList = (TObjArray*) cutFile->FindObjectAny("cutList"); cutList = (TObjArray*) cutFile->FindObjectAny("cutList");
numCut = cutList->GetLast()+1; numCut = cutList->GetLast()+1;
@ -75,8 +132,8 @@ void DecayFinder::Begin(TTree * /*tree*/){
hvetoF = new TH1F *[numCut]; hvetoF = new TH1F *[numCut];
for( int i = 0; i < numCut; i++){ for( int i = 0; i < numCut; i++){
hDecay[i] = new TH1F(Form("hDecay%02d",i), Form("Decay cut-%02d ; [ms]; count", i), 50, -backwardTime, forwardTime); hDecay[i] = new TH1F(Form("hDecay%02d%s",isotopes[i][0], isoSym(isotopes[i][1]).c_str()), Form("Decay cut-%02d ; [ms]; count", i), 200, -backwardTime, forwardTime);
hvetoF[i] = new TH1F(Form("hvetoF%02d",i), Form("veto-F cut-%02d ; [ms]; count", i), 100, 0, 6000); hvetoF[i] = new TH1F(Form("hvetoF%02d%s",isotopes[i][0], isoSym(isotopes[i][1]).c_str()), Form("veto-F cut-%02d ; [ms]; count", i), 100, 0, 6000);
} }
printf("=====================================\n"); printf("=====================================\n");
@ -90,10 +147,10 @@ double dist(double x1, double x2, double y1, double y2){
//===================== Process //===================== Process
Bool_t DecayFinder::Process(Long64_t entry){ Bool_t DecayFinder::Process(Long64_t entry){
if( entry > (Long64_t)1e6 ) return kTRUE; //if( entry > (Long64_t)1000 ) return kTRUE;
b_TOF->GetEntry(entry); b_Z->GetEntry(entry);
b_energy->GetEntry(entry); b_AoQ->GetEntry(entry);
b_crossTime->GetEntry(entry); b_crossTime->GetEntry(entry);
b_crossEnergy->GetEntry(entry); b_crossEnergy->GetEntry(entry);
b_flag->GetEntry(entry); b_flag->GetEntry(entry);
@ -111,28 +168,34 @@ Bool_t DecayFinder::Process(Long64_t entry){
if( debug && entry % 10 == 0 ) printf("------------- %llu\n", entry); if( debug && entry % 10 == 0 ) printf("------------- %llu\n", entry);
if( flag & 1 ) hPID->Fill(TOF, energy); if( (flag == 3) && (vetoFlag & 2) == 0 ) hPID->Fill(AoQ, Z);
int cutID = -1; int cutID = -1;
for(int i = 0; i < numCut; i++){ for(int i = 0; i < numCut; i++){
cut = (TCutG*) cutList->At(i); cut = (TCutG*) cutList->At(i);
if( cut->IsInside(TOF, energy) ) cutID = i; if( cut->IsInside(AoQ, Z) ) cutID = i;
} }
if( cutID == -1 ) return kTRUE; if( cutID == -1 ) return kTRUE;
if( (flag & 1) == 0 ) return kTRUE; /// no beam
//if( (flag & 2) == 0 ) return kTRUE; /// no Ions if( flag >= 4 ) return kTRUE;
if( (flag & 4) == 0 ) return kTRUE; /// no beta ///if( flag != 3 ) return kTRUE; ///
//if( (flag & 1) == 0 ) return kTRUE; /// no beam
//if( (flag & 2) == 0 ) return kTRUE; /// no Ions
//if( (flag & 4) == 4 ) return kTRUE; /// has beta
//if( veto_f > 0 ) return kTRUE; //if( veto_f > 0 ) return kTRUE;
if( veto_r > 0 ) return kTRUE; //if( veto_r > 0 ) return kTRUE;
if( dyIonsTime[0] == 0 ) return kTRUE; //if( dyIonsTime[0] == 0 ) return kTRUE;
ImplantEntry = entry; ImplantEntry = entry;
ImplantTime = dyIonsTime[0]; ImplantTime = dyIonsTime[0];
ImplantX = xIons; ImplantX = xIons;
ImplantY = yIons; ImplantY = yIons;
if( TMath::IsNaN(xIons) ) return kTRUE;
if( TMath::IsNaN(yIons) ) return kTRUE;
hIonsPos->Fill( xIons, yIons); hIonsPos->Fill( xIons, yIons);
if( debug ) printf("===========%8lld, %13llu, %f, %f\n", ImplantEntry, ImplantTime, ImplantX, ImplantY); if( debug ) printf("===========%8lld, %13llu, %f, %f\n", ImplantEntry, ImplantTime, ImplantX, ImplantY);
@ -151,12 +214,13 @@ Bool_t DecayFinder::Process(Long64_t entry){
if ( k == ImplantEntry ) continue; if ( k == ImplantEntry ) continue;
//if( crossTime > 0 ) continue; //if( crossTime > 0 ) continue;
if( crossEnergy > 0 ) continue; //if( crossEnergy > 0 ) continue;
//if( (flag & 2) == 2 ) continue; /// has Ions if( (flag & 2) == 2 ) continue; /// has Ions
//if( flag != 4 ) continue; /// no Beta
if( (flag & 4) == 0 ) continue; /// no Beta if( (flag & 4) == 0 ) continue; /// no Beta
if( !TMath::IsNaN(veto_f) ) continue; if( !TMath::IsNaN(veto_f) ) continue;
if( !TMath::IsNaN(veto_r) ) continue; if( !TMath::IsNaN(veto_r) ) continue;
if( dyBetaTime[0] == 0) continue; //if( dyBetaTime[0] == 0) continue;
if( dist(ImplantX, xBeta, ImplantY, yBeta) > distThreshold ) continue; if( dist(ImplantX, xBeta, ImplantY, yBeta) > distThreshold ) continue;
if( k < ImplantEntry && ( ImplantTime - dyBetaTime[0] ) * tick2ms > backwardTime) continue; if( k < ImplantEntry && ( ImplantTime - dyBetaTime[0] ) * tick2ms > backwardTime) continue;
@ -165,7 +229,7 @@ Bool_t DecayFinder::Process(Long64_t entry){
clock.Stop("timer"); clock.Stop("timer");
Double_t time = clock.GetRealTime("timer"); Double_t time = clock.GetRealTime("timer");
clock.Start("timer"); clock.Start("timer");
printf(" %llu[%.2f%%], searched next %lld events | match case %d | expect : %5.2f min\r", printf(" %llu[%4.2f%%], searched next %lld events | match case %d | expect : %5.2f min\r",
ImplantEntry, ImplantEntry*100./totNumEntry, k - ImplantEntry, ImplantEntry, ImplantEntry*100./totNumEntry, k - ImplantEntry,
numMatch, totNumEntry*time/(ImplantEntry+1.)/60.); numMatch, totNumEntry*time/(ImplantEntry+1.)/60.);
@ -184,7 +248,7 @@ Bool_t DecayFinder::Process(Long64_t entry){
} }
hDecay[cutID]->Fill(dT); hDecay[cutID]->Fill(dT);
hvetoF[cutID]->Fill(crossEnergy); //hvetoF[cutID]->Fill(crossEnergy);
hBetaPos->Fill(xBeta, yBeta); hBetaPos->Fill(xBeta, yBeta);
} }
@ -212,12 +276,12 @@ void DecayFinder::Terminate(){
int padID = 1; cDecay->cd(padID); int padID = 1; cDecay->cd(padID);
hPID->Draw("colz"); hPID->Draw("colz");
for( int i = 0; i < numCut; i++){ for( int i = 0; i < numCut; i++){
cut = (TCutG*) cutList->At(i); cut = (TCutG*) cutList->At(i);
cut->SetLineColor(i+2); cut->SetLineColor(i+2);
cut->Draw("same"); cut->Draw("same");
} }
padID++; cDecay->cd(padID); //cDecay->cd(padID)->SetLogy(); padID++; cDecay->cd(padID); cDecay->cd(padID)->SetLogy();
double yMax = 0; double yMax = 0;
for( int i = 0; i < numCut; i++){ for( int i = 0; i < numCut; i++){
if( hDecay[i]->GetMaximum() > yMax ) yMax = hDecay[i]->GetMaximum(); if( hDecay[i]->GetMaximum() > yMax ) yMax = hDecay[i]->GetMaximum();
@ -225,7 +289,7 @@ void DecayFinder::Terminate(){
for( int i = 0; i < numCut; i++){ for( int i = 0; i < numCut; i++){
hDecay[i]->SetLineColor(i+2); hDecay[i]->SetLineColor(i+2);
hDecay[i]->SetMaximum(yMax *1.2); hDecay[i]->SetMaximum(yMax *1.2);
hDecay[i]->Draw(i == 0 ? "" : "same"); hDecay[i]->Draw(i == 0 ? "" : "same");
} }
//padID++; cDecay->cd(padID); //padID++; cDecay->cd(padID);
@ -242,5 +306,20 @@ void DecayFinder::Terminate(){
padID++; cDecay->cd(padID); padID++; cDecay->cd(padID);
hBetaPos->Draw("colz"); hBetaPos->Draw("colz");
//================ Save historgram
TString fHistRootName = "zzz_all_decay.root";
if( true ){
TFile * fHist = new TFile(fHistRootName, "recreate");
fHist->cd();
for( int i = 0 ; i < numCut; i++){
hDecay[i]->Write();
}
fHist->Close();
printf("---- Save PID histogram as %s \n", fHistRootName.Data());
}
} }

25
DecayFinder.h
View File

@ -24,9 +24,11 @@ public :
// Declaration of leaf types // Declaration of leaf types
ULong64_t eventID; ULong64_t eventID;
//Int_t runID; Int_t runID;
Double_t TOF; //Double_t TOF;
Double_t energy; //Double_t energy;
Double_t Z;
Double_t AoQ;
UInt_t crossEnergy; UInt_t crossEnergy;
ULong64_t crossTime; ULong64_t crossTime;
Short_t flag; Short_t flag;
@ -42,9 +44,11 @@ public :
// List of branches // List of branches
TBranch *b_eventID; //! TBranch *b_eventID; //!
//TBranch *b_runID; //! TBranch *b_runID; //!
TBranch *b_TOF; //! //TBranch *b_TOF; //!
TBranch *b_energy; //! //TBranch *b_energy; //!
TBranch *b_Z; //!
TBranch *b_AoQ; //!
TBranch *b_crossTime; //! TBranch *b_crossTime; //!
TBranch *b_crossEnergy; //! TBranch *b_crossEnergy; //!
TBranch *b_flag; //! TBranch *b_flag; //!
@ -94,9 +98,12 @@ void DecayFinder::Init(TTree *tree){
fChain->SetMakeClass(1); fChain->SetMakeClass(1);
fChain->SetBranchAddress("eventID", &eventID, &b_eventID); fChain->SetBranchAddress("eventID", &eventID, &b_eventID);
//fChain->SetBranchAddress("runID", &runID, &b_runID); fChain->SetBranchAddress("runID", &runID, &b_runID);
fChain->SetBranchAddress("TOF", &TOF, &b_TOF); //fChain->SetBranchAddress("TOF", &TOF, &b_TOF);
fChain->SetBranchAddress("energy", &energy, &b_energy); //fChain->SetBranchAddress("energy", &energy, &b_energy);
fChain->SetBranchAddress("Z", &Z, &b_Z);
fChain->SetBranchAddress("AoQ", &AoQ, &b_AoQ);
fChain->SetBranchAddress("crossTime", &crossTime, &b_crossTime);
fChain->SetBranchAddress("crossTime", &crossTime, &b_crossTime); fChain->SetBranchAddress("crossTime", &crossTime, &b_crossTime);
fChain->SetBranchAddress("crossEnergy", &crossEnergy, &b_crossEnergy); fChain->SetBranchAddress("crossEnergy", &crossEnergy, &b_crossEnergy);
fChain->SetBranchAddress("flag", &flag, &b_flag); fChain->SetBranchAddress("flag", &flag, &b_flag);

71
Hist_Style/PID_paper.C Normal file
View File

@ -0,0 +1,71 @@
void PID_paper(double gamma, int start= 1){
TFile * f1 = new TFile("pid_total.root");
TH2F * hPID = (TH2F*) f1->Get("hPID2");
gStyle->SetOptStat(""); /// no state
gStyle->SetOptTitle(0); /// no title
TCanvas * ccc = new TCanvas("ccc", "ccc", 1000, 1000);
///ccc->Divide(2,1);
///========= set Z-color
///---- custom color
///Int_t colors[] = {0, 1, 2, 3, 4, 5, 6};
///gStyle->SetPalette((sizeof(colors)/sizeof(Int_t)), colors);
///----- build-in color
int color = kBeach;
//int color = kGreenBrownTerrain;
//int color = kColorPrintableOnGrey;
//int color = kSunset;
///gStyle->SetPalette(kBlueRedYellow);
///gStyle->SetPalette(kRainBow);
gStyle->SetPalette(color);TColor::InvertPalette();
///======== set contours
TGraph * ggg = new TGraph(); /// this is for displaying the gamma correction
const int nCon = 200; /// number of contours
double contours[nCon];
double zMax = hPID->GetMaximum();
for( int i = 0; i < nCon; i++){
double xpos = i*1.0/nCon;
double ypos = TMath::Power(xpos, gamma) *0.95;
contours[i] = start + ypos*(zMax - start);
ggg->AddPoint(xpos, ypos);
}
ccc->cd();ccc->cd()->SetGrid(0,0);
hPID->GetYaxis()->SetRangeUser(1.5, 16);
hPID->SetContour(nCon, contours);
hPID->Draw("col");
TEllipse * si42 = new TEllipse(3.01, 14, 0.05, 0.55);
si42->SetFillStyle(0);
si42->SetLineColor(2);
si42->SetLineWidth(3);
si42->Draw("same");
TLatex text;
text.SetTextColor(2);
text.SetTextSize(0.05);
text.DrawLatex(3.05, 14.2, "^{42}Si");
///------ in case want to see the gamma correction
///ccc->cd(2);
///ggg->Draw("AP*");
///====== save as pdf
ccc->SaveAs(Form("haha_%3.1f_%d_%d.pdf", gamma, start, color));
}

109
Hist_Style/ipid.C Normal file
View File

@ -0,0 +1,109 @@
{
TFile * f1 = new TFile("PIDFile.root");
TH2F * hPID = (TH2F*) f1->Get("ipid");
gStyle->SetOptStat(""); /// no state
gStyle->SetOptTitle(0); /// no title
TCanvas * ccc = new TCanvas("ccc", "ccc", 1000, 1000);
const int nBinX = hPID->GetNbinsX();
const int nBinY = hPID->GetNbinsY();
const double xBinWidth = hPID->GetXaxis()->GetBinWidth(1);
const double yBinWidth = hPID->GetYaxis()->GetBinWidth(1);
double skewPar[6] = {10.9179,
0.00416034,
-233.306,
3607.16,
1538.37,
0.000430609};
const double c = 299.792458; /// mm/ns
const double FL = 56423 ; /// mm
const double Brho = 9.04; /// T.mm
const double ma = 931.5;
const double beta0 = Brho * c / TMath::Sqrt( 9*ma*ma + Brho*Brho*c*c); // for AoQ = 3
const double tofOffset = 504.69; /// ns
TH2F * hPID2 = new TH2F("hPID2", "PID; ToF; dE", 500, -280, -210, 500, 0, 5500);
TH2F * hZAoQ = new TH2F("hZAoQ", "PID; A/Q; Z", 400, 2.0, 4.0, 400, 1, 17);
TH1F * hZ = new TH1F("hZ", "Z; Z; count", 500, 1, 17);
for( int i = 1; i < nBinX ; i++){
for( int j = 1; j < nBinY ; j++){
double xpos = hPID->GetXaxis()->GetBinCenter(i) + (gRandom->Rndm()-0.5) * xBinWidth;
double ypos = hPID->GetYaxis()->GetBinCenter(j) + (gRandom->Rndm()-0.5) * yBinWidth;
double zpos = hPID->GetBinContent(i,j);
///===== Skew correction
double ToF = xpos - (skewPar[0] * exp( - skewPar[1] * ypos ) + skewPar[2] - skewPar[5] * (skewPar[3] - ypos) * exp( - pow( (skewPar[3]-ypos)/skewPar[4],2) ) ) - 234 + 12;
hPID2->Fill( ToF, ypos, zpos);
double beta = FL/c/(ToF + tofOffset);
double gamma = 1./TMath::Sqrt(1.-beta*beta);
double Z = sqrt((ypos + 11.9473) / 23.097) * TMath::Power(beta / beta0, 1.3) ;
double AoQ = c * Brho / gamma/ beta / ma;
Z = -0.732942 + 1.47488 * Z - 0.100088 * Z*Z + 0.009867 * Z*Z*Z - 0.000476 * Z*Z*Z*Z + 8.90519e-6 * Z*Z*Z*Z*Z;
hZAoQ->Fill(AoQ, Z, zpos);
hZ->Fill(Z, zpos);
}
}
double gamma = 4.0;
int start = 0;
int color = kColorPrintableOnGrey;
gStyle->SetPalette(color);TColor::InvertPalette();
double zMax = hZAoQ->GetMaximum();
const int nCon = 200; /// number of contours
double contours[nCon];
for( int i = 0; i < nCon; i++){
double xpos = i*1.0/nCon;
double ypos = TMath::Power(xpos, gamma) *0.95;
contours[i] = start + ypos*(zMax - start);
}
///ccc->Divide(2,1);
ccc->cd(1); ///ccc->cd(1)->SetLogz();
hZAoQ->SetContour(nCon, contours);
hZAoQ->Draw("col");
TEllipse * si42 = new TEllipse(3.01, 14, 0.05, 0.55);
si42->SetFillStyle(0);
si42->SetLineColor(2);
si42->SetLineWidth(3);
si42->Draw("same");
TLatex text;
text.SetTextColor(2);
text.SetTextSize(0.05);
text.DrawLatex(3.05, 14.2, "^{42}Si");
//ccc->cd(2);
//hZ->Draw("hist");
ccc->SaveAs(Form("iPID_%3.1f_%d_%d.pdf", gamma, start, color));
}

8
PIDCutCreator.C
View File

@ -45,7 +45,7 @@ void PIDCutCreator(TH2F* hist,
for( int k = 0; k < numCut; k++){ for( int k = 0; k < numCut; k++){
if( cutList->At(k) != NULL ){ if( cutList->At(k) != NULL ){
printf("found a cut at %2d \n", k); printf("found a cut at %2d \n", k);
cut = (TCutG*) cutList->At(k); cut = (TCutG*) cutList->At(k);
cut->Draw("same"); cut->Draw("same");
}else{ }else{
printf(" No cut at %2d \n", k); printf(" No cut at %2d \n", k);
@ -68,10 +68,10 @@ void PIDCutCreator(TH2F* hist,
break; break;
} }
TString name; name.Form("cut%dd", countCut); TString name; name.Form("cut%2d", countCut);
cut->SetName(name); cut->SetName(name);
cut->SetVarX("TOF"); cut->SetVarX("AoQ");
cut->SetVarY("dE"); cut->SetVarY("Z");
cut->SetTitle(Form("cut%2d", countCut)); cut->SetTitle(Form("cut%2d", countCut));
cut->SetLineColor(countCut+1); cut->SetLineColor(countCut+1);
printf(" cut-%d \n", countCut); printf(" cut-%d \n", countCut);

122
armory/AutoFit.C
View File

@ -1897,40 +1897,91 @@ void fitNGaussPolSub(TH1F * hist, int degPol, TString fitFile = "AutoFit_para.t
} }
//######################################## //########################################
//###### fit A * Exp(-t/T) + B //###### fit B + A * Exp(-t/Th1) + A * Th2/(Th1-Th2)*( Exp(-t/Th1) - Exp(-t/Th2) )
//######################################## //########################################
void fitDecay(TH1F * hist){ Double_t decayFunc(Double_t *x, Double_t *par){
///"[0]*TMath::Exp(-x/[1])+ [2]* [3]/([1]-[3]) * ( Exp[-x/[1]] - Exp[-x/[3]]) + [4]"
double ln2 = par[5];
if ( x[0] < 0 ) {
return par[4];
}else{
double d1 = par[0]*TMath::Exp(-x[0]/par[1]*ln2);
double d2 = par[2]*(TMath::Exp(-x[0]/par[1]*ln2) - TMath::Exp(-x[0]/par[3]*ln2)) * par[3] / (par[1]-par[3]);
return d1 + d2 +par[4];
}
}
void fitDecay(TH1F * hist, double decayConst1 = 0, double decayConst2 = 0, bool halfLive = false){
printf("===============================================================\n");
printf("================= fit decay + conseqent decay ================\n");
printf(" A1 Exp(-t/T1) + A2 * T2/(T1-T2) * ( Exp(-t/T1) - Exp(-t/T2) ) + B ");
printf(" decayConst1 = T1 \n");
printf(" decayConst2 = T2 : when 0, auto guess, when < 0, set A2 = 0\n");
printf(" halfLive = false by default, can switch on to change T1, T2 to be half-live\n");
printf(" \n");
printf("================================================================\n");
gStyle->SetOptStat(""); gStyle->SetOptStat("");
TCanvas * cFitDecay = NULL; TCanvas * cFitDecay = NULL;
if( gROOT->FindObjectAny("cFitDecay") == NULL ){ if( gROOT->FindObjectAny("cFitDecay") == NULL ){
cFitDecay = new TCanvas("cFitDecay", Form("fit A * Exp(-t/T) + B | fitDecay"), 800, 1000); cFitDecay = new TCanvas("cFitDecay", Form("fit consqeunce decay | fitDecay"), 800, 1000);
}else{ }else{
delete gROOT->FindObjectAny("cFitDecay") ; delete gROOT->FindObjectAny("cFitDecay") ;
cFitDecay = new TCanvas("cFitDecay", Form("fit A * Exp(-t/T) + B | fitDecay"), 800, 1000); cFitDecay = new TCanvas("cFitDecay", Form("fit consqeunce decay | fitDecay"), 800, 1000);
} }
cFitDecay->Divide(1, 2); cFitDecay->Divide(1, 2);
double xMin = hist->GetXaxis()->GetXmin();
double xMax = hist->GetXaxis()->GetXmax();
double yMax = hist->GetBinContent(hist->GetMaximumBin()); int xBinYMax = hist->GetMaximumBin();
double yMax = hist->GetBinContent(xBinYMax);
double bg = hist->GetBinContent(1);
hist->SetMaximum(yMax*1.2); hist->SetMaximum(yMax*1.2);
cFitDecay->cd(1); hist->Draw(); cFitDecay->cd(1); hist->Draw();
double estHalf = hist->GetBinCenter(hist->FindLastBinAbove(yMax/2.)) *1.3; double estHalf1 = (decayConst1 == 0 ? hist->GetBinCenter(hist->FindLastBinAbove((yMax-bg)/2.)) *1.3 : decayConst1);
double estHalf2 = (decayConst2 == 0 ? estHalf1 * 3 : decayConst2);
if( decayConst2 >= 0 ){
printf("estimator : A1 = A2 = %.2f, T1 = %.2f, T2 = %.2f, bg = %.2f \n", yMax, estHalf1, estHalf2, bg);
}else{
printf("estimator : A1 = %.2f, T1 = %.2f, A2 = 0, T2 = 1, bg = %.2f \n", yMax, estHalf1, bg);
}
//printf("estimator : A = %.2f, T = %.2f \n", yMax, estHalf); TF1 * fit = new TF1("fit", decayFunc, xMin, xMax, 6 );
TF1 * fit = new TF1("fit", "[0]*TMath::Exp(-x/[1]/TMath::Log(2))+[2]");
fit->SetLineWidth(3); fit->SetLineWidth(3);
fit->SetLineColor(1); fit->SetLineColor(1);
fit->SetNpx(1000); fit->SetNpx(1000);
fit->SetParameter(0, yMax ); fit->SetParameter(0, yMax );
fit->SetParameter(1, estHalf ); fit->SetParameter(1, estHalf1 );
fit->SetParameter(2, 1); fit->SetParameter(2, yMax );
fit->SetParameter(3, estHalf2 );
fit->SetParLimits(0, 0, 1e9); fit->SetParameter(4, bg);
if( halfLive == false ) {
fit->FixParameter(5, 1 );
}else{
fit->FixParameter(5, TMath::Log(2) );
}
fit->SetParLimits(0, 1, 1e9);
fit->SetParLimits(1, 0, 1e9); fit->SetParLimits(1, 0, 1e9);
fit->SetParLimits(2, 0, 1e9); fit->SetParLimits(4, 0, 1e9);
///======= put A2 = 0 for decayConst < 0
if( decayConst2 < 0 ) {
fit->FixParameter(2, 0);
fit->FixParameter(3, 1);
}else{
fit->SetParLimits(2, 0, 1e9);
fit->SetParLimits(3, 0, 1e9);
}
gStyle->SetOptFit(000000); gStyle->SetOptFit(000000);
hist->Fit("fit", "q"); hist->Fit("fit", "q");
@ -1943,32 +1994,53 @@ void fitDecay(TH1F * hist){
text.SetTextFont(82); text.SetTextFont(82);
text.SetTextSize(0.04); text.SetTextSize(0.04);
text.SetTextColor(1); text.SetTextColor(1);
if( decayConst2 >= 0 ){
text.DrawLatex(0.6, 0.8, "A e^{#left(-#frac{t}{T Ln(2)}#right)}+B"); text.DrawLatex(0.3, 0.8, "A1 e^{#left(-#frac{t}{T1}#right)} + A2 #frac{T2}{T1-T2} (e^{#left(-#frac{t}{T1}#right)}-e^{#left(-#frac{t}{T2}#right)})+B");
text.DrawLatex(0.6, 0.75, Form("A : %.f(%.f)", paraA[0], paraE[0])); }else{
text.DrawLatex(0.6, 0.7, Form("T : %.f(%.f)", paraA[1], paraE[1])); text.DrawLatex(0.3, 0.8, "A1 e^{#left(-#frac{t}{T1}#right)} +B");
text.DrawLatex(0.6, 0.65, Form("B : %.f(%.f)", paraA[2], paraE[2])); }
text.DrawLatex(0.3, 0.75, Form("A1 : %.2f(%.2f)", paraA[0], paraE[0]));
text.DrawLatex(0.6, 0.75, Form("T1 : %.2f(%.2f)", paraA[1], paraE[1]));
if( decayConst2 >= 0 ){
text.DrawLatex(0.3, 0.70, Form("A2 : %.2f(%.2f)", paraA[2], paraE[2]));
text.DrawLatex(0.6, 0.70, Form("T2 : %.2f(%.2f)", paraA[3], paraE[3]));
}
text.DrawLatex(0.6, 0.65, Form(" B : %.2f(%.2f)", paraA[4], paraE[4]));
double chi2 = fit->GetChisquare(); double chi2 = fit->GetChisquare();
int ndf = fit->GetNDF(); int ndf = fit->GetNDF();
text.DrawLatex(0.6, 0.6, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf)); double sigma = TMath::Sqrt(paraA[4]);
text.DrawLatex(0.6, 0.60, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf/sigma));
TF1 * fite = new TF1("fite", "[0]*TMath::Exp(-x/[1])" ); TF1 * fite = new TF1("fite", "[0]*TMath::Exp(-x/[1])", 0, xMax );
fite->SetLineColor(4);
fite->SetLineWidth(2);
fite->SetParameter(0, paraA[0]); fite->SetParameter(0, paraA[0]);
fite->SetParameter(1, paraA[1]); fite->SetParameter(1, paraA[1]);
fite->Draw("same"); fite->Draw("same");
if( decayConst2 >= 0 ){
TF1 * fitd = new TF1("fitd", "[0]*(TMath::Exp(-x/[1])-TMath::Exp(-x/[2]))*[2]/([1]-[2])", 0, xMax );
fitd->SetLineColor(6);
fitd->SetLineWidth(2);
fitd->SetParameter(0, paraA[2]);
fitd->SetParameter(1, paraA[1]);
fitd->SetParameter(2, paraA[3]);
fitd->Draw("same");
}
TF1 * fitbg = new TF1("fitbg", "pol0" ); TF1 * fitbg = new TF1("fitbg", "pol0", xMin, xMax );
fitbg->SetParameter(0, paraA[2]); fitbg->SetLineColor(1);
fitbg->SetLineWidth(1);
fitbg->SetParameter(0, paraA[4]);
fitbg->Draw("same"); fitbg->Draw("same");
cFitDecay->cd(2); cFitDecay->cd(2);
PlotResidual(hist, fit); PlotResidual(hist, fit);
} }
//################################################# //#################################################
//#### fit N-Gauss with pol-n BG using mouse click //#### fit N-Gauss with pol-n BG using mouse click
//################################################# //#################################################

261
peachCake.C
View File

@ -40,10 +40,10 @@ TH1F * hDecay_veto;
TH1F * hFlag; TH1F * hFlag;
TH1F * hvetoFlag; TH1F * hvetoFlag;
TH1F * hZ3; ///TH1F * hZ3;
TH1F * hZ6; ///TH1F * hZ6;
TH1F * hZ10; ///TH1F * hZ10;
TH1F * hZ14; ///TH1F * hZ14;
//============ parameters //============ parameters
//TString cutFileName="PIDCuts.root"; //TString cutFileName="PIDCuts.root";
@ -54,9 +54,8 @@ TH1F * hZ14;
vector<vector<double>> eCorr; vector<vector<double>> eCorr;
double TOFCorrection(double x){ double TOFCorrection(double x){
///========= for run 250
// for run 250
double par[6] = {10.9179, double par[6] = {10.9179,
0.00416034, 0.00416034,
-233.306, -233.306,
@ -64,7 +63,6 @@ double TOFCorrection(double x){
1538.37, 1538.37,
0.000430609}; 0.000430609};
return par[0] * exp( - par[1] * x ) + par[2] - par[5] * (par[3] - x) * exp( - pow( (par[3]-x)/par[4],2) ); return par[0] * exp( - par[1] * x ) + par[2] - par[5] * (par[3] - x) * exp( - pow( (par[3]-x)/par[4],2) );
} }
@ -83,15 +81,14 @@ void peachCake::Begin(TTree * /*tree*/){
hZ = createTH1F("hZ", "Z; Z; count", 500, 0, 16); hZ = createTH1F("hZ", "Z; Z; count", 500, 0, 16);
hImplantIons = createTH2F("hImplantIons", "Implat low gain; x ; y", 400, 0, 1, 400, 0, 1); hImplantIons = createTH2F("hImplantIons", "Implant low gain; x ; y", 400, 0, 1, 400, 0, 1);
hImplantBeta = createTH2F("hImplantBeta", "Implat high gain; x ; y", 400, 0, 1, 400, 0, 1); hImplantBeta = createTH2F("hImplantBeta", "Implant high gain; x ; y", 400, 0, 1, 400, 0, 1);
hImplantX = createTH2F("hImplantX", "Implat X; low ; high", 400, 0, 1, 400, 0, 1); hImplantX = createTH2F("hImplantX", "Implant X; low ; high", 400, 0, 1, 400, 0, 1);
hImplantY = createTH2F("hImplantY", "Implat Y; low ; high", 400, 0, 1, 400, 0, 1); hImplantY = createTH2F("hImplantY", "Implant Y; low ; high", 400, 0, 1, 400, 0, 1);
hImplantDyIons = createTH2F("hImplantDyIonslow", "Implant low; sum_a; dy", 400, 0, 16000, 400, 0, 16000);
hImplantDyIons = createTH2F("hImplantDyIonsow", "Implat low; sum_a; dy", 400, 0, 80000, 400, 0, 80000); hImplantDyBeta = createTH2F("hImplantDyBetahigh", "Implant high; sum_a; dy", 400, 0, 80000, 400, 0, 80000);
hImplantDyBeta = createTH2F("hImplantDyBetaigh", "Implat high; sum_a; dy", 400, 0, 8000, 400, 0, 8000);
hDecay = createTH1F("hDecay", "Decay (beta.time - ion.time) ; [ticks]; counts", 100, 0, 2000); hDecay = createTH1F("hDecay", "Decay (beta.time - ion.time) ; [ticks]; counts", 100, 0, 2000);
hDecay_veto = createTH1F("hDecay_veto", "Decay (beta.time - ion.time) [veto]; [ticks]; counts", 100, 0, 2000); hDecay_veto = createTH1F("hDecay_veto", "Decay (beta.time - ion.time) [veto]; [ticks]; counts", 100, 0, 2000);
@ -99,10 +96,10 @@ void peachCake::Begin(TTree * /*tree*/){
hFlag = createTH1F("hFlag", "Flag ( 1 = beam, 2 = Ions, 4 = Beta)", 8, 0, 8); hFlag = createTH1F("hFlag", "Flag ( 1 = beam, 2 = Ions, 4 = Beta)", 8, 0, 8);
hvetoFlag = createTH1F("hvetoFlag", "vetoFlag ( 1 = front, 4 = rear)", 8, 0, 8); hvetoFlag = createTH1F("hvetoFlag", "vetoFlag ( 1 = front, 4 = rear)", 8, 0, 8);
hZ3 = createTH1F("hZ3", "Z=3; A ; count ", 200, 6, 12); ///hZ3 = createTH1F("hZ3", "Z=3; A ; count ", 200, 6, 12);
hZ6 = createTH1F("hZ6", "Z=6; A ; count ", 200, 14, 21); ///hZ6 = createTH1F("hZ6", "Z=6; A ; count ", 200, 14, 21);
hZ10 = createTH1F("hZ10", "Z=10; A ; count ", 200, 25, 33); ///hZ10 = createTH1F("hZ10", "Z=10; A ; count ", 200, 25, 33);
hZ14 = createTH1F("hZ14", "Z=14; A ; count ", 200, 38, 44); ///hZ14 = createTH1F("hZ14", "Z=14; A ; count ", 200, 38, 44);
eCorr = LoadCorrectionParameters("correction_e.dat"); eCorr = LoadCorrectionParameters("correction_e.dat");
@ -127,14 +124,18 @@ Bool_t peachCake::Process(Long64_t entry){
b_runID->GetEntry(entry); b_runID->GetEntry(entry);
//========= initialization //========= initialization
energy = TMath::QuietNaN(); dE = TMath::QuietNaN();
TOF = TMath::QuietNaN();
Z = TMath::QuietNaN();
AoQ = TMath::QuietNaN();
Long64_t startTimeL = 0; Long64_t startTimeL = 0;
Long64_t startTimeR = 0; Long64_t startTimeR = 0;
Long64_t stopTime = 0;
double cfdL = TMath::QuietNaN(); double cfdL = TMath::QuietNaN();
double cfdR = TMath::QuietNaN(); double cfdR = TMath::QuietNaN();
double cfd2 = TMath::QuietNaN(); double cdB2 = TMath::QuietNaN();
Long64_t stopTime = 0;
Long64_t startTime40 = 0; Long64_t startTime40 = 0;
double cfd40 = TMath::QuietNaN(); double cfd40 = TMath::QuietNaN();
@ -181,40 +182,40 @@ Bool_t peachCake::Process(Long64_t entry){
} }
} }
//======= Scanning the event //======= Fill detectors
for( int i = 0; i < multi; i++){ for( int i = 0; i < multi; i++){
//----- clover //----- clover
if( 0 <= detID[i] && detID[i] < 100 ){ ///if( 0 <= detID[i] && detID[i] < 100 ){
if( e[i] > 0 ) { /// if( e[i] > 0 ) {
int id = detID[i]; /// int id = detID[i];
double eCal = 0; /// double eCal = 0;
for( int k = 0; k < int(eCorr[id].size()); k++){ /// for( int k = 0; k < int(eCorr[id].size()); k++){
eCal += eCorr[id][k] * TMath::Power(e[i], k); /// eCal += eCorr[id][k] * TMath::Power(e[i], k);
} /// }
hCloverID->Fill(detID[i], eCal); /// hCloverID->Fill(detID[i], eCal);
} /// }
} ///}
//----- dy Beta //----- dy Ions (low-gain)
if( detID[i] == 200 ) { if( detID[i] == 200 ) {
if ( multiDyBeta < 4 ){ if ( multiDyBeta < 4 ){
dyBetaTime[multiDyBeta] = e_t[i]; dyIonsTime[multiDyBeta] = e_t[i];
dyBetaEnergy[multiDyBeta] = e[i]; dyIonsEnergy[multiDyBeta] = e[i];
multiDyBeta ++;
}
}
if( 201 <= detID[i] && detID[i] < 250){
a_H[detID[i] - 200] = e[i];
}
//----- dy Ions
if( detID[i] == 250 ) {
if ( multiDyIons < 4 ){
dyIonsTime[multiDyIons] = e_t[i];
dyIonsEnergy[multiDyIons] = e[i];
multiDyIons ++; multiDyIons ++;
} }
} }
if( 201 <= detID[i] && detID[i] < 250){
a_L[detID[i] - 200] = e[i];
}
//----- dy Beta (high-gain)
if( detID[i] == 250 ) {
if ( multiDyIons < 4 ){
dyBetaTime[multiDyIons] = e_t[i];
dyBetaEnergy[multiDyIons] = e[i];
multiDyBeta ++;
}
}
if( 251 <= detID[i] && detID[i] < 260){ if( 251 <= detID[i] && detID[i] < 260){
a_L[detID[i] - 250] = e[i]; a_H[detID[i] - 250] = e[i];
} }
//----- veto_front //----- veto_front
if( detID[i] == 290 ) { if( detID[i] == 290 ) {
@ -249,13 +250,13 @@ Bool_t peachCake::Process(Long64_t entry){
//----- cross B2 logic //----- cross B2 logic
if( detID[i] == 307 ) { if( detID[i] == 307 ) {
stopTime = e_t[i]; stopTime = e_t[i];
cfd2 = cfd[i]/16384.; cdB2 = cfd[i]/16384.;
if( (vetoFlag & 4) == 0) vetoFlag += 4; if( (vetoFlag & 4) == 0) vetoFlag += 4;
} }
//----- MSX40 //----- MSX40
///if( detID[i] == 308 ) energy = e[i]; ///if( detID[i] == 308 ) dE = e[i];
//----- MSX100 //----- MSX100
if( detID[i] == 309 ) energy = e[i]; if( detID[i] == 309 ) dE = e[i];
//----- MSX40 logic //----- MSX40 logic
if( detID[i] == 310 ) { if( detID[i] == 310 ) {
startTime40 = e_t[i]; startTime40 = e_t[i];
@ -267,38 +268,33 @@ Bool_t peachCake::Process(Long64_t entry){
Long64_t startTime = startTimeL; Long64_t startTime = startTimeL;
double startCFD = cfdL; double startCFD = cfdL;
TOF = TMath::QuietNaN(); if( startTime > 0 && stopTime > 0 && dE > 0 ){
if( startTime > 0 && stopTime > 0 && energy > 0 ){
if( stopTime > startTime ){ if( stopTime > startTime ){
TOF = double(stopTime - startTime) - startCFD + cfd2; TOF = double(stopTime - startTime) - startCFD + cdB2;
}else{ }else{
TOF = -double(startTime - stopTime) - startCFD + cfd2 ; TOF = -double(startTime - stopTime) - startCFD + cdB2 ;
} }
TOF = TOF * 8; // to ns TOF = TOF * 8; /// to ns
///============== PID correction ///============== PID correction
if( pidCorrFileName != "" ){ if( pidCorrFileName != "" ){
for( int k = 0; k < (int) pidCorr.size() ; k++){ for( int k = 0; k < (int) pidCorr.size() ; k++){
if( int(pidCorr[k][0]) != runID/100 ) continue; if( int(pidCorr[k][0]) != runID/100 ) continue;
TOF = pidCorr[k][1] + pidCorr[k][2] * TOF; TOF = pidCorr[k][1] + pidCorr[k][2] * TOF;
energy = pidCorr[k][3] + pidCorr[k][4] * energy + pidCorr[k][5] * energy * energy; dE = pidCorr[k][3] + pidCorr[k][4] * dE + pidCorr[k][5] * dE * dE;
} }
} }
/// TOF correction, A/Q = 3 to have a fixed TOF. /// TOF correction, A/Q = 3 to have a fixed TOF.
if( dE > 160./44.*(TOF + 260) +180 ) {
if( energy > 160./44.*(TOF + 260) +180 ) { hPID0->Fill(TOF, dE);
hPID0->Fill(TOF, energy);
hTOF->Fill(TOF); hTOF->Fill(TOF);
hdE->Fill(energy); hdE->Fill(dE);
} }
double newTOF = TOF - TOFCorrection(energy) - 234; ///--------- slew correction
double newTOF = TOF - TOFCorrection(dE) - 234;
hPID->Fill(newTOF, energy); hPID->Fill(newTOF, dE);
flag += 1; flag += 1;
@ -313,18 +309,18 @@ Bool_t peachCake::Process(Long64_t entry){
double beta = FL/c/(newTOF+tofOffset); double beta = FL/c/(newTOF+tofOffset);
double gamma = 1./TMath::Sqrt(1.-beta*beta); double gamma = 1./TMath::Sqrt(1.-beta*beta);
double Z = sqrt((energy + 11.9473) / 23.097) * TMath::Power(beta / beta0, 1.3) ; Z = sqrt((dE + 11.9473) / 23.097) * TMath::Power(beta / beta0, 1.3) ;
double AoQ = c * Brho / gamma/ beta / ma; AoQ = c * Brho / gamma/ beta / ma;
Z = 0.0153343 + 1.00339 * Z; Z = 0.0153343 + 1.00339 * Z;
hPID2->Fill(AoQ, Z); hPID2->Fill(AoQ, Z);
hZ->Fill(Z); hZ->Fill(Z);
if( 3.5 > Z && Z > 2.5 ) hZ3->Fill( AoQ * 3); ///if( 3.5 > Z && Z > 2.5 ) hZ3->Fill( AoQ * 3);
if( 6.5 > Z && Z > 5.5 ) hZ6->Fill( AoQ * 6); ///if( 6.5 > Z && Z > 5.5 ) hZ6->Fill( AoQ * 6);
if( 10.5 > Z && Z > 9.5 ) hZ10->Fill( AoQ * 10); ///if( 10.5 > Z && Z > 9.5 ) hZ10->Fill( AoQ * 10);
if( 14.5 > Z && Z > 13.5 ) hZ14->Fill( AoQ * 14); ///if( 14.5 > Z && Z > 13.5 ) hZ14->Fill( AoQ * 14);
} }
@ -341,8 +337,10 @@ Bool_t peachCake::Process(Long64_t entry){
xBeta = (a_H[3]+a_H[2])/sumA_H; xBeta = (a_H[3]+a_H[2])/sumA_H;
yBeta = (a_H[1]+a_H[2])/sumA_H; yBeta = (a_H[1]+a_H[2])/sumA_H;
if( 0 < xIons && xIons < 1 && 0 < yIons && yIons < 1 ) flag += 2;
if( 0 < xBeta && xBeta < 1 && 0 < yBeta && yBeta < 1 ) flag += 4;
if( (flag & 1) == 0 ) { if( (flag & 1) == 0 ) { ///has beam
hImplantIons->Fill(xIons, yIons); hImplantIons->Fill(xIons, yIons);
hImplantBeta->Fill(xBeta, yBeta); hImplantBeta->Fill(xBeta, yBeta);
hImplantDyIons->Fill(sumA_L, dyIonsEnergy[0]); hImplantDyIons->Fill(sumA_L, dyIonsEnergy[0]);
@ -353,17 +351,15 @@ Bool_t peachCake::Process(Long64_t entry){
} }
if( !TMath::IsNaN(veto_f) && !TMath::IsNaN(veto_r) && crossTime == 0 ){ if( !TMath::IsNaN(veto_f) && !TMath::IsNaN(veto_r) && crossTime == 0 ){
//hImplantIons->Fill(xIons, yIons); ///hImplantIons->Fill(xIons, yIons);
//hImplantBeta->Fill(xBeta, yBeta); ///hImplantBeta->Fill(xBeta, yBeta);
//hImplantDyIons->Fill(sumA_L, dyIonsEnergy[0]); ///hImplantDyIons->Fill(sumA_L, dyIonsEnergy[0]);
//hImplantDyBeta->Fill(sumA_H, dyBetaEnergy[0]); ///hImplantDyBeta->Fill(sumA_H, dyBetaEnergy[0]);
// ///
//hImplantX->Fill(xIons, xBeta); ///hImplantX->Fill(xIons, xBeta);
//hImplantY->Fill(yIons, yBeta); ///hImplantY->Fill(yIons, yBeta);
} }
if( 0 < xIons && xIons < 1 && 0 < yIons && yIons < 1 ) flag += 2;
if( 0 < xBeta && xBeta < 1 && 0 < yBeta && yBeta < 1 ) flag += 4;
hFlag->Fill(flag); hFlag->Fill(flag);
hvetoFlag->Fill(vetoFlag); hvetoFlag->Fill(vetoFlag);
@ -375,9 +371,9 @@ Bool_t peachCake::Process(Long64_t entry){
printf("################# %lld, multi:%d\n", entry, multi); printf("################# %lld, multi:%d\n", entry, multi);
printf("----- beam Line \n"); printf("----- beam Line \n");
printf(" MSX100 eneergy : %f \n", energy); printf(" MSX100 eneergy : %f \n", dE);
printf(" startTime : %llu, CFD : %f \n", startTime, startCFD); printf(" startTime : %llu, CFD : %f \n", startTime, startCFD);
printf(" stopTime : %llu, CFD : %f \n", stopTime, cfd2); printf(" stopTime : %llu, CFD : %f \n", stopTime, cdB2);
printf(" TOF : %f ns \n", TOF); printf(" TOF : %f ns \n", TOF);
printf("----- Ions, low gain \n"); printf("----- Ions, low gain \n");
@ -426,6 +422,7 @@ Bool_t peachCake::Process(Long64_t entry){
//============================= Rejection //============================= Rejection
if ( flag == 0 ) return kTRUE; if ( flag == 0 ) return kTRUE;
if ( Z < 0 ) return kTRUE;
//============================= Save Tree //============================= Save Tree
if( saveNewTree ){ if( saveNewTree ){
@ -456,80 +453,50 @@ void peachCake::Terminate(){
printf("|%s|\n", canvasTitle.Data()); printf("|%s|\n", canvasTitle.Data());
TCanvas * c1 = new TCanvas("c1", canvasTitle, 700 * div[0], 700 * div[1]); TCanvas * c1 = new TCanvas("c1", canvasTitle, 700 * div[0], 700 * div[1]);
c1->Divide(div[0], div[1]); c1->Divide(div[0], div[1]);
int padID=0;
///padID ++;
///c1->cd(padID); //c1->cd(padID)->SetLogz();
///hPID0->Draw("colz");
int padID=1; ///padID ++;
///c1->cd(padID); //c1->cd(padID)->SetLogz();
c1->cd(padID); //c1->cd(padID)->SetLogz(); ///hPID->Draw("colz");
hPID0->Draw("colz");
//padID ++;
//c1->cd(padID); //c1->cd(padID)->SetLogz();
//hTOF->Draw();
//
//padID ++;
//c1->cd(padID); //c1->cd(padID)->SetLogz();
//hdE->Draw();
padID ++;
c1->cd(padID); //c1->cd(padID)->SetLogz();
hPID->Draw("colz");
padID ++; padID ++;
c1->cd(padID); c1->cd(padID)->SetLogz(); c1->cd(padID); c1->cd(padID)->SetLogz();
hPID2->Draw("colz"); hPID2->Draw("colz");
padID ++;
c1->cd(padID); //c1->cd(padID)->SetLogz();
hZ->Draw("colz");
padID ++;
c1->cd(padID); //c1->cd(padID)->SetLogz();
hZ3->Draw("colz");
padID ++;
c1->cd(padID); //c1->cd(padID)->SetLogz();
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 ++; //padID ++;
//c1->cd(padID); c1->cd(padID)->SetGrid(); c1->cd(padID)->SetLogz(); //c1->cd(padID); c1->cd(padID)->SetGrid(); c1->cd(padID)->SetLogz();
//hCloverID->Draw("colz"); hCloverID->SetNdivisions(-13); //hCloverID->Draw("colz"); hCloverID->SetNdivisions(-13);
/*
padID ++;
c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantIons->Draw("colz");
padID ++; padID ++;
c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantBeta->Draw("colz"); c1->cd(padID); c1->cd(padID)->SetLogz();
hImplantIons->Draw("colz");
//padID ++;
//c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantBeta_veto->Draw("colz"); padID ++;
c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantBeta->Draw("colz");
padID ++;
c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantDyIons->Draw("colz"); //padID ++;
//c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantBeta_veto->Draw("colz");
padID ++;
c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantDyBeta->Draw("colz"); padID ++;
c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantDyIons->Draw("colz");
padID ++;
c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantX->Draw("colz"); padID ++;
c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantDyBeta->Draw("colz");
padID ++;
c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantY->Draw("colz"); padID ++;
c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantX->Draw("colz");
padID ++;
c1->cd(padID); hFlag->Draw(""); padID ++;
hvetoFlag->SetLineColor(2); hvetoFlag->Draw("same"); c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantY->Draw("colz");
*/ padID ++;
c1->cd(padID); hFlag->Draw("");
hvetoFlag->SetLineColor(2); hvetoFlag->Draw("same");
} }

96
peachCake.h
View File

@ -97,8 +97,8 @@ public :
Int_t eventID; Int_t eventID;
UInt_t crossEnergy; UInt_t crossEnergy;
ULong64_t crossTime; ULong64_t crossTime;
Double_t TOF; Double_t TOF, dE;
Double_t energy; Double_t Z, AoQ;
Short_t flag; /// 1 = has beam, 2 = has Ions, 4 = hasBeta, default 0; Short_t flag; /// 1 = has beam, 2 = has Ions, 4 = hasBeta, default 0;
@ -161,51 +161,61 @@ void peachCake::Init(TTree *tree){
printf(" is save tree? :%d \n", saveNewTree); printf(" is save tree? :%d \n", saveNewTree);
if( saveNewTree ){ if( saveNewTree ){
saveFileName = "test"; saveFileName = "zzz_";
int oldSeg = -100; int oldRun = -1;
int numFile = fChain->GetListOfFiles()->GetLast()+1; bool flag1 = false;
for( int i = 0; i < numFile; i++){ int numFile = fChain->GetListOfFiles()->GetLast()+1;
TString name = fChain->GetListOfFiles()->At(i)->GetTitle(); for( int i = 0; i < numFile; i++){
int pos = name.Last('/'); TString name = fChain->GetListOfFiles()->At(i)->GetTitle();
name = name.Remove(0,pos+1); //this should be run-XXXX-XX.root int pos = name.Last('/');
if( i == 0 ) { name = name.Remove(0,pos+1); ///this should be run-XXXX-XX.root
saveFileName = name; pos = name.Last('-');
pos = saveFileName.Last('-'); name = name.Remove(pos);
saveFileName = saveFileName.Remove(pos); //run-XXXX pos = name.Last('-');
name = name.Remove(0, pos+2);
int runNum = atoi(name);
if( runNum == oldRun + 1 ) {
int len = saveFileName.Sizeof();
if( flag1 == false) {
saveFileName += "-";
}
if( flag1 == true) saveFileName.Remove(len-4);
flag1 = true;
}
if( runNum > oldRun + 1 ) {
flag1 = false;
saveFileName += "_";
}
if( runNum > oldRun ){
saveFileName += Form("%03d", runNum );
oldRun = runNum;
}
} }
pos = name.Last('-'); saveFileName = saveFileName + ".root";
int segNum = atoi(name.Remove(0, pos+1).Remove(2));
if( oldSeg == -100 ) oldSeg = segNum;
saveFileName = saveFileName + Form("_%02d",segNum);
if( segNum == oldSeg + 1 && i != numFile -1 ){
int len = saveFileName.Length();
saveFileName = saveFileName.Remove(len - 3);
oldSeg = segNum;
}
}
saveFileName = saveFileName + ".root";
printf("=========== saveFile : %s \n", saveFileName.Data()); printf("=========== saveFile : %s \n", saveFileName.Data());
saveFile = new TFile(saveFileName, "recreate"); saveFile = new TFile(saveFileName, "recreate");
newTree = new TTree("tree", "tree"); newTree = new TTree("tree", "tree");
newTree->Branch("eventID", &eventID, "eventID/l"); newTree->Branch("eventID", &eventID, "eventID/l");
//newTree->Branch("runID", &runID, "runID/I"); newTree->Branch("runID", &runID, "runID/I");
newTree->Branch("TOF", &TOF, "TOF/D"); ///newTree->Branch("TOF", &TOF, "TOF/D");
newTree->Branch("energy", &energy, "energy/D"); ///newTree->Branch("dE", &dE, "dE/D");
newTree->Branch("crossTime", &crossTime, "crossTime/l"); newTree->Branch("AoQ", &AoQ, "AoQ/D");
newTree->Branch("crossEnergy", &crossEnergy, "crossEnergy/i"); newTree->Branch("Z", &Z, "Z/D");
newTree->Branch("flag", &flag, "flag/S"); newTree->Branch("crossTime", &crossTime, "crossTime/l");
newTree->Branch("vetoFlag", &vetoFlag, "vetoFlag/S"); newTree->Branch("crossEnergy", &crossEnergy, "crossEnergy/i");
newTree->Branch("xIons", &xIons, "xIons/D"); newTree->Branch("flag", &flag, "flag/S");
newTree->Branch("yIons", &yIons, "yIons/D"); newTree->Branch("vetoFlag", &vetoFlag, "vetoFlag/S");
newTree->Branch("xBeta", &xBeta, "xBeta/D"); newTree->Branch("xIons", &xIons, "xIons/D");
newTree->Branch("yBeta", &yBeta, "yBeta/D"); newTree->Branch("yIons", &yIons, "yIons/D");
newTree->Branch("dyIonsTime", dyIonsTime, "dyIonsTime[4]/l"); newTree->Branch("xBeta", &xBeta, "xBeta/D");
newTree->Branch("dyBetaTime", dyBetaTime, "dyBetaTime[4]/l"); newTree->Branch("yBeta", &yBeta, "yBeta/D");
newTree->Branch("veto_f", &veto_f, "veto_f/D"); newTree->Branch("dyIonsTime", dyIonsTime, "dyIonsTime[4]/l");
newTree->Branch("veto_r", &veto_r, "veto_r/D"); newTree->Branch("dyBetaTime", dyBetaTime, "dyBetaTime[4]/l");
newTree->Branch("veto_f", &veto_f, "veto_f/D");
newTree->Branch("veto_r", &veto_r, "veto_r/D");
} }
//==== clock //==== clock

2
plotPID.C
View File

@ -31,7 +31,7 @@ void plotPID() {
if( nFile == 0 ) continue; if( nFile == 0 ) continue;
selector->SetHistRootFileName(Form("PID_%03d.root", run)); selector->SetHistRootFileName(Form("PID_runID/PID_%03d.root", run));
chain->Process(selector, ""); chain->Process(selector, "");

44
script.C
View File

@ -17,31 +17,31 @@ void script() {
//chain->Add("root_data/run-0241-*.root"); //chain->Add("root_data/run-0241-*.root");
// new beam // new beam
//chain->Add("root_data/run-024[6-8]*.root"); chain->Add("root_data/run-024[6-8]*.root");
//chain->Add("root_data/run-0250*.root"); chain->Add("root_data/run-0250*.root");
//chain->Add("root_data/run-025[1-4].root"); chain->Add("root_data/run-025[1-4].root");
//chain->Add("root_data/run-025[6-9].root"); chain->Add("root_data/run-025[6-9].root");
//chain->Add("root_data/run-026[1-4]*.root"); chain->Add("root_data/run-026[1-4]*.root");
//chain->Add("root_data/run-0269*.root"); chain->Add("root_data/run-0269*.root");
//chain->Add("root_data/run-027[0-1]*.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[3-4]*.root");
//chain->Add("root_data/run-028[6-8]*.root"); chain->Add("root_data/run-028[6-8]*.root");
//chain->Add("root_data/run-0292*.root"); chain->Add("root_data/run-0292*.root");
//chain->Add("root_data/run-029[4-6]*.root"); chain->Add("root_data/run-029[4-6]*.root");
//chain->Add("root_data/run-029[8-9]*.root"); chain->Add("root_data/run-029[8-9]*.root");
//chain->Add("root_data/run-031*.root"); chain->Add("root_data/run-031*.root");
//chain->Add("root_data/run-032*.root"); chain->Add("root_data/run-032*.root");
chain->Add("root_data/run-0246-00.root"); //chain->Add("root_data/run-0246-00.root");
chain->Add("root_data/run-0246-01.root"); //chain->Add("root_data/run-0246-01.root");
chain->Add("root_data/run-0247-00.root"); //chain->Add("root_data/run-0247-00.root");
chain->Add("root_data/run-0247-01.root"); //chain->Add("root_data/run-0247-01.root");
chain->Add("root_data/run-0248-00.root"); //chain->Add("root_data/run-0248-00.root");
chain->Add("root_data/run-0250-00.root"); //chain->Add("root_data/run-0250-00.root");
chain->Add("root_data/run-0250-01.root"); //chain->Add("root_data/run-0250-01.root");
bool isSaveNewTree = false; bool isSaveNewTree = true;
TString histRootFileName = "";// Form("PID_%03d.root", runNum); TString histRootFileName = "";// Form("PID_%03d.root", runNum);