#define e121_analysis_cxx
// INPUT FILES
#define INPUT1 "/data.local1/E121/root-files/part1/e121_run_0253.root"
//#define INPUT2 ""
//#define INPUT3 ""
//#define INPUT4 ""
//#define INPUT5 ""
//#define INPUT6 ""
//#define INPUT7 ""
//#define INPUT8 ""
//#define INPUT9 ""
//#define INPUT10 ""
//OUTPUT FILE
#define OUTPUT "run253.ana.root"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <cmath>
#include <string>
#include <cstdlib>
#include <cstdio>
#include <ctime>
#include <fstream>
#include <iostream>
#include <stdint.h>
#include "TROOT.h"
#include "TAttText.h"
#include "TAxis.h"
#include "TCanvas.h"
#include "TChain.h"
#include "TCut.h"
#include "TF1.h"
#include "TFile.h"
#include "TGraph.h"
#include "TGraphAsymmErrors.h"
#include "TGraphErrors.h"
#include "TH1.h"
#include "THistPainter.h"
#include "TKey.h"
#include "TLatex.h"
#include "TLegend.h"
#include "TMath.h"
#include "TMatrixD.h"
#include "TMinuit.h"
#include "TMultiGraph.h"
#include "TNtuple.h"
#include "TPave.h"
#include "TPaveText.h"
#include "TPoint.h"
#include "TRandom.h"
#include "TRint.h"
#include "TString.h"
#include "TTree.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TSystem.h"
#include "TProfile.h"
#include "TVirtualFitter.h"
#include "TLegend.h"
#include "TColor.h"
#include "TBranch.h"
#include "TList.h"
using namespace std;
void loop(TChain *fChain)
{
if (fChain == 0) return;
//#include "e121_raw.h" // use this for raw mapping of channels
#include "e121_mapped.h" // use this for mapped channel names
///////////////////////////
// histogram definitions //
///////////////////////////
// energy
TList *e_histo= new TList();
TH1D *h_e_pad_1[7];
TH1D *h_e_pad_2[7];
TH1D *h_e_pad_3[7];
TH1D *h_e_pad_4[7];
TH1D *h_e_pad_5[7];
TH1D *h_e_pad_6[7];
TH1D *h_e_pad_n[7];
TH1D *h_e_pad_n_lo[7];
TH1D *h_e_pad_n_hi[7];
char chnumber[300];
for (int i=0;i<7;i++)
{
sprintf(chnumber,"e_pad_1_strip_%d",i);
h_e_pad_1[i]=new TH1D(chnumber,chnumber,4096,0.5,4096.5);
e_histo->Add(h_e_pad_1[i]);
sprintf(chnumber,"e_pad_2_strip_%d",i);
h_e_pad_2[i]=new TH1D(chnumber,chnumber,4096,0.5,4096.5);
e_histo->Add(h_e_pad_2[i]);
sprintf(chnumber,"e_pad_3_strip_%d",i);
h_e_pad_3[i]=new TH1D(chnumber,chnumber,4096,0.5,4096.5);
e_histo->Add(h_e_pad_3[i]);
sprintf(chnumber,"e_pad_4_strip_%d",i);
h_e_pad_4[i]=new TH1D(chnumber,chnumber,4096,0.5,4096.5);
e_histo->Add(h_e_pad_4[i]);
sprintf(chnumber,"e_pad_5_strip_%d",i);
h_e_pad_5[i]=new TH1D(chnumber,chnumber,4096,0.5,4096.5);
e_histo->Add(h_e_pad_5[i]);
sprintf(chnumber,"e_pad_6_strip_%d",i);
h_e_pad_6[i]=new TH1D(chnumber,chnumber,4096,0.5,4096.5);
e_histo->Add(h_e_pad_6[i]);
sprintf(chnumber,"e_pad_%d_n-side",i);
h_e_pad_n[i]=new TH1D(chnumber,chnumber,4096,0.5,4096.5);
e_histo->Add(h_e_pad_n[i]);
sprintf(chnumber,"e_pad_%d_n-side_lo",i);
h_e_pad_n_lo[i]=new TH1D(chnumber,chnumber,4096,0.5,4096.5);
e_histo->Add(h_e_pad_n_lo[i]);
sprintf(chnumber,"e_pad_%d_n-side_hi",i);
h_e_pad_n_hi[i]=new TH1D(chnumber,chnumber,4096,0.5,4096.5);
e_histo->Add(h_e_pad_n_hi[i]);
}
TH1D *h_e_dssd_left=new TH1D("e_dssd_left","e_dssd_left",4096,0.5,4096.5);
e_histo->Add(h_e_dssd_left);
TH1D *h_e_dssd_right=new TH1D("e_dssd_right","e_dssd_right",4096,0.5,4096.5);
e_histo->Add(h_e_dssd_right);
TH1D *h_e_dssd_top=new TH1D("e_dssd_top","e_dssd_top",4096,0.5,4096.5);
e_histo->Add(h_e_dssd_top);
TH1D *h_e_dssd_bottom=new TH1D("e_dssd_bottom","e_dssd_bottom",4096,0.5,4096.5);
e_histo->Add(h_e_dssd_bottom);
TH1D *h_e_csi[2];
h_e_csi[0]=new TH1D("e_csi_0","e_csi_0",4096,0.5,4096.5);
e_histo->Add(h_e_csi[0]);
h_e_csi[1]=new TH1D("e_csi_1","e_csi_1",4096,0.5,4096.5);
e_histo->Add(h_e_csi[1]);
// position
// time
// scaler
TList *sc_histo= new TList();
Int_t sc_bins = 500;
TGraph *g_sc_si=new TGraph();
g_sc_si->SetNameTitle("g_sc_si", "rate Si OR");
sc_histo->Add(g_sc_si);
TGraph *g_sc_mwpc=new TGraph();
g_sc_mwpc->SetNameTitle("g_sc_mwpc", "rate MWPC anode");
sc_histo->Add(g_sc_mwpc);
TGraph *g_sc_target=new TGraph();
g_sc_target->SetNameTitle("g_sc_target", "target density");
sc_histo->Add(g_sc_target);
TGraph *g_sc_trafo=new TGraph();
g_sc_trafo->SetNameTitle("g_sc_trafo", "ESR ion current");
sc_histo->Add(g_sc_trafo);
TGraph *g_sc_inhibit=new TGraph();
g_sc_inhibit->SetNameTitle("g_sc_inhibit", "target inhibit");
sc_histo->Add(g_sc_inhibit);
TGraph *g_sc_clock=new TGraph();
g_sc_clock->SetNameTitle("g_sc_clock", "1.5 MHz clock");
sc_histo->Add(g_sc_clock);
ULong64_t sc_si_sum = 0;
ULong64_t sc_mwpc_sum = 0;
ULong64_t sc_clock_sum = 0;
ULong64_t sc_target_sum = 0;
ULong64_t sc_trafo_sum = 0;
ULong64_t sc_inhibit_sum = 0;
Int_t sc_interval = 1; //interval time [sec] for scaler plotting
Double_t clock_ref = 1.5e6; // clock reference rate = 1.5 MHz
Double_t t_diff = 0;
Double_t t_elapsed = 0;
// get number of events to process
Long64_t nentries = fChain->GetEntries();
Long64_t nbytes = 0, nb = 0;
//////////////////////
// EVENT LOOP START //
//////////////////////
for (Long64_t i=0; i<nentries;i++)
// for (Long64_t i=0; i<10000;i++)
{
nb = fChain->GetEntry(i); nbytes += nb;
// event countdown
if ((float(i)/100000.)==int(i/100000)){cout << "event: " << i << " \tof " << nentries << endl;}
//////////////////////
//scaler processing //
//////////////////////
sc_si_sum += sc_silicon_or;
sc_mwpc_sum += sc_mwpc_anode;
sc_clock_sum += sc_clock;
sc_target_sum += sc_target;
sc_trafo_sum += sc_trafo;
sc_inhibit_sum += sc_inhibit;
t_elapsed += (double)(sc_clock/clock_ref);
t_diff = (double)(sc_clock_sum/clock_ref);
// fill the scaler sums to histo for very interval
if( !((int)t_diff % sc_interval) && (int)t_diff > 0 )
//if( true )
{
if ( (int)t_diff <= 2*sc_interval ) //exclude large time intervals >> crap data
{
// cout << "time: " << t_elapsed << " trafo:" << sc_trafo_sum;
// cout << " diff: " << t_diff << " modulo: " << ((int)t_diff % sc_interval) << endl;
g_sc_si->SetPoint(g_sc_si->GetN() ,t_elapsed, sc_si_sum);
g_sc_mwpc->SetPoint(g_sc_mwpc->GetN() ,t_elapsed, sc_mwpc_sum);
g_sc_target->SetPoint(g_sc_target->GetN() ,t_elapsed, sc_target_sum);
g_sc_trafo->SetPoint(g_sc_trafo->GetN() ,t_elapsed, sc_trafo_sum);
g_sc_inhibit->SetPoint(g_sc_inhibit->GetN(),t_elapsed, sc_inhibit_sum);
g_sc_clock->SetPoint(g_sc_clock->GetN() ,t_elapsed, sc_clock_sum);
}
sc_si_sum = 0;
sc_mwpc_sum = 0;
sc_clock_sum = 0;
sc_target_sum = 0;
sc_trafo_sum = 0;
sc_inhibit_sum = 0;
}
////////////////////
// energy spectra //
////////////////////
if(TRIGGER==1)
{
for (int i=0;i<7;i++)
{
h_e_pad_1[i]->Fill(e_pad_1[i]);
h_e_pad_2[i]->Fill(e_pad_2[i]);
h_e_pad_3[i]->Fill(e_pad_3[i]);
h_e_pad_4[i]->Fill(e_pad_4[i]);
h_e_pad_5[i]->Fill(e_pad_5[i]);
h_e_pad_6[i]->Fill(e_pad_6[i]);
if ( i > 0 )
{
h_e_pad_n[i]->Fill(e_pad_n[i]);
if ( t_elapsed >= 30 && t_elapsed <= 100 ) h_e_pad_n_hi[i]->Fill(e_pad_n[i]);
if ( t_elapsed >= 250 && t_elapsed <= 350 ) h_e_pad_n_lo[i]->Fill(e_pad_n[i]);
}
}
}
}
////////////////////
// EVENT LOOP END //
////////////////////
// write output file
TFile *outfile = TFile::Open(OUTPUT, "RECREATE");
e_histo->Write("energy", TObject::kSingleKey);
sc_histo->Write("scaler", TObject::kSingleKey);
outfile->Close();
cout << "\033[0;32m" << OUTPUT << " is created!\033[0m" << endl;
// direct plotting (optional)
TCanvas *c_test = new TCanvas("test","test", 1000, 1200);
c_test->Divide(1,5);
c_test->cd(1);
g_sc_si->SetFillColor(15);
g_sc_si->Draw("AB");
c_test->cd(2);
g_sc_mwpc->SetFillColor(15);
g_sc_mwpc->Draw("AB");
c_test->cd(3);
g_sc_target->SetFillColor(15);
g_sc_target->Draw("AB");
c_test->cd(4);
g_sc_trafo->SetFillColor(15);
g_sc_trafo->Draw("AB");
c_test->cd(5);
... 108 more lines ...
|