247 lines
6.7 KiB
C
247 lines
6.7 KiB
C
#define DecayFinder_cxx
|
|
|
|
#include "DecayFinder.h"
|
|
#include <TH2.h>
|
|
#include <TH1.h>
|
|
#include <TMath.h>
|
|
#include <TCanvas.h>
|
|
#include <TStyle.h>
|
|
#include <TCutG.h>
|
|
#include <TObjArray.h>
|
|
#include <TRandom.h>
|
|
|
|
//===================== Settings
|
|
Bool_t debug = false;
|
|
|
|
TString pidCutFileName = "PIDCuts.root";
|
|
|
|
double distThreshold = 0.2;
|
|
double forwardTime = 250; //ms
|
|
double backwardTime = 100; //ms
|
|
|
|
//===================== histograms
|
|
|
|
TH2F * hPID;
|
|
TH1F ** hDecay ;
|
|
TH2F * hDist;
|
|
TH1F ** hvetoF;
|
|
|
|
TH2F * hIonsPos;
|
|
TH2F * hBetaPos;
|
|
|
|
//===================== Parameters
|
|
|
|
int tick2ns = 8;
|
|
double tick2ms = tick2ns / 1e6;
|
|
|
|
TFile * cutFile;
|
|
TCutG * cut = NULL;
|
|
TObjArray * cutList;
|
|
int numCut = 0;
|
|
|
|
int numMatch = 0;
|
|
|
|
//===================== Begins
|
|
void DecayFinder::Begin(TTree * /*tree*/){
|
|
|
|
TString option = GetOption();
|
|
|
|
hPID = new TH2F("hPID", "PID; TOF; dE", 400, -260, -220, 400, 0, 6500);
|
|
|
|
hDist = new TH2F("hDist", "dist; dx; dy", 100, -distThreshold, distThreshold, 100, -distThreshold, distThreshold);
|
|
|
|
hIonsPos = new TH2F("hIonsPos", "Ions Pos; x; y", 200, 0, 1, 200, 0, 1);
|
|
hBetaPos = new TH2F("hBetaPos", "Beta Pos; x; y", 200, 0, 1, 200, 0, 1);
|
|
|
|
//-------------- GetCut;
|
|
cutFile = new TFile(pidCutFileName, "UPDATE");
|
|
bool listExist = cutFile->GetListOfKeys()->Contains("cutList");
|
|
if( !listExist ) {
|
|
cutList = new TObjArray();
|
|
}else{
|
|
cutList = (TObjArray*) cutFile->FindObjectAny("cutList");
|
|
numCut = cutList->GetLast()+1;
|
|
printf("------------ found %d cuts in %s \n", numCut, pidCutFileName.Data());
|
|
for( int k = 0; k < numCut; k++){
|
|
if( cutList->At(k) != NULL ){
|
|
printf("found a cut at %2d \n", k);
|
|
}else{
|
|
printf(" No cut at %2d \n", k);
|
|
}
|
|
}
|
|
}
|
|
|
|
hDecay = new TH1F *[numCut];
|
|
hvetoF = new TH1F *[numCut];
|
|
|
|
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);
|
|
hvetoF[i] = new TH1F(Form("hvetoF%02d",i), Form("veto-F cut-%02d ; [ms]; count", i), 100, 0, 6000);
|
|
}
|
|
|
|
printf("=====================================\n");
|
|
|
|
}
|
|
|
|
double dist(double x1, double x2, double y1, double y2){
|
|
return sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
|
|
}
|
|
|
|
//===================== Process
|
|
Bool_t DecayFinder::Process(Long64_t entry){
|
|
|
|
if( entry > (Long64_t)1e6 ) return kTRUE;
|
|
|
|
b_TOF->GetEntry(entry);
|
|
b_energy->GetEntry(entry);
|
|
b_crossTime->GetEntry(entry);
|
|
b_crossEnergy->GetEntry(entry);
|
|
b_flag->GetEntry(entry);
|
|
b_vetoFlag->GetEntry(entry);
|
|
b_xIons->GetEntry(entry);
|
|
b_yIons->GetEntry(entry);
|
|
b_dyIonsTime->GetEntry(entry);
|
|
b_veto_r->GetEntry(entry);
|
|
b_veto_f->GetEntry(entry);
|
|
|
|
ULong64_t ImplantEntry = 0;
|
|
ULong64_t ImplantTime = 0;
|
|
double ImplantX = TMath::QuietNaN();
|
|
double ImplantY = TMath::QuietNaN();
|
|
|
|
if( debug && entry % 10 == 0 ) printf("------------- %llu\n", entry);
|
|
|
|
if( flag & 1 ) hPID->Fill(TOF, energy);
|
|
|
|
int cutID = -1;
|
|
for(int i = 0; i < numCut; i++){
|
|
cut = (TCutG*) cutList->At(i);
|
|
if( cut->IsInside(TOF, energy) ) cutID = i;
|
|
}
|
|
|
|
if( cutID == -1 ) return kTRUE;
|
|
if( (flag & 1) == 0 ) return kTRUE; /// no beam
|
|
//if( (flag & 2) == 0 ) return kTRUE; /// no Ions
|
|
if( (flag & 4) == 0 ) return kTRUE; /// no beta
|
|
|
|
//if( veto_f > 0 ) return kTRUE;
|
|
if( veto_r > 0 ) return kTRUE;
|
|
if( dyIonsTime[0] == 0 ) return kTRUE;
|
|
|
|
ImplantEntry = entry;
|
|
ImplantTime = dyIonsTime[0];
|
|
ImplantX = xIons;
|
|
ImplantY = yIons;
|
|
|
|
hIonsPos->Fill( xIons, yIons);
|
|
|
|
if( debug ) printf("===========%8lld, %13llu, %f, %f\n", ImplantEntry, ImplantTime, ImplantX, ImplantY);
|
|
|
|
for( ULong64_t k = ImplantEntry - 10000; k < totNumEntry ; k++){
|
|
b_flag->GetEntry(k);
|
|
b_vetoFlag->GetEntry(k);
|
|
b_dyBetaTime->GetEntry(k);
|
|
b_crossTime->GetEntry(k);
|
|
b_crossEnergy->GetEntry(k);
|
|
b_veto_f->GetEntry(k);
|
|
b_veto_r->GetEntry(k);
|
|
b_xBeta->GetEntry(k);
|
|
b_yBeta->GetEntry(k);
|
|
|
|
if ( k == ImplantEntry ) continue;
|
|
|
|
//if( crossTime > 0 ) continue;
|
|
if( crossEnergy > 0 ) continue;
|
|
//if( (flag & 2) == 2 ) continue; /// has Ions
|
|
if( (flag & 4) == 0 ) continue; /// no Beta
|
|
if( !TMath::IsNaN(veto_f) ) continue;
|
|
if( !TMath::IsNaN(veto_r) ) continue;
|
|
if( dyBetaTime[0] == 0) continue;
|
|
if( dist(ImplantX, xBeta, ImplantY, yBeta) > distThreshold ) continue;
|
|
|
|
if( k < ImplantEntry && ( ImplantTime - dyBetaTime[0] ) * tick2ms > backwardTime) continue;
|
|
if( k > ImplantEntry && ( dyBetaTime[0] - ImplantTime ) * tick2ms > forwardTime ) break;
|
|
|
|
clock.Stop("timer");
|
|
Double_t time = clock.GetRealTime("timer");
|
|
clock.Start("timer");
|
|
printf(" %llu[%.2f%%], searched next %lld events | match case %d | expect : %5.2f min\r",
|
|
ImplantEntry, ImplantEntry*100./totNumEntry, k - ImplantEntry,
|
|
numMatch, totNumEntry*time/(ImplantEntry+1.)/60.);
|
|
|
|
if( debug ) printf(" %8lld, %13llu, %f, %f\n", k, dyBetaTime[0], xBeta, yBeta);
|
|
if( debug ) printf("+++++++++++%8s, %13llu | %.f ms\n", "+", dyBetaTime[0], (dyBetaTime[0] - ImplantTime)*tick2ms);
|
|
|
|
hDist->Fill(ImplantX - xBeta, ImplantY - yBeta);
|
|
|
|
numMatch ++;
|
|
double dT = 0;
|
|
if( dyBetaTime[0] > ImplantTime ){
|
|
dT = (dyBetaTime[0] - ImplantTime)*tick2ms;
|
|
}else{
|
|
dT = (ImplantTime - dyBetaTime[0])*tick2ms;
|
|
dT = -dT;
|
|
}
|
|
|
|
hDecay[cutID]->Fill(dT);
|
|
hvetoF[cutID]->Fill(crossEnergy);
|
|
hBetaPos->Fill(xBeta, yBeta);
|
|
|
|
}
|
|
|
|
return kTRUE;
|
|
}
|
|
|
|
|
|
void DecayFinder::Terminate(){
|
|
|
|
//============== Save the results
|
|
TFile * haha = new TFile("results.root", "recreate");
|
|
haha->cd();
|
|
for( int i = 0; i < numCut; i++){
|
|
hDecay[i]->Write();
|
|
}
|
|
hPID->Write();
|
|
haha->Close();
|
|
|
|
//============== Plot result
|
|
gStyle->SetOptStat("neiou");
|
|
TCanvas * cDecay = new TCanvas("cDecay", "Decay", 1000, 1000);
|
|
cDecay->Divide(2,2);
|
|
|
|
int padID = 1; cDecay->cd(padID);
|
|
hPID->Draw("colz");
|
|
for( int i = 0; i < numCut; i++){
|
|
cut = (TCutG*) cutList->At(i);
|
|
cut->SetLineColor(i+2);
|
|
cut->Draw("same");
|
|
}
|
|
|
|
padID++; cDecay->cd(padID); //cDecay->cd(padID)->SetLogy();
|
|
double yMax = 0;
|
|
for( int i = 0; i < numCut; i++){
|
|
if( hDecay[i]->GetMaximum() > yMax ) yMax = hDecay[i]->GetMaximum();
|
|
}
|
|
for( int i = 0; i < numCut; i++){
|
|
hDecay[i]->SetLineColor(i+2);
|
|
hDecay[i]->SetMaximum(yMax *1.2);
|
|
hDecay[i]->Draw(i == 0 ? "" : "same");
|
|
}
|
|
|
|
//padID++; cDecay->cd(padID);
|
|
//hDist->Draw("colz");
|
|
//
|
|
//padID++; cDecay->cd(padID);
|
|
//for( int i = 0; i < numCut; i++){
|
|
// hvetoF[i]->SetLineColor(i+2);
|
|
// hvetoF[i]->Draw(i == 0 ? "" : "same");
|
|
//}
|
|
|
|
padID++; cDecay->cd(padID);
|
|
hIonsPos->Draw("colz");
|
|
|
|
padID++; cDecay->cd(padID);
|
|
hBetaPos->Draw("colz");
|
|
|
|
}
|