From 4c7938ed6dd7929238b36ec499ba687c4a00111b Mon Sep 17 00:00:00 2001 From: Thomas Prokscha Date: Fri, 19 Dec 2008 09:53:10 +0000 Subject: [PATCH] Added analysis.C root macro for analyzing geant4 output files. --- geant4/LEMuSR/run/data/analysis.C | 228 ++++++++++++++++++++++++++++++ 1 file changed, 228 insertions(+) create mode 100644 geant4/LEMuSR/run/data/analysis.C diff --git a/geant4/LEMuSR/run/data/analysis.C b/geant4/LEMuSR/run/data/analysis.C new file mode 100644 index 0000000..fca0e41 --- /dev/null +++ b/geant4/LEMuSR/run/data/analysis.C @@ -0,0 +1,228 @@ +#define analysis_cxx +#include "analysis.h" +#include +#include +#include + +void analysis::Loop() +{ +// In a ROOT session, you can do: +// Root > .L analysis.C +// Root > analysis t +// Root > t.GetEntry(12); // Fill t data members with entry number 12 +// Root > t.Show(); // Show values of entry 12 +// Root > t.Show(16); // Read and show values of entry 16 +// Root > t.Loop(); // Loop on all entries <- Very Important! +// Root > .ls // Shows what histos are available for plotting +// +// This is the loop skeleton where: +// jentry is the global entry number in the chain +// ientry is the entry number in the current Tree +// Note that the argument to GetEntry must be: +// jentry for TChain::GetEntry +// ientry for TTree::GetEntry and TBranch::GetEntry +// +// To read only selected branches, insert statements as: +// METHOD1: +// fChain->SetBranchStatus("*",0); // disable all branches +// fChain->SetBranchStatus("branchname",1); // activate branchname +// METHOD2: replace line +// fChain->GetEntry(jentry); //read all branches +// b_branchname->GetEntry(ientry); //read only this branch +// +// +// To create TEMPLATE files for analysis proceed as follows: +// 1. Load the root file containing all the variables, e.g.: +// TFile* f=new TFile("data/lem4_1049.root"); +// 2. Quickly check the Tree and Variable names: +// f->ls(); // Shows tree name (i.e. t1) +// t1->Print(); // Shows variable names +// 3. Create the template files .C and .h (t1 is the Tree name) +// t1->MakeClass("myAnalysis") +// +// +// How to write data to an EXTERNAL file (use cout instead of outdata for printing on screen): +// See also: www.bgsu.edu/departments/compsci/docs/index.html +// +// ofstream outdata; +// outdata.open("export.dat"); +// for (Int_t i=0; iGetNbinsX(); i++) {outdata<GetBinContent(i)<Sumw2(); + hEdepositCF->Sumw2(); + hEdepTrig->Sumw2(); + + Long64_t nentries = fChain->GetEntriesFast(); + //nentries=1000; + + Long64_t nbytes = 0, nb = 0; + for (Long64_t jentry=0; jentryGetEntry(jentry); nbytes += nb; + + // FILL IN HISTOGRAMS + hEdepositCF->Fill(save_ke[0]/1000.); + hBeamSpot->Fill(save_x[0], save_y[0]); + + //for (Int_t i=0; iFill(save_ke[i]);} + // } + + hTof->Fill(muTargetTime); + if (muTargetTime>0.23) tofFlag = 1; + //tofFlag = 1.; + + // Coincidence between Up scintillators + for (Int_t i=0; iFill(det_edep[i]); + for (Int_t j=0; jeCut && det_edep[j]>eCut && tofFlag){ + hDettUp->Fill(det_time_start[j]); + hDetz->Fill(det_z[i]); + hEventID->Fill(eventID); + break; //fill only once + } + } + } + } + } + + // Coincidence between Down scintillators + for (Int_t i=0; iFill(det_edep[i]); + for (Int_t j=0; jeCut && det_edep[j]>eCut && tofFlag){ + hDettDown->Fill(det_time_start[j]); + hDetz->Fill(det_z[i]); + hEventID->Fill(eventID); + break; + } + } + } + } + } + + // Coincidence between Right scintillators + for (Int_t i=0; iFill(det_edep[i]); + for (Int_t j=0; jeCut && det_edep[j]>eCut && tofFlag){ + hDettRight->Fill(det_time_start[j]); + hDetz->Fill(det_z[i]); + hEventID->Fill(eventID); + break; + } + } + } + } + } + + // Coincidence between Left scintillators + for (Int_t i=0; iFill(det_edep[i]); + for (Int_t j=0; jeCut && det_edep[j]>eCut && tofFlag){ + hDettLeft->Fill(det_time_start[j]); + hDetz->Fill(det_z[i]); + hEventID->Fill(eventID); + break; + } + } + } + } + } + + + // Testing event coincidence with SciFi + Bool_t twentyone= false; + Bool_t twozero7 = false; + Bool_t twozero8 = false; + Int_t ikk=0; + for (Int_t i=0; i0.0)) {twentyone = true; ikk = i;} + if ((det_ID[i]==207)&&(det_edep[i]>0.0)) twozero7 = true; + if ((det_ID[i]==208)&&(det_edep[i]>0.0)) twozero8 = true; + } + //if (twentyone) hEdepoTest->Fill(det_edep[ikk]); // Fiber hits + //if (twozero7&& twozero8) hEdepoTest->Fill(det_edep[ikk]); // Coincidence + if (twentyone&&twozero7&& twozero8) hEdepoTest->Fill(det_edep[ikk]); // Fibre hits + coincidence + + } + + // PLOT HISTOGRAMS + TCanvas* c1= new TCanvas("c1","canvas 1"); + + // hEdeposited->Scale(1/(number_of_generated_events*hEdeposited->GetBinWidth(1))); + TH1F *htemp1 = (TH1F*) hDettUp->Clone(); htemp1->SetName("htemp1"); + TH1F *htemp2 = (TH1F*) hDettUp->Clone(); htemp2->SetName("htemp2"); + + TH1F *htemp3 = (TH1F*) hDettUp->Clone(); htemp3->SetName("htemp3"); + TH1F *htemp4 = (TH1F*) hDettUp->Clone(); htemp4->SetName("htemp4"); + + TH1F *htemp5 = (TH1F*) hDettUp->Clone(); htemp5->SetName("htemp5"); + htemp5->SetTitle("Muon decay asymmetry U-D; Time (#mus); Asymmetry"); + htemp5->Sumw2(); + + TH1F *htemp6 = (TH1F*) hDettUp->Clone(); htemp6->SetName("htemp6"); + htemp6->SetTitle("Muon decay asymmetry L-R; Time (#mus); Asymmetry"); + htemp6->Sumw2(); + + htemp1->Add(hDettUp,hDettDown,1.,-1.); + htemp2->Add(hDettUp,hDettDown,1., 1.); + htemp5->Divide(htemp1,htemp2,1.,1.); + + htemp3->Add(hDettLeft,hDettRight,1.,-1.); + htemp4->Add(hDettLeft,hDettRight,1., 1.); + htemp6->Divide(htemp3,htemp4,1.,1.); + + htemp5->SetLineWidth(2); + htemp5->SetLineColor(2); + htemp5->SetAxisRange(-0.8,0.8,"Y"); + htemp5->Draw("hist"); + + htemp6->SetLineWidth(2); + htemp6->SetLineColor(4); + htemp6->SetAxisRange(-0.8,0.8,"Y"); + htemp6->Draw("same"); + + hEdepositCF->SetLineWidth(2); + hEdepositCF->SetLineColor(2); + hEdepositCF->Draw("hist"); + //hEdepositCF->Fit("gaus"); + + TCanvas* c2= new TCanvas("c2","canvas 2"); + hEdeposited->Draw(); + + // sprintf(filenamePrint,"musrSrClass_%s.ps",fileNr); + // c1->Print(filenamePrint); +}