WCSim/sample-root-scripts/plot_pmts.C at develop · WCSim/WCSim · GitHub

History

335 lines (306 loc) · 9.72 KB

Raw

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

#include <iostream>

#include <stdio.h>

#include <stdlib.h>

#include <algorithm>

#include "TTree.h"

#include "TGraph.h"

#include "TEllipse.h"

#include "TH1F.h"

#include "THStack.h"

#include "TH1D.h"

#include "TStyle.h"

#include "TROOT.h"

#include "TSystem.h"

#include "TCanvas.h"

#include "TFile.h"

#include "WCSimRootOptions.hh"

#include "WCSimRootGeom.hh"

#include "WCSimRootEvent.hh"

Produces X-Y and R-Z displays for each PMT type

- In the X-Y displays, the detector edge is shown as a circle

- Note that both top & bottom cap PMTs are displayed on this plot

By construction in WCSim, PMTs are placed in the same location on each end cap

- In the R-Z displays, the detector edge is the edge of the image

Output

- display_geo.pdf has the 2D plots

- There is one page per PMT type (20" / OD / mPMT)

- display_geo.root has

- The TGraph's used to create the 2D plots

- Additionally two THStack's (one for R, one for Z), that contain the 1D information

Running

From the command line, something like

root -b -q $WCSIM_SOURCE_DIR/sample-root-scripts/plot_pmts.C+g'(1,"wcsim.root")'

or however you prefer to run your root macros

int display(int verbose, //how verbose should the printout be? 0 = lowest, higher number = more

const char *filename); //WCSim file. Should work for all geometries

TGraph * MakeGraph(vector<pair<double, double> > v, int col, int style, const float size=1.5)

{

cout << "Making graph from vector of size " << v.size() << endl;

TGraph * g = new TGraph(v.size());

for(size_t i = 0; i < v.size(); i++) {

auto pair = v.at(i);

g->SetPoint(i, pair.first, pair.second);

//cout << pair.first << "\t" << pair.second << endl;

}

g->SetMarkerColor(col);

g->SetMarkerStyle(style);

g->SetMarkerSize(size);

return g;

}

double GetAzimuth(double x, double y)

{

double t = TMath::ATan(y / x);

if(x < 0) {

if(y < 0)

t -= TMath::Pi();

else

t += TMath::Pi();

}

t = TMath::ATan2(y, x);

return t;

}

void PrintPMT(const WCSimRootPMT * tube, const int ipmt)

{

double x, y, z, r, t;

x = tube->GetPosition(0);

y = tube->GetPosition(1);

z = tube->GetPosition(2);

r = TMath::Sqrt(x*x + y*y);

t = GetAzimuth(x, y);

cout << std::setw(6) << ipmt

<< " " << std::setw(8) << tube->GetTubeNo()

<< " " << std::setw(8) << tube->GetmPMTNo()

<< " " << std::setw(3) << tube->GetmPMT_PMTNo()

<< " " << std::setw(9) << x

<< " " << std::setw(9) << y

<< " " << std::setw(9) << z

<< " " << std::setw(9) << r

<< " " << std::setw(9) << t

<< " " << std::setw(2) << tube->GetCylLoc()

<< endl;

}

std::string CycLocs[9] = {"ID top cap",

"ID barrel",

"ID bottom cap",

"OD bottom cap",

"OD barrel",

"OD top cap",

"mPMT top cap",

"mPMT barrel",

"mPMT bottom cap"};

void DumpGeoTree(WCSimRootGeom * geo,

TCanvas * c,

int verbose)

{

const float detR = geo->GetWCCylRadius();

const float detZ = geo->GetWCCylLength() / 2;

TClonesArray * pmts = geo->GetPMTs();

cout << "N 20\" PMTs in TClonesArray: " << pmts->GetEntries() << endl

<< "N 20\" PMTs from rootgeom : " << geo->GetWCNumPMT() << endl;

TClonesArray * pmts2 = geo->GetPMTs(true);

cout << "N mPMTs in TClonesArray: " << pmts2->GetEntries() << endl

<< "N mPMTs from rootgeom : " << geo->GetWCNumPMT(true) << endl;

TClonesArray * pmtsod = geo->GetODPMTs();

cout << "N OD PMTs in TClonesArray: " << pmtsod->GetEntries() << endl;

cout << "N OD PMTs from rootgeom : " << geo->GetODWCNumPMT() << endl;

const int narrays = 9;

vector<pair<double, double> > pos[narrays];

TH1D * Rs[narrays], * Zs[narrays];

vector<double> Rsv[narrays], Zsv[narrays];

THStack * stackR = new THStack("stackR", ";R (cm);Number in bin");

THStack * stackZ = new THStack("stackZ", ";Z (cm);Number in bin");

for(int i = 0; i < narrays; i++) {

Rs[i] = new TH1D(TString::Format("R%d", i), ";R (cm);Number in bin", TMath::Ceil(detR*8), 0, detR);

Zs[i] = new TH1D(TString::Format("Z%d", i), ";Z (cm);Number in bin", TMath::Ceil(detZ*8), -detZ, +detZ);

stackR->Add(Rs[i]);

stackZ->Add(Zs[i]);

}

double x, y, z, r, t;

//Get the 20" PMT information

if(verbose > 1)

cout << endl << "20\" PMT list" << endl;

for(int ipmt = 0; ipmt < pmts->GetEntries() - 1; ipmt++) {

const WCSimRootPMT * tube = geo->GetPMTPtr(ipmt);

//get the tube position. Depends on cycloc

x = tube->GetPosition(0);

y = tube->GetPosition(1);

z = tube->GetPosition(2);

r = TMath::Sqrt(x*x + y*y);

t = GetAzimuth(x, y);

const int loc = tube->GetCylLoc();

switch(loc) {

//endcap

case 0:

case 2:

pos[loc].push_back(std::make_pair(x, y));

break;

//wall

case 1:

pos[loc].push_back(std::make_pair(t, z));

break;

default:

cout << "Unknown case " << tube->GetCylLoc() << endl;

break;

}

Rs[loc]->Fill(r);

Zs[loc]->Fill(z);

Rsv[loc].push_back(r);

Zsv[loc].push_back(z);

if(verbose > 1)

PrintPMT(tube, ipmt);

}//20" PMTs

//and for OD pmts

if(verbose > 1)

cout << endl << "OD PMT list" << endl;

for(int ipmt = 0; ipmt < pmtsod->GetEntries() - 1; ipmt++) {

const WCSimRootPMT * tube = geo->GetODPMTPtr(ipmt);

//get the tube position. Depends on cycloc

x = tube->GetPosition(0);

y = tube->GetPosition(1);

z = tube->GetPosition(2);

r = TMath::Sqrt(x*x + y*y);

t = GetAzimuth(x, y);

const int loc = tube->GetCylLoc();

switch(loc) {

//endcap

case 3:

case 5:

pos[loc].push_back(std::make_pair(x, y));

break;

//wall

case 4:

pos[loc].push_back(std::make_pair(t, z));

break;

default:

cout << "Unknown case " << tube->GetCylLoc() << endl;

break;

}

Rs[loc]->Fill(r);

Zs[loc]->Fill(z);

Rsv[loc].push_back(r);

Zsv[loc].push_back(z);

if(verbose > 1)

PrintPMT(tube, ipmt);

}//OD PMTs

//and for mPMTs

if(verbose > 1)

cout << endl << "mPMT list" << endl;

for(int ipmt = 0; ipmt < pmts2->GetEntries() - 1; ipmt++) {

const WCSimRootPMT * tube = geo->GetPMTPtr(ipmt, true);

//get the tube position. Depends on cycloc

x = tube->GetPosition(0);

y = tube->GetPosition(1);

z = tube->GetPosition(2);

r = TMath::Sqrt(x*x + y*y);

t = GetAzimuth(x, y);

const int loc = tube->GetCylLoc();

switch(loc) {

//endcap

case 6:

case 8:

pos[loc].push_back(std::make_pair(x, y));

break;

//wall

case 7:

pos[loc].push_back(std::make_pair(t, z));

break;

default:

cout << "Unknown case " << tube->GetCylLoc() << endl;

break;

}

Rs[loc]->Fill(r);

Zs[loc]->Fill(z);

Rsv[loc].push_back(r);

Zsv[loc].push_back(z);

if(verbose > 1)

PrintPMT(tube, ipmt);

}//mPMTs

//now setup the plotting

c->SaveAs("display_geo.pdf[");

TGraph * g[narrays];

Color_t cols[narrays] = {kBlack, kBlue, kRed,

kBlack, kBlue, kRed,

kBlack, kBlue, kRed};

Style_t styles[narrays] = {29, 29, 29,

20, 20, 20,

21, 21, 21};

for(int i = 0; i < narrays; i++) {

cout << "Vector " << i << " (" << CycLocs[i]

<< ") contains " << pos[i].size() << " entries" << endl;

g[i] = MakeGraph(pos[i], cols[i], styles[i], 0.2);

}

TEllipse circle(0, 0, detR);

for(int i = 0; i < narrays; i+=3) {

c->cd(1)->DrawFrame(-detR,-detR,+detR,+detR,";X (cm);Y (cm)")->GetYaxis()->SetTitleOffset(1);

circle.Draw();

if(g[i]->GetN())

g[i]->Draw("P");

if(g[i+2]->GetN())

g[i+2]->Draw("P");

c->cd(2)->DrawFrame(-TMath::Pi(),-detZ,+TMath::Pi(),+detZ,";Azimuthal angle;Z (cm)")->GetYaxis()->SetTitleOffset(1);

if(g[i+1]->GetN())

g[i+1]->Draw("P");

c->SaveAs("display_geo.pdf");

}

c->SaveAs("display_geo.pdf]");

TFile fout("display_geo.root", "RECREATE");

for(int i = 0; i < narrays; i++)

g[i]->Write(TString::Format("graph%d", i));

stackR->Write();

stackZ->Write();

fout.Close();

for(int i = 0; i < narrays; i++) {

cout << i << " " << CycLocs[i] << endl;

if(!Rsv[i].size()) {

cout << " No PMTs in this cycloc" << endl;

continue;

}

cout << " MIN R: " << *std::min_element(Rsv[i].begin(), Rsv[i].end())

<< "\t MAX R: " << *std::max_element(Rsv[i].begin(), Rsv[i].end()) << endl;

cout << " MIN Z: " << *std::min_element(Zsv[i].begin(), Zsv[i].end())

<< "\t MAX Z: " << *std::max_element(Zsv[i].begin(), Zsv[i].end()) << endl;

}

cout << "Detector half z: " << detZ << endl

<< "Detector radius: " << detR << endl;

}

/////////////////////

int plot_pmts(int verbose=0,

const char *filename="wcsim.root")

{

// Clear global scope

//gROOT->Reset();

// Open the file

TFile * file = new TFile(filename,"read");

if (!file->IsOpen()){

cout << "Error, could not open input file: " << filename << endl;

return -1;

}

// Geometry tree - only need 1 "event"

TTree *geotree = (TTree*)file->Get("wcsimGeoT");

WCSimRootGeom *geo = 0;

geotree->SetBranchAddress("wcsimrootgeom", &geo);

if(verbose) std::cout << "Geotree has: " << geotree->GetEntries() << " entries (1 expected)" << std::endl;

if (geotree->GetEntries() == 0) {

exit(9);

}

geotree->GetEntry(0);

// Options tree - only need 1 "event"

TTree *opttree = (TTree*)file->Get("wcsimRootOptionsT");

WCSimRootOptions *opt = 0;

opttree->SetBranchAddress("wcsimrootoptions", &opt);

if(verbose) std::cout << "Options tree has: " << opttree->GetEntries() << " entries (1 expected)" << std::endl;

if (opttree->GetEntries() == 0) {

exit(9);

}

opttree->GetEntry(0);

opt->Print();

float win_scale = 0.75;

int n_wide(2);

int n_high(2);

TCanvas* c1 = new TCanvas("c1", "First canvas", 1000*n_wide*win_scale, 500*n_high*win_scale);

c1->Divide(2,1);

DumpGeoTree(geo, c1, verbose);

return 0;

}

WCSim/sample-root-scripts/plot_pmts.C at develop · WCSim/WCSim · GitHub

Uh oh!

Expand file tree

Breadcrumbs

plot_pmts.C

Latest commit

History

Breadcrumbs

plot_pmts.C

File metadata and controls