replace splitpole.C/h to SplitPolePlotter.C, to TreeReader

This commit is contained in:
Ryan Tang 2024-07-19 15:00:29 -04:00
parent 00ddf3dcf5
commit 3fc682f918
6 changed files with 36 additions and 312 deletions

View File

@ -83,6 +83,7 @@ int main(int argc, char **argv) {
printf(" Include Trace = %s\n", traceOn ? "Yes" : "No");
printf(" Debug level = %d\n", debug);
printf(" Max multiplity = %d hits/event (hard coded)\n", MAX_MULTI);
if( traceOn ) printf(" Max Trace Length = %d (hard coded)\n", MAX_TRACE_LENGTH);
printf("========================================= Grouping files\n");
std::vector<std::vector<FileInfo>> fileGroupList; // fileName and ID = SN * 1000 + index

View File

@ -47,7 +47,8 @@ TH1F * hXavg;
TH2F * hFocal;
TH2F * hXavgVQ;
TH2F * hXavg_Q;
TH2F * hXavg_Theta;
TH2F * haha;
@ -56,7 +57,7 @@ TH1F * hEx;
ULong64_t t1, t2;
#define XMIN -20
#define XMAX 100
#define XMAX 150
//^###########################################
@ -82,7 +83,7 @@ void SplitPolePlotter(TChain *tree){
//*====================================================== histograms
PID = new TH2F("hPID", "PID; Scin_X ; AnodeB", 200, 0, 30000, 100, 0, 70000);
PID = new TH2F("hPID", "PID; Scin_X ; AnodeB", 200, 0, 2000, 100, 0, 4000);
coin = new TH2F("hCoin", "Coincident ", 16, 0, 16, 16, 0, 16);
hMulti = new TH1F("hMulti", "Multiplicity", 16, 0, 16);
@ -92,14 +93,15 @@ void SplitPolePlotter(TChain *tree){
hXavg = new TH1F("hAvg", "Xavg", 600, XMIN, XMAX);
hFocal = new TH2F("hFocal", "Front vs Back ", 200, XMIN, XMAX, 200, XMIN, XMAX);
hXavgVQ = new TH2F("hXavgVQ", "Xavg vs Q ", 200, XMIN, XMAX, 200, 0, 40000);
hXavg_Q = new TH2F("hXavg_Q", "Xavg vs Q ", 200, XMIN, XMAX, 200, 0, 5000);
hXavg_Theta = new TH2F("hXavg_Theta", "Xavg vs Theta ", 200, XMIN, XMAX, 200, 2.5, 3);
haha = new TH2F("haha", "", 400, XMIN, XMAX, 400, -50, 50);
hEx = new TH1F("hEx", "Ex; Ex [MeV]; count/100 keV", 250, -5, 20);
hit.SetMassTablePath("../analyzers/mass20.txt");
hit.CalConstants("12C", "12C", "4He", 80, 5); // 80MeV, 5 deg
hit.CalConstants("36S", "d", "p", 80, 12); // 80MeV, 5 deg
hit.CalZoffset(1.41); // 1.41 T
t1 = 0;
@ -142,21 +144,23 @@ void SplitPolePlotter(TChain *tree){
hit.ClearData();
hMulti->Fill(sn.GetSize());
// if( multi.Get()[0] != 9 ) continue;
for( int i = 0; i < sn.GetSize(); i++){
t2 = e_t[i];
if( t2 < t1 ) printf("entry %lld-%d, timestamp is not in order. %llu, %llu\n", processedEntries, i, t2, t1);
if( i == 0 ) t1 = e_t[i];
if( ch[i] == ChMap::ScinR ) {hit.eSR = e[i]; hit.tSR = e_t[i]/4 + e_f[i]/1000.;}
if( ch[i] == ChMap::ScinL ) {hit.eSL = e[i]; hit.tSL = e_t[i]/4 + e_f[i]/1000.;}
if( ch[i] == ChMap::dFR ) {hit.eFR = e[i]; hit.tFR = e_t[i]/4 + e_f[i]/1000.;}
if( ch[i] == ChMap::dFL ) {hit.eFL = e[i]; hit.tFL = e_t[i]/4 + e_f[i]/1000.;}
if( ch[i] == ChMap::dBR ) {hit.eBR = e[i]; hit.tBR = e_t[i]/4 + e_f[i]/1000.;}
if( ch[i] == ChMap::dBL ) {hit.eBL = e[i]; hit.tBL = e_t[i]/4 + e_f[i]/1000.;}
if( ch[i] == ChMap::Cathode ) {hit.eCath = e[i]; hit.tCath = e_t[i]/4 + e_f[i]/1000.;}
if( ch[i] == ChMap::AnodeF ) {hit.eAF = e[i]; hit.tAF = e_t[i]/4 + e_f[i]/1000.;}
if( ch[i] == ChMap::AnodeB ) {hit.eAB = e[i]; hit.tAB = e_t[i]/4 + e_f[i]/1000.;}
if( ch[i] == ChMap::ScinR ) {hit.eSR = e[i]; hit.tSR = e_t[i] + e_f[i]/1000.*4;}
if( ch[i] == ChMap::ScinL ) {hit.eSL = e[i]; hit.tSL = e_t[i] + e_f[i]/1000.*4;}
if( ch[i] == ChMap::dFR ) {hit.eFR = e[i]; hit.tFR = e_t[i] + e_f[i]/1000.*4;}
if( ch[i] == ChMap::dFL ) {hit.eFL = e[i]; hit.tFL = e_t[i] + e_f[i]/1000.*4;}
if( ch[i] == ChMap::dBR ) {hit.eBR = e[i]; hit.tBR = e_t[i] + e_f[i]/1000.*4;}
if( ch[i] == ChMap::dBL ) {hit.eBL = e[i]; hit.tBL = e_t[i] + e_f[i]/1000.*4;}
if( ch[i] == ChMap::Cathode ) {hit.eCath = e[i]; hit.tCath = e_t[i] + e_f[i]/1000.*4;}
if( ch[i] == ChMap::AnodeF ) {hit.eAF = e[i]; hit.tAF = e_t[i] + e_f[i]/1000.*4;}
if( ch[i] == ChMap::AnodeB ) {hit.eAB = e[i]; hit.tAB = e_t[i] + e_f[i]/1000.*4;}
for( int j = i+1; j < sn.GetSize(); j++){
coin->Fill(ch[i], ch[j]);
@ -174,7 +178,7 @@ void SplitPolePlotter(TChain *tree){
// if( hit.eCath == 0 ) return kTRUE;
// if( hit.eCath > 13000 ) return kTRUE;
hit.CalData();
hit.CalData(4);
if( !TMath::IsNaN(hit.x1) || !TMath::IsNaN(hit.x2) ) {
hFocal->Fill(hit.x1, hit.x2);
@ -182,7 +186,8 @@ void SplitPolePlotter(TChain *tree){
hB->Fill(hit.x2);
hXavg->Fill(hit.xAvg);
hXavgVQ->Fill(hit.xAvg, dQ);
hXavg_Q->Fill(hit.xAvg, dQ);
hXavg_Theta->Fill( hit.xAvg, hit.theta);
for( int i = 0; i < 400; i++){
double y = -50 + 100/400.*i;
@ -215,11 +220,11 @@ void SplitPolePlotter(TChain *tree){
//^###########################################################
//^ * Plot
//^###########################################################
TCanvas * canvas = new TCanvas("cc", "Split-Pole", 1800, 1200);
TCanvas * canvas = new TCanvas("cc", "Split-Pole", 1600, 1200);
gStyle->SetOptStat("neiou");
canvas->Divide(3, 3);
canvas->Divide(4, 3);
canvas->cd(1); PID->Draw("colz");
//canvas->cd(2); coin->Draw("colz");
@ -228,7 +233,7 @@ void SplitPolePlotter(TChain *tree){
canvas->cd(3); hF->Draw();
canvas->cd(4); hB->Draw();
canvas->cd(5); hXavgVQ->Draw("colz");
canvas->cd(5); hXavg_Q->Draw("colz");
canvas->cd(6); hXavg->Draw("colz");
@ -238,4 +243,6 @@ void SplitPolePlotter(TChain *tree){
canvas->cd(9); canvas->cd(9)->SetLogy(); hMulti->Draw();
canvas->cd(10); hXavg_Theta->Draw("colz");
}

View File

@ -7,7 +7,8 @@ void script(){
TChain * chain = new TChain("tree");
chain->Add("data/temp_002_336_1000.root");
// chain->Add("data/temp_002_336_1000.root");
chain->Add("run123_100000.root");
SplitPolePlotter(chain);

View File

@ -1,184 +0,0 @@
#define splitpole_cxx
#include "splitpole.h"
#include <TH2.h>
#include <TH1.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include "../analyzers/SplitPoleHit.h"
namespace ChMap{
const short ScinR = 0;
const short ScinL = 1;
const short dFR = 8;
const short dFL = 9;
const short dBR = 10;
const short dBL = 11;
const short Cathode = 7;
const short AnodeF = 13;
const short AnodeB = 15;
};
const double c = 299.792458; // mm/ns
const double pi = M_PI;
const double deg2rad = pi/180.;
SplitPoleHit hit;
TH2F * PID;
TH2F * coin;
TH1F * hMulti;
TH1F * hF;
TH1F * hB;
TH1F * hXavg;
TH2F * hFocal;
TH2F * hXavgVQ;
TH2F * haha;
TH1F * hEx;
ULong64_t t1, t2;
#define XMIN -20
#define XMAX 100
void splitpole::Begin(TTree * /*tree*/){
TString option = GetOption();
PID = new TH2F("hPID", "PID; Scin_X ; AnodeB", 200, 0, 30000, 100, 0, 70000);
coin = new TH2F("hCoin", "Coincident ", 16, 0, 16, 16, 0, 16);
hMulti = new TH1F("hMulti", "Multiplicity", 16, 0, 16);
hF = new TH1F("hF", "Front delay line position", 600, XMIN, XMAX);
hB = new TH1F("hB", "Back delay line position", 600, XMIN, XMAX);
hXavg = new TH1F("hAvg", "Xavg", 600, XMIN, XMAX);
hFocal = new TH2F("hFocal", "Front vs Back ", 200, XMIN, XMAX, 200, XMIN, XMAX);
hXavgVQ = new TH2F("hXavgVQ", "Xavg vs Q ", 200, XMIN, XMAX, 200, 0, 40000);
haha = new TH2F("haha", "", 400, XMIN, XMAX, 400, -50, 50);
hEx = new TH1F("hEx", "Ex; Ex [MeV]; count/100 keV", 250, -5, 20);
hit.SetMassTablePath("../analyzers/mass20.txt");
hit.CalConstants("12C", "12C", "4He", 80, 5); // 80MeV, 5 deg
hit.CalZoffset(1.41); // 1.41 T
t1 = 0;
t2 = 0;
}
Bool_t splitpole::Process(Long64_t entry){
b_multi->GetEntry(entry);
b_ch->GetEntry(entry);
b_e->GetEntry(entry);
b_e2->GetEntry(entry);
b_e_t->GetEntry(entry);
b_e_f->GetEntry(entry);
//if( multi < 9) return kTRUE;
hit.ClearData();
hMulti->Fill(multi);
for( int i = 0; i < multi; i++){
t2 = e_t[i];
if( t2 < t1 ) printf("entry %lld-%d, timestamp is not in order. %llu, %llu\n", entry, i, t2, t1);
if( i == 0 ) t1 = e_t[i];
if( ch[i] == ChMap::ScinR ) {hit.eSR = e[i]; hit.tSR = e_t[i] + e_f[i]/1000.;}
if( ch[i] == ChMap::ScinL ) {hit.eSL = e[i]; hit.tSL = e_t[i] + e_f[i]/1000.;}
if( ch[i] == ChMap::dFR ) {hit.eFR = e[i]; hit.tFR = e_t[i] + e_f[i]/1000.;}
if( ch[i] == ChMap::dFL ) {hit.eFL = e[i]; hit.tFL = e_t[i] + e_f[i]/1000.;}
if( ch[i] == ChMap::dBR ) {hit.eBR = e[i]; hit.tBR = e_t[i] + e_f[i]/1000.;}
if( ch[i] == ChMap::dBL ) {hit.eBL = e[i]; hit.tBL = e_t[i] + e_f[i]/1000.;}
if( ch[i] == ChMap::Cathode ) {hit.eCath = e[i]; hit.tCath = e_t[i] + e_f[i]/1000.;}
if( ch[i] == ChMap::AnodeF ) {hit.eAF = e[i]; hit.tAF = e_t[i] + e_f[i]/1000.;}
if( ch[i] == ChMap::AnodeB ) {hit.eAB = e[i]; hit.tAB = e_t[i] + e_f[i]/1000.;}
for( int j = i+1; j < multi; j++){
coin->Fill(ch[i], ch[j]);
}
}
unsigned int dQ = hit.eAB; // delta Q
unsigned int Qt = hit.eSR; // total Q
if( Qt > 0 && dQ > 0 ) {
PID->Fill(Qt, dQ);
}
// if( hit.eAF < 50000 ) return kTRUE;
// if( hit.eCath == 0 ) return kTRUE;
// if( hit.eCath > 13000 ) return kTRUE;
hit.CalData();
if( !TMath::IsNaN(hit.x1) || !TMath::IsNaN(hit.x2) ) {
hFocal->Fill(hit.x1, hit.x2);
hF->Fill(hit.x1);
hB->Fill(hit.x2);
hXavg->Fill(hit.xAvg);
hXavgVQ->Fill(hit.xAvg, dQ);
for( int i = 0; i < 400; i++){
double y = -50 + 100/400.*i;
double x = (y/42.8625 + 0.5)* ( hit.x2-hit.x1) + hit.x1;
haha->Fill(x,y);
}
double ex = hit.Rho2Ex( (hit.xAvg/100 + 0.363) );
//if( XMIN < hit.xAvg && hit.xAvg < XMAX) printf("x1 : %6.2f, x2 : %6.2f, xAvg %6.2f cm , ex : %f \n", hit.x1, hit.x2, hit.xAvg, ex);
hEx->Fill(ex);
}
return kTRUE;
}
void splitpole::Terminate(){
TCanvas * canvas = new TCanvas("cc", "Split-Pole", 1800, 1200);
gStyle->SetOptStat("neiou");
canvas->Divide(3, 3);
canvas->cd(1); PID->Draw("colz");
//canvas->cd(2); coin->Draw("colz");
canvas->cd(2); haha->Draw("colz");
canvas->cd(3); hF->Draw();
canvas->cd(4); hB->Draw();
canvas->cd(5); hXavgVQ->Draw("colz");
canvas->cd(6); hXavg->Draw("colz");
canvas->cd(7); hEx->Draw();
canvas->cd(8); coin->Draw("colz");
canvas->cd(9); canvas->cd(9)->SetLogy(); hMulti->Draw();
}

View File

@ -1,101 +0,0 @@
//////////////////////////////////////////////////////////
// This class has been automatically generated on
// Wed Jun 14 15:51:03 2023 by ROOT version 6.26/04
// from TTree tree/temp_036.root
// found on file: temp_036.root
//////////////////////////////////////////////////////////
#ifndef splitpole_h
#define splitpole_h
#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
#include <TSelector.h>
// Header file for the classes stored in the TTree if any.
#define MaxMulti 100
class splitpole : public TSelector {
public :
TTree *fChain; //!pointer to the analyzed TTree or TChain
// Fixed size dimensions of array or collections stored in the TTree if any.
// Declaration of leaf types
ULong64_t evID;
UShort_t multi;
UShort_t bd[MaxMulti]; //[multi]
UShort_t ch[MaxMulti]; //[multi]
UShort_t e[MaxMulti]; //[multi]
UShort_t e2[MaxMulti]; //[multi]
ULong64_t e_t[MaxMulti]; //[multi]
UShort_t e_f[MaxMulti]; //[multi]
// List of branches
TBranch *b_event_ID; //!
TBranch *b_multi; //!
TBranch *b_bd; //!
TBranch *b_ch; //!
TBranch *b_e; //!
TBranch *b_e2; //!
TBranch *b_e_t; //!
TBranch *b_e_f; //!
splitpole(TTree * /*tree*/ =0) : fChain(0) { }
virtual ~splitpole() { }
virtual Int_t Version() const { return 2; }
virtual void Begin(TTree *tree);
virtual void SlaveBegin(TTree *tree);
virtual void Init(TTree *tree);
virtual Bool_t Notify();
virtual Bool_t Process(Long64_t entry);
virtual Int_t GetEntry(Long64_t entry, Int_t getall = 0) { return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; }
virtual void SetOption(const char *option) { fOption = option; }
virtual void SetObject(TObject *obj) { fObject = obj; }
virtual void SetInputList(TList *input) { fInput = input; }
virtual TList *GetOutputList() const { return fOutput; }
virtual void SlaveTerminate();
virtual void Terminate();
ClassDef(splitpole,0);
};
#endif
#ifdef splitpole_cxx
void splitpole::Init(TTree *tree){
// Set branch addresses and branch pointers
if (!tree) return;
fChain = tree;
fChain->SetMakeClass(1);
fChain->SetBranchAddress("evID", &evID, &b_event_ID);
fChain->SetBranchAddress("multi", &multi, &b_multi);
fChain->SetBranchAddress("bd", bd, &b_bd);
fChain->SetBranchAddress("ch", ch, &b_ch);
fChain->SetBranchAddress("e", e, &b_e);
fChain->SetBranchAddress("e2", e2, &b_e2);
fChain->SetBranchAddress("e_t", e_t, &b_e_t);
fChain->SetBranchAddress("e_f", e_f, &b_e_f);
}
Bool_t splitpole::Notify(){
return kTRUE;
}
void splitpole::SlaveBegin(TTree * /*tree*/){
TString option = GetOption();
}
void splitpole::SlaveTerminate(){
}
#endif // #ifdef splitpole_cxx

View File

@ -194,14 +194,14 @@ public:
isConstantCal = false;
}
void CalData(){
void CalData(float scale = 2.){
if( eSR > 0 && eSL > 0 ) eSAvg = (eSR + eSL)/2;
if( eSR > 0 && eSL == 0 ) eSAvg = eSR;
if( eSR == 0 && eSL > 0 ) eSAvg = eSL;
if( tFR > 0 && tFL > 0 ) x1 = (tFL - tFR)/2./2.1;
if( tBR > 0 && tBL > 0 ) x2 = (tBL - tBR)/2./1.98;
if( tFR > 0 && tFL > 0 ) x1 = (tFL - tFR)/scale/2.1;
if( tBR > 0 && tBL > 0 ) x2 = (tBL - tBR)/scale/1.98;
if( !std::isnan(x1) && !std::isnan(x2)) {