EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
qaplots.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file qaplots.cxx
1 // Load and smear a file,
2 // Then create a variety of QA plots
3 // This implementation mirrors the "standard" work flow
4 // Build -> Smear -> Load and befriend -> Analyze
5 
6 #include "qaplots.hh"
7 
8 #include <eicsmear/functions.h>
10 
15 #include "eicsmear/smear/EventS.h"
16 #include "eicsmear/smear/Smear.h"
19 
20 #include "eicsmeardetectors.hh"
21 
22 // Note: The remaining includes are not necessary for eic-smear usage
23 #include <TSystem.h>
24 #include <TDatabasePDG.h>
25 #include <TFile.h>
26 #include <TChain.h>
27 #include <TCanvas.h>
28 #include <TProfile.h>
29 #include <TLegend.h>
30 #include <TText.h>
31 #include <TStyle.h>
32 
33 #include <string>
34 #include <iomanip>
35 #include <cctype>
36 #include <exception>
37 #include <vector>
38 #include <map>
39 
40 // Convenience only
41 using std::cout;
42 using std::cerr;
43 using std::endl;
44 using std::string;
45 using std::vector;
46 using std::map;
47 
48 // some helpers
49 static const TString getrootname(const qaparameters& qapars );
50 static void initializepidqabook(const qaparameters& qapars, map<int,pidqacollection>& qabook );
51 static void initializeeventqa(const qaparameters& qapars, eventqacollection& eventqa );
52 static void FillEventQA( eventqacollection& eventqa , const erhic::EventMC* const inEvent, const Smear::Event* const inEventS );
53 static void FillParticleQA( map<int,pidqacollection>& qabook, const Particle* const inParticle, const Smear::ParticleMCS* const inParticleS );
54 static void PlotQA ( const qaparameters& qapars, eventqacollection& eventqa, map<int,pidqacollection>& qabook );
55 
56 int main(int argc, char* argv[]){
57 
58  // Parse arguments
59  // ---------------
60  // Defaults are set in qaplots.h
61  qaparameters qapars = ParseArguments ( argc, argv );
62 
63  // Set up output name
64  TString rootname = getrootname(qapars);
65 
66  // First try to instantiate the detector
67  // to avoid pointlessly transforming if that doesn't work
68  // ------------------------------------------------------
69  Smear::Detector detector;
70  if ( qapars.beam_mom_nn < 0 ) detector = BuildByName(qapars.detstring);
71  else detector = BuildByName(qapars.detstring, qapars.beam_mom_nn);
72 
73  if ( detector.GetNDevices() == 0 ) {
74  cerr << "Detector sepcified as " << qapars.detstring
75  << " not recognized or empty." << endl;
76  return -1;
77  }
78 
79 
80  // Convert input file to tree
81  // --------------------------
82  auto buildresult = BuildTree(qapars.txtfilename.c_str(), qapars.outpath.c_str(), qapars.nevents);
83  if ( buildresult !=0 ){
84  cerr << "Failed to build a tree from " << qapars.txtfilename << endl;
85  return buildresult;
86  }
87 
88  // Smear the tree
89  // --------------
90  TString smearedname = rootname;
91  smearedname.ReplaceAll (".root",".smeared.root" );
92  // Can disable warnings here.
93  // Many warnings are harmless ( x or y can be smeared to values >1)
94  // But it is recommended to leave it on and follow up on "inf", "nan" etc. if you test a new detector
95  // erhic::DisKinematics::BoundaryWarning=false;
96  SmearTree( detector, rootname.Data(), smearedname.Data());
97 
98  // -------------
99  // Load the tree
100  // -------------
101  TChain* inTree = new TChain("EICTree");
102  inTree->Add(rootname);
103  inTree->AddFriend("Smeared",smearedname);
104 
105  // Setup Input Event Buffer
106  erhic::EventMC* inEvent(NULL);
107  Smear::Event* inEventS(NULL);
108  inTree->SetBranchAddress("event",&inEvent);
109  inTree->SetBranchAddress("eventS",&inEventS);
110 
111  // Open histo root file and book histograms
112  // ----------------------------------------
113  TFile * outfile = new TFile ( qapars.outfilebase + qapars.detstring + ".root", "RECREATE");
114 
115  // A collection of event-wise qa plots, like kinematics
116  // eventqacollection is defined in the header
117  eventqacollection eventqa = {}; // this syntax initializes everything to 0
118  initializeeventqa ( qapars, eventqa );
119 
120  // We'll also have a collection of histos and maybe other info for every pid
121  // pidqacollection is defined in the header
122  map<int,pidqacollection> qabook;
123  // By default, use the standard particles
124  if ( qapars.pids.size() == 0 ) qapars.pids = { 11, 211, 321, 2212, 2112 }; // e, pi, K, p, n
125  initializepidqabook ( qapars, qabook );
126 
127  // -------------
128  // Analysis loop
129  // -------------
130  for(long iEvent=0; iEvent<inTree->GetEntries(); iEvent++){
131 
132  //Read next Event
133  if(inTree->GetEntry(iEvent) <=0) break;
134  if(iEvent%10000 == 0) cout << "Event " << iEvent << endl;
135 
136  // -------------
137  // event-wise QA
138  // -------------
139  // Following function just contains many statements of the form
140  // if ( inEventS->GetQ2()>0 ) eventqa.Q2_NM->Fill ( std::log10(inEvent->GetQ2()), std::log10(inEventS->GetQ2()));
141  FillEventQA( eventqa, inEvent, inEventS );
142 
143  // -------------------
144  // Loop over Particles
145  // -------------------
146  for(int j=0; j<inEventS->GetNTracks(); j++){
147  // Skip beam
148  if ( j<3 ) continue;
149 
150  const Smear::ParticleMCS* inParticleS = inEventS->GetTrack(j); // Smeared Particle
151  const Particle* inParticle = inEvent->GetTrack(j); // Unsmeared Particle
152 
153  // Skip non-final particles.
154  if ( inParticle->GetStatus() != 1 ) continue;
155 
156  // Particle was not smeared
157  if(inParticleS == NULL) continue;
158 
159  // ----------------
160  // Particle-wise QA
161  // ----------------
162  // Following function just contains many statements of the form
163  // coll.dEta_p->Fill(inParticle->GetP(), inParticle->GetEta() - inParticleS->GetEta());
164  FillParticleQA( qabook, inParticle, inParticleS );
165  }
166  }
167 
168  // NOTE: The remainder of this long file is tedious and explicit creation and filling of histograms
169  // The Fill*QA functions demonstrate how to access values in detail, but this is the end of the
170  // basic work flow to use eic-smear from start to finish!
171 
172  qapars.usedevents = inTree->GetEntries();
173  PlotQA( qapars, eventqa, qabook );
174  outfile->Write();
175 
176  return 0;
177 }
178 
179 // -----------------------------------------------------------------------------
180 // -----------------------------------------------------------------------------
181 
182 void FillEventQA( eventqacollection& eventqa , const erhic::EventMC* const inEvent, const Smear::Event* const inEventS ){
183 
184  // Q2
185  if ( eventqa.Q2_NM && inEventS->GetQ2()>0) eventqa.Q2_NM->Fill ( std::log10(inEvent->GetQ2()), std::log10(inEventS->GetQ2()));
186  else eventqa.missedQ2_NM++;
187  if ( eventqa.Q2_JB && inEventS->GetQ2JacquetBlondel()>0 ) eventqa.Q2_JB->Fill ( std::log10(inEvent->GetQ2()), std::log10(inEventS->GetQ2JacquetBlondel()));
188  else eventqa.missedQ2_JB++;
189  if ( eventqa.Q2_DA && inEventS->GetQ2DoubleAngle()>0 ) eventqa.Q2_DA->Fill ( std::log10(inEvent->GetQ2()), std::log10(inEventS->GetQ2DoubleAngle()));
190  else eventqa.missedQ2_DA++;
191 
192  if ( eventqa.delQ2_NM && inEventS->GetQ2()>0) eventqa.delQ2_NM->Fill ( std::log10(inEvent->GetQ2()), std::abs(inEventS->GetQ2()-inEvent->GetQ2())/ inEvent->GetQ2());
193  if ( eventqa.delQ2_JB && inEventS->GetQ2JacquetBlondel()>0 ) eventqa.delQ2_JB->Fill ( std::log10(inEvent->GetQ2()), std::abs(inEventS->GetQ2JacquetBlondel()-inEvent->GetQ2())/ inEvent->GetQ2());
194  if ( eventqa.delQ2_DA && inEventS->GetQ2DoubleAngle()>0 ) eventqa.delQ2_DA->Fill ( std::log10(inEvent->GetQ2()), (inEventS->GetQ2DoubleAngle()-inEvent->GetQ2())/ inEvent->GetQ2());
195 
196  // y
197  if ( eventqa.y_NM && inEventS->GetY()>0) eventqa.y_NM->Fill ( inEvent->GetY(), inEventS->GetY());
198  else eventqa.missedy_NM++;
199  if ( eventqa.y_JB && inEventS->GetYJacquetBlondel()>0 ) eventqa.y_JB->Fill ( inEvent->GetY(), inEventS->GetYJacquetBlondel());
200  else eventqa.missedy_JB++;
201  if ( eventqa.y_DA && inEventS->GetYDoubleAngle()>0 ) eventqa.y_DA->Fill ( inEvent->GetY(), inEventS->GetYDoubleAngle());
202  else eventqa.missedy_DA++;
203 
204  if ( eventqa.dely_NM && inEventS->GetY()>0) eventqa.dely_NM->Fill ( inEvent->GetY(), std::abs(inEventS->GetY()-inEvent->GetY())/ inEvent->GetY());
205  if ( eventqa.dely_JB && inEventS->GetYJacquetBlondel()>0 ) eventqa.dely_JB->Fill ( inEvent->GetY(), std::abs(inEventS->GetYJacquetBlondel()-inEvent->GetY())/ inEvent->GetY());
206  if ( eventqa.dely_DA && inEventS->GetYDoubleAngle()>0 ) eventqa.dely_DA->Fill ( inEvent->GetY(), std::abs(inEventS->GetYDoubleAngle()-inEvent->GetY())/ inEvent->GetY());
207 
208  // x
209  if ( eventqa.x_NM && inEventS->GetX()>0) eventqa.x_NM->Fill ( std::log10(inEvent->GetX()), std::log10(inEventS->GetX()));
210  else eventqa.missedx_NM++;
211  if ( eventqa.x_JB && inEventS->GetXJacquetBlondel()>0 ) eventqa.x_JB->Fill ( std::log10(inEvent->GetX()), std::log10(inEventS->GetXJacquetBlondel()));
212  else eventqa.missedx_JB++;
213  if ( eventqa.x_DA && inEventS->GetXDoubleAngle()>0 ) eventqa.x_DA->Fill ( std::log10(inEvent->GetX()), std::log10(inEventS->GetXDoubleAngle()));
214  else eventqa.missedx_DA++;
215 
216  if ( eventqa.delx_NM && inEventS->GetX()>0) eventqa.delx_NM->Fill ( std::log10(inEvent->GetX()), std::abs(inEventS->GetX()-inEvent->GetX())/ inEvent->GetX());
217  if ( eventqa.delx_JB && inEventS->GetXJacquetBlondel()>0 ) eventqa.delx_JB->Fill ( std::log10(inEvent->GetX()), std::abs(inEventS->GetXJacquetBlondel()-inEvent->GetX())/ inEvent->GetX());
218  if ( eventqa.delx_DA && inEventS->GetXDoubleAngle()>0 ) eventqa.delx_DA->Fill ( std::log10(inEvent->GetX()), std::abs(inEventS->GetXDoubleAngle()-inEvent->GetX())/ inEvent->GetX());
219 
220 }
221 // ---------------------------------------------------------------
222 void FillParticleQA( map<int,pidqacollection>& qabook, const Particle* const inParticle, const Smear::ParticleMCS* const inParticleS ){
223 
224  // If any component is smeared, all others are either smeared or 0 (meaning "not detected").
225  // Could detect the latter case (with a given accuracy):
226  // const double epsilon = 1e-9;
227 
228  // Fill histograms
229  for ( auto& pidcoll : qabook ){
230  auto& pid = pidcoll.first;
231  auto& coll = pidcoll.second;
232  if ( pid==0 || inParticle->GetPdgCode() == pid ){
233 
234  if ( inParticleS->IsPSmeared() ){
235  auto delP = (inParticle->GetP() - inParticleS->GetP()) / inParticle->GetP();
236  coll.DelP_th->Fill(inParticle->GetTheta(), delP);
237  coll.DelP_eta->Fill(inParticle->GetEta(), delP);
238  }
239 
240  if ( inParticleS->IsESmeared() ){
241  auto delE = (inParticle->GetE() - inParticleS->GetE()) / inParticle->GetE();
242  coll.DelE_E->Fill(inParticle->GetE(), delE);
243  coll.DelE_th->Fill(inParticle->GetTheta(), delE);
244  coll.DelE_eta->Fill(inParticle->GetEta(), delE);
245  }
246 
247  if ( inParticleS->IsThetaSmeared() ){
248  coll.dTh_p->Fill(inParticle->GetP(), inParticle->GetTheta() - inParticleS->GetTheta());
249  coll.dEta_p->Fill(inParticle->GetP(), inParticle->GetEta() - inParticleS->GetEta());
250  }
251 
252  if ( inParticleS->IsPhiSmeared() ){
253  coll.dPhi_p->Fill(inParticle->GetP(), inParticle->GetPhi() - inParticleS->GetPhi());
254  }
255  }
256  }
257 }
258 
259 // -----------------------------------------------------------------------------
260 
261 qaparameters ParseArguments ( int argc, char* argv[] ){
262  vector<string> arguments(argv + 1, argv + argc);
263  bool argsokay=true;
264  qaparameters qapars;
265  try{
266  for ( auto parg = arguments.begin() ; parg!=arguments.end() ; ++parg){
267  string arg=*parg;
268  if ( arg == "-h" ){
269  argsokay=false;
270  break;
271  } else if ( arg == "-o" ){
272  if (++parg == arguments.end() ){ argsokay=false; break; }
273  qapars.outfilebase=*parg;
274  } else if ( arg == "-i" ){
275  if (++parg ==arguments.end() ){ argsokay=false; break; }
276  qapars.txtfilename=*parg;
277  } else if ( arg == "-N" ){
278  if (++parg==arguments.end() ){ argsokay=false; break; }
279  qapars.nevents=std::stoi(parg->data());
280  } else if ( arg == "-addpid" ){
281  if ( ++parg == arguments.end() ){ argsokay=false; break; }
282  qapars.pids.push_back(std::stoi(parg->data()));
283  } else if ( arg == "-det" ){
284  if (++parg == arguments.end() ){ argsokay=false; break; }
285  qapars.detstring=*parg;
286  for (auto & c: qapars.detstring) c = toupper(c);
287  if ( TString(qapars.detstring).Contains("MATRIX") && TString(qapars.detstring).Contains("FF")){
288  if (++parg == arguments.end() ){ argsokay=false; break; }
289  qapars.beam_mom_nn = std::stoi(parg->data());
290  }
291  if ( TString(qapars.detstring).Contains("CORE") && !TString(qapars.detstring).Contains("B")){
292  if (++parg == arguments.end() ){ argsokay=false; break; }
293  qapars.beam_mom_nn = std::stod(parg->data());
294  }
295  } else {
296  argsokay=false;
297  break;
298  }
299  }
300  } catch ( const std::exception& e){
301  cerr << "Caught exception during argument parsing: "
302  << e.what() << endl;
303  argsokay=false;
304  }
305 
306  if ( !argsokay ) {
307  cerr << "usage: " << argv[0] << endl
308  << " [-i txtfilename] (Lund-style file)" << endl
309  << " [-o OutFileBase] (extension will be added)" << endl
310  << " [-N Nevents] (<0 for all)" << endl
311  << " [-addpid pid] (can be called multiple times)" << endl
312  << " [-det detstring] matrix02B3, matrix02B15, matrix01, matrixff [beam_mom_nn], handbook, perfect, beast, ephenix, zeus, jleic (capitalization does not matter.)" << endl
313  << endl;
314  throw std::runtime_error("Not a valid list of options");
315  }
316  for (auto & c: qapars.detstring) c = toupper(c);
317 
318  return qapars;
319 }
320 
321 // ---------------------------------------------------------------
322 const TString getrootname(const qaparameters& qapars ){
323  // The root file name is created by replacing the extension by ".root"
324 
325  TString rootname = qapars.txtfilename;
326 
327  // Remove zip extension, if there is one.
328  if ( rootname.EndsWith(".gz", TString::kIgnoreCase) ||
329  rootname.EndsWith(".zip", TString::kIgnoreCase) )
330  rootname.Replace(rootname.Last('.'), rootname.Length(), "");
331 
332  // Remove the remaining extension, if there is one.
333  if (rootname.Last('.') > -1) {
334  rootname.Replace(rootname.Last('.'), rootname.Length(), "");
335  } // if
336 
337  // BuildTree includes event number in partial transformation
338  if ( qapars.nevents>=0 ) {
339  rootname += ".";
340  rootname += qapars.nevents;
341  rootname += "event";
342  }
343  rootname += ".root";
344  rootname = gSystem->BaseName( rootname );
345  rootname.Prepend( qapars.outpath );
346  cout << " ======================= " << endl;
347  cout << " Transforming input file " << endl
348  << qapars.txtfilename << endl
349  << " into root file " << endl
350  << rootname << endl;
351  return rootname;
352 }
353 // ---------------------------------------------------------------
354 void initializepidqabook(const qaparameters& qapars, map<int,pidqacollection>& qabook ){
355  gStyle->SetHistLineColor(kRed); // for Profiles
356 
357  TString s;
358  float pmin = 0;
359  float pmax = 20;
360  int pbins = 80;
361 
362  float dpmin = -0.15;
363  float dpmax = 0.15;
364  int dpbins = 100;
365 
366  float emin = 0;
367  float emax = 20;
368  int ebins = 80;
369 
370  float demin = -1.5;
371  float demax = 1.5;
372  int debins = 100;
373 
374  float thmin = 0;
375  float thmax = TMath::Pi();
376  int thbins = 64;
377 
378  float dthmin = -0.1;
379  float dthmax = 0.1;
380  int dthbins = 100;
381 
382  float etamin = -5;
383  float etamax = 5;
384  int etabins = 100;
385 
386  float detamin = -0.1;
387  float detamax = 0.1;
388  int detabins = 100;
389 
390  float phimin = 0;
391  float phimax = TMath::TwoPi();
392  int phibins = 64;
393 
394  float dphimin = -0.1;
395  float dphimax = 0.1;
396  int dphibins = 100;
397 
398  for ( auto pid : qapars.pids ){
399  pid = abs ( pid ); // ignoring charge
400 
401  s = qapars.detstring + "_DelE_E_"; s += pid;
402  qabook[pid].DelE_E = new TH2D( s,s+";E;#Delta E/E", ebins, emin, emax, debins, demin, demax);
403 
404  s = qapars.detstring + "_dPhi_p_"; s += pid;
405  qabook[pid].dPhi_p = new TH2D( s,s+";p;#Delta#phi", pbins, pmin, pmax, dphibins, dphimin, dphimax );
406 
407  s = qapars.detstring + "_DelP_th_"; s += pid;
408  qabook[pid].DelP_th = new TH2D( s,s+";#theta;#Delta p/p", thbins, thmin, thmax, dpbins, dpmin, dpmax);
409 
410  s = qapars.detstring + "_DelE_th_"; s += pid;
411  qabook[pid].DelE_th = new TH2D( s,s+";#theta;#Delta E/E", thbins, thmin, thmax, debins, demin, demax);
412 
413  s = qapars.detstring + "_dTh_p_"; s += pid;
414  qabook[pid].dTh_p = new TH2D( s,s+";p;#Delta#theta", pbins, pmin, pmax, dthbins, dthmin, dthmax );
415 
416  s = qapars.detstring + "_DelP_eta_"; s += pid;
417  qabook[pid].DelP_eta = new TH2D( s,s+";#eta;#Delta p/p", etabins, etamin, etamax, dpbins, dpmin, dpmax);
418 
419  s = qapars.detstring + "_DelE_eta_"; s += pid;
420  qabook[pid].DelE_eta = new TH2D( s,s+";#eta;#Delta E/E", etabins, etamin, etamax, debins, demin, demax);
421 
422  s = qapars.detstring + "_dEta_p_"; s += pid;
423  qabook[pid].dEta_p = new TH2D( s,s+";p;#Delta#eta", pbins, pmin, pmax, detabins, detamin, detamax );
424 
425  }
426 
427 }
428 
429 // ---------------------------------------------------------------
430 void initializeeventqa(const qaparameters& qapars, eventqacollection& eventqa ){
431  // recording log of x, y, Q2
432 
433  TString s;
434 
435  float Q2min = 1e-3;
436  float Q2max = 1e5;
437  int logQ2bins = 240;
438 
439  s = qapars.detstring + "_LogQ2_NM";
440  eventqa.Q2_NM = new TH2D( s,s+";log Q^{2};log Q^{2}_{NM}", logQ2bins, std::log10(Q2min), std::log10(Q2max), logQ2bins, std::log10(Q2min), std::log10(Q2max));
441 
442  s = qapars.detstring + "_LogQ2_JB";
443  eventqa.Q2_JB = new TH2D( s,s+";log Q^{2};log Q^{2}_{JB}", logQ2bins, std::log10(Q2min), std::log10(Q2max), logQ2bins, std::log10(Q2min), std::log10(Q2max));
444 
445  s = qapars.detstring + "_LogQ2_DA";
446  eventqa.Q2_DA = new TH2D( s,s+";log Q^{2};log Q^{2}_{DA}", logQ2bins, std::log10(Q2min), std::log10(Q2max), logQ2bins, std::log10(Q2min), std::log10(Q2max));
447 
448  int delQ2bins = 100;
449  s = qapars.detstring + "_delQ2_NM";
450  eventqa.delQ2_NM = new TH2D( s,s+";log Q^{2};(Q^{2}_{NM}-Q^{2})/Q2", logQ2bins, std::log10(Q2min), std::log10(Q2max), delQ2bins, 0, 1.2);
451 
452  s = qapars.detstring + "_delQ2_JB";
453  eventqa.delQ2_JB = new TH2D( s,s+";log Q^{2};(Q^{2}_{JB}-Q^{2})/Q2", logQ2bins, std::log10(Q2min), std::log10(Q2max), delQ2bins, 0, 1.2);
454 
455  s = qapars.detstring + "_delQ2_DA";
456  eventqa.delQ2_DA = new TH2D( s,s+";log Q^{2};(Q^{2}_{DA}-Q^{2})/Q2", logQ2bins, std::log10(Q2min), std::log10(Q2max), delQ2bins, 0, 1.2);
457 
458 
459  float ymin = 0;
460  float ymax = 1.2; // see y>1 as well
461  int ybins = 100;
462 
463  s = qapars.detstring + "_y_NM";
464  eventqa.y_NM = new TH2D( s,s+";y;y_{NM}", ybins, ymin, ymax, ybins, ymin, ymax);
465 
466  s = qapars.detstring + "_y_JB";
467  eventqa.y_JB = new TH2D( s,s+";y;y_{JB}", ybins, ymin, ymax, ybins, ymin, ymax);
468 
469  s = qapars.detstring + "_y_DA";
470  eventqa.y_DA = new TH2D( s,s+";y;y_{DA}", ybins, ymin, ymax, ybins, ymin, ymax);
471 
472  int delybins = 100;
473  s = qapars.detstring + "_dely_NM";
474  eventqa.dely_NM = new TH2D( s,s+";y;(y_{NM}-y)/y", ybins, ymin, ymax, delybins, 0, 1.2);
475 
476  s = qapars.detstring + "_dely_JB";
477  eventqa.dely_JB = new TH2D( s,s+";y;(y_{JB}-y)/y", ybins, ymin, ymax, delybins, 0, 1.2);
478 
479  s = qapars.detstring + "_dely_DA";
480  eventqa.dely_DA = new TH2D( s,s+";y;(y_{DA}-y)/y", ybins, ymin, ymax, delybins, 0, 1.2);
481 
482  // float ymin = 1e-4;
483  // float ymax = 1e1; // see y>1 as well
484  // int logybins = 200;
485  // s = qapars.detstring + "_Logy_NM";
486  // eventqa.y_NM = new TH2D( s,s+";log y;log y_{NM}", logybins, std::log10(ymin), std::log10(ymax), logybins, std::log10(ymin), std::log10(ymax));
487 
488  // s = qapars.detstring + "_Logy_JB";
489  // eventqa.y_JB = new TH2D( s,s+";log y;log y_{JB}", logybins, std::log10(ymin), std::log10(ymax), logybins, std::log10(ymin), std::log10(ymax));
490 
491  // s = qapars.detstring + "_Logy_DA";
492  // eventqa.y_DA = new TH2D( s,s+";log y;log y_{DA}", logybins, std::log10(ymin), std::log10(ymax), logybins, std::log10(ymin), std::log10(ymax));
493 
494  // int delybins = 100;
495  // s = qapars.detstring + "_dely_NM";
496  // eventqa.dely_NM = new TH2D( s,s+";log y;(y_{NM}-y)/y", logybins, std::log10(ymin), std::log10(ymax), delybins, 0, 1.2);
497 
498  // s = qapars.detstring + "_dely_JB";
499  // eventqa.dely_JB = new TH2D( s,s+";log y;(y_{JB}-y)/y", logybins, std::log10(ymin), std::log10(ymax), delybins, 0, 1.2);
500 
501  // s = qapars.detstring + "_dely_DA";
502  // eventqa.dely_DA = new TH2D( s,s+";log y;(y_{DA}-y)/y", logybins, std::log10(ymin), std::log10(ymax), delybins, 0, 1.2);
503 
504  float xmin = 1e-4;
505  float xmax = 1e1; // see x>1 as well
506  int logxbins = 200;
507 
508  s = qapars.detstring + "_Logx_NM";
509  eventqa.x_NM = new TH2D( s,s+";log x;log x_{NM}", logxbins, std::log10(xmin), std::log10(xmax), logxbins, std::log10(xmin), std::log10(xmax));
510 
511  s = qapars.detstring + "_Logx_JB";
512  eventqa.x_JB = new TH2D( s,s+";log x;log x_{JB}", logxbins, std::log10(xmin), std::log10(xmax), logxbins, std::log10(xmin), std::log10(xmax));
513 
514  s = qapars.detstring + "_Logx_DA";
515  eventqa.x_DA = new TH2D( s,s+";log x;log x_{DA}", logxbins, std::log10(xmin), std::log10(xmax), logxbins, std::log10(xmin), std::log10(xmax));
516 
517  int delxbins = 100;
518  s = qapars.detstring + "_delx_NM";
519  eventqa.delx_NM = new TH2D( s,s+";log x;(x_{NM}-x)/x", logxbins, std::log10(xmin), std::log10(xmax), delxbins, 0, 1.2);
520 
521  s = qapars.detstring + "_delx_JB";
522  eventqa.delx_JB = new TH2D( s,s+";log x;(x_{JB}-x)/x", logxbins, std::log10(xmin), std::log10(xmax), delxbins, 0, 1.2);
523 
524  s = qapars.detstring + "_delx_DA";
525  eventqa.delx_DA = new TH2D( s,s+";log x;(x_{DA}-x)/x", logxbins, std::log10(xmin), std::log10(xmax), delxbins, 0, 1.2);
526 }
527 
528 // ---------------------------------------------------------------
529 void PlotQA ( const qaparameters& qapars, eventqacollection& eventqa, map<int,pidqacollection>& qabook ){
530 
531  // Stat position and size
532  // -----------------------
533  gStyle->SetStatX(0.25);
534  gStyle->SetStatW(0.15);
535  gStyle->SetStatY(0.9);
536  gStyle->SetStatH(0.15);
537 
538  // Position of the "Missed: " box
539  float missx = 0.55;
540  float missy = 0.2;
541  float missy2 = 0.8;
542  TText t;
543  t.SetNDC();
544 
545  // Suppress obnoxious flood of
546  // "Current canvas added to pdf file"
547  gErrorIgnoreLevel = kWarning;
548 
549  // prep a pdf collection
550  new TCanvas;
551  gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf[" );
552 
553  // event-wise qa
554 
555  // response-style
556  // NM
557  if ( eventqa.y_NM ) {
558  eventqa.y_NM->Draw("colz");
559  t.DrawText( missx,missy, Form("Missed: %ld / %ld",eventqa.missedy_NM, qapars.usedevents));
560  gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf" );
561  }
562  if ( eventqa.x_NM ) {
563  eventqa.x_NM->Draw("colz");
564  t.DrawText( missx,missy, Form("Missed: %ld / %ld",eventqa.missedx_NM, qapars.usedevents));
565  gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf" );
566  }
567  if ( eventqa.Q2_NM ) {
568  eventqa.Q2_NM->Draw("colz");
569  t.DrawText( missx,missy, Form("Missed: %ld / %ld",eventqa.missedQ2_NM, qapars.usedevents));
570  gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf" );
571  }
572  // DA
573  if ( eventqa.y_DA ) {
574  eventqa.y_DA->Draw("colz");
575  t.DrawText( missx,missy, Form("Missed: %ld / %ld",eventqa.missedy_DA, qapars.usedevents));
576  gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf" );
577  }
578  if ( eventqa.x_DA ) {
579  eventqa.x_DA->Draw("colz");
580  t.DrawText( missx,missy, Form("Missed: %ld / %ld",eventqa.missedx_DA, qapars.usedevents));
581  gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf" );
582  }
583  if ( eventqa.Q2_DA ) {
584  eventqa.Q2_DA->Draw("colz");
585  t.DrawText( missx,missy, Form("Missed: %ld / %ld",eventqa.missedQ2_DA, qapars.usedevents));
586  gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf" );
587  }
588  // JB
589  if ( eventqa.y_JB ) {
590  eventqa.y_JB->Draw("colz");
591  t.DrawText( missx,missy, Form("Missed: %ld / %ld",eventqa.missedy_JB, qapars.usedevents));
592  gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf" );
593  }
594  if ( eventqa.x_JB ) {
595  eventqa.x_JB->Draw("colz");
596  t.DrawText( missx,missy, Form("Missed: %ld / %ld",eventqa.missedx_JB, qapars.usedevents));
597  gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf" );
598  }
599  if ( eventqa.Q2_JB ) {
600  eventqa.Q2_JB->Draw("colz");
601  t.DrawText( missx,missy, Form("Missed: %ld / %ld",eventqa.missedQ2_JB, qapars.usedevents));
602  gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf" );
603  }
604 
605  // // resolution-style
606  // // NM
607  // if ( eventqa.dely_NM ) {
608  // eventqa.dely_NM->Draw("colz");
609  // t.DrawText( missx,missy2, Form("Missed: %ld / %ld",eventqa.missedy_NM, qapars.usedevents));
610  // gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf" );
611  // }
612  // if ( eventqa.delx_NM ) {
613  // eventqa.delx_NM->Draw("colz");
614  // t.DrawText( missx,missy2, Form("Missed: %ld / %ld",eventqa.missedx_NM, qapars.usedevents));
615  // gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf" );
616  // }
617  // if ( eventqa.delQ2_NM ) {
618  // eventqa.delQ2_NM->Draw("colz");
619  // t.DrawText( missx,missy2, Form("Missed: %ld / %ld",eventqa.missedQ2_NM, qapars.usedevents));
620  // gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf" );
621  // }
622  // // DA
623  // if ( eventqa.dely_DA ) {
624  // eventqa.dely_DA->Draw("colz");
625  // t.DrawText( missx,missy2, Form("Missed: %ld / %ld",eventqa.missedy_DA, qapars.usedevents));
626  // gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf" );
627  // }
628  // if ( eventqa.delx_DA ) {
629  // eventqa.delx_DA->Draw("colz");
630  // t.DrawText( missx,missy2, Form("Missed: %ld / %ld",eventqa.missedx_DA, qapars.usedevents));
631  // gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf" );
632  // }
633  // if ( eventqa.delQ2_DA ) {
634  // eventqa.delQ2_DA->Draw("colz");
635  // t.DrawText( missx,missy2, Form("Missed: %ld / %ld",eventqa.missedQ2_DA, qapars.usedevents));
636  // gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf" );
637  // }
638  // // JB
639  // if ( eventqa.dely_JB ) {
640  // eventqa.dely_JB->Draw("colz");
641  // t.DrawText( missx,missy2, Form("Missed: %ld / %ld",eventqa.missedy_JB, qapars.usedevents));
642  // gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf" );
643  // }
644  // if ( eventqa.delx_JB ) {
645  // eventqa.delx_JB->Draw("colz");
646  // t.DrawText( missx,missy2, Form("Missed: %ld / %ld",eventqa.missedx_JB, qapars.usedevents));
647  // gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf" );
648  // }
649  // if ( eventqa.delQ2_JB ) {
650  // eventqa.delQ2_JB->Draw("colz");
651  // t.DrawText( missx,missy2, Form("Missed: %ld / %ld",eventqa.missedQ2_JB, qapars.usedevents));
652  // gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf" );
653  // }
654 
655  // particle QA
656  // -----------
657  gStyle->SetStatX(0.55); // reposition stat box
658  for ( auto& pidcoll : qabook ){
659  auto& pid = pidcoll.first;
660  auto& coll = pidcoll.second;
661 
662  // option "s" in Profile shows rms
663 
664  coll.DelP_th->Draw("colz");
665  gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf" );
666 
667  coll.DelP_eta->Draw("colz");
668  coll.DelP_eta->ProfileX("_px",1,-1,"s")->Draw("same");
669  gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf" );
670 
671  gStyle->SetStatX(0.25); // reposition stat box
672  coll.DelE_th->Draw("colz");
673  gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf" );
674 
675  coll.DelE_eta->Draw("colz");
676  coll.DelE_eta->ProfileX("_px",1,-1,"s")->Draw("same");
677  gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf" );
678  gStyle->SetStatX(0.55); // reposition stat box
679 
680  coll.DelE_E->Draw("colz");
681  coll.DelE_E->ProfileX("_px",1,-1,"s")->Draw("same");
682  gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf" );
683 
684  coll.dTh_p->Draw("colz");
685  gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf" );
686 
687  coll.dEta_p->Draw("colz");
688  coll.dEta_p->ProfileX("_px",1,-1,"s")->Draw("same");
689  gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf" );
690 
691  coll.dPhi_p->Draw("colz");
692  coll.dPhi_p->ProfileX("_px",1,-1,"s")->Draw("same");
693  gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf" );
694  }
695 
696  // return to standard warning level
697  gErrorIgnoreLevel = kInfo;
698 
699  // close the pdf collection
700  gPad->SaveAs( qapars.outfilebase + qapars.detstring + ".pdf]" );
701 }
702 // ---------------------------------------------------------------