ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
prttools.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file prttools.C
1 // prttools - useful functions for hld*, prt*
2 // original author: Roman Dzhygadlo - GSI Darmstadt
3 
4 #ifndef prttools_h
5 #define prttools_h 1
6 
7 #include "TROOT.h"
8 #include "TSystem.h"
9 #include "TStyle.h"
10 #include "TCanvas.h"
11 #include "TPad.h"
12 #include "TH1.h"
13 #include "TH2.h"
14 #include "TGraph.h"
15 #include "TMultiGraph.h"
16 #include "TSpline.h"
17 #include "TF1.h"
18 #include "TFile.h"
19 #include "TTree.h"
20 #include "TClonesArray.h"
21 #include "TVector3.h"
22 #include "TMath.h"
23 #include "TChain.h"
24 #include "TGaxis.h"
25 #include "TColor.h"
26 #include "TString.h"
27 #include "TArrayD.h"
28 #include "TSpectrum.h"
29 #include "TSpectrum2.h"
30 #include "Math/Minimizer.h"
31 #include "Math/Factory.h"
32 #include "Math/Functor.h"
33 #include "TRandom2.h"
34 #include "TError.h"
35 #include "TPaveStats.h"
36 #include "TObjString.h"
37 #include "TApplication.h"
38 #include <TLegend.h>
39 #include <TAxis.h>
40 #include <TPaletteAxis.h>
41 #include <TRandom.h>
42 #include <TCutG.h>
43 #include <TKey.h>
44 #include "TPRegexp.h"
45 #include "TFitResult.h"
46 
47 
48 #include <iostream>
49 #include <fstream>
50 #include <sstream>
51 
52 #if defined(prt__sim) || defined(prt__beam) || defined(eic__sim)
53 
54 class PrtRun;
55 class PrtEvent;
56 class PrtHit;
57 PrtEvent* prt_event = 0;
58 #endif
59 
60 #if defined(prt__sim) || defined(prt__beam)
61 #include "datainfo.C"
62 DataInfo prt_data_info;
63 #endif
64 
65 #if defined(eic__sim)
66 const int prt_nmcp = 4 * 7;
67 const int prt_npix = 16 * 16;
68 #else
69 const int prt_nmcp = 24;// 12; // 12;
70 const int prt_npix = 16 * 16; // 64
71 #endif
72 
73 // const int prt_ntdc=16;
74 // TString prt_tdcsid[prt_ntdc] ={"10","11","12","13",
75 // "20","21","22","23",
76 // "780","781","782","783",
77 // "840","841","842","843"
78 // };
79 
80 //may2015
81 const int prt_ntdc_may2015=41;
82 TString prt_tdcsid_may2015[] ={"2000","2001","2002","2003","2004","2005","2006","2007","2008","2009",
83  "200a","200b","200c","200d","200e","200f","2010","2011","2012","2013",
84  "2014","2015","2016","2018","2019","201a","201c","2020","2023","2024",
85  "2025","2026","2027","2028","2029","202a","202b","202c","202d","202e","202f"
86 };
87 
88 //jun2015
89 const int prt_ntdc_jun2015=30;
90 TString prt_tdcsid_jun2015[] ={"2000","2001","2002","2003","2004","2005","2006","2007","2008","2009",
91  "200a","200b","200c","200d","200e","200f","2010","2011","2012","2013",
92  "2014","2015","2016","2018","2019","201a","201c","201d","202c","202d"
93 };
94 
95 //oct2016
96 const int prt_ntdc_oct2016=20;
97 TString prt_tdcsid_oct2016[] ={"2000","2001","2002","2003","2004","2005","2006","2007","2008","2009",
98  "200a","200b","200c",
99  "2018","201b","201c","201f","202c","202d","202d"
100 };
101 
102 //aug2017 jul2018
103 const int prt_ntdc_jul2018 = 32;
104 TString prt_tdcsid_jul2018 [] ={"2000","2001","2002","2003","2004","2005","2006","2007","2008","2009","200a","200b","200c","200d","200e","200f","2010","2011","2012","2013",
105  "2014","2015","2016","2017","2018","2019","201a","201b",
106  "201c","201d","201e","201f"
107 };
108 
109 //jul2019
110 const int prt_ntdc_jul2019 = 21;
111 TString prt_tdcsid_jul2019[] ={"2014","2015","2016","2017","2000","2001","2002","2003","2004","2005",
112  "2006","2007","2008","2009","200a","200b","200c","200d","200e","200f",
113  "2010"};
114 
115 const int prt_maxdircch = prt_nmcp*prt_npix;
116 const int prt_maxnametdc=10000;
117 
118 TRandom prt_rand;
119 TChain* prt_ch = 0;
120 int prt_entries(0),prt_particle(0),prt_geometry(2023),prt_beamx(0),prt_beamz(0);
121 int prt_last_layoutId,prt_last_max,prt_last_maxz,prt_last_minz;
122 double prt_theta(0), prt_test1(0),prt_test2(0),prt_mom(0),prt_phi(0);
123 TString prt_savepath(""),prt_info("");
124 TH2F* prt_hdigi[prt_nmcp];
125 TSpectrum *prt_spect = new TSpectrum(2);
126 
127 const int prt_ntdcm = 80; //41
128 int prt_ntdc = prt_ntdcm;
129 int prt_maxch = prt_ntdc*48;
130 TString prt_tdcsid[prt_ntdcm];
131 const int prt_maxchm = prt_ntdcm*48;
132 int map_tdc[prt_maxnametdc];
133 int map_mpc[prt_maxchm/prt_npix][prt_npix];
134 int map_mcp[prt_maxchm];
135 int map_pix[prt_maxchm];
136 int map_row[prt_maxchm];
137 int map_col[prt_maxchm];
138 
139 int prt_pid(0), prt_pdg[]={11,13,211,321,2212};
140 double prt_mass[]={0.000511,0.1056584,0.139570,0.49368,0.9382723};
141 TString prt_name[]={"e","muon","pion","kaon","proton"};
142 TString prt_lname[]={"e","#mu","#pi","K","p"};
143 int prt_color[]={kOrange+6,kCyan+1,kBlue+1,kRed+1,kRed+1};
144 double prt_particleArray[3000];
145 
146 TF1 *prt_gaust;
147 TVector3 prt_fit(TH1 *h, double range = 3, double threshold=20, double limit=2, int peakSearch=1,int bkg = 0, TString opt="MQ"){
148  int binmax = h->GetMaximumBin();
149  double xmax = h->GetXaxis()->GetBinCenter(binmax);
150  if(bkg==0) prt_gaust = new TF1("prt_gaust","[0]*exp(-0.5*((x-[1])/[2])^2)",xmax-range,xmax+range);
151  else if(bkg==1) prt_gaust = new TF1("prt_gaust","[0]*exp(-0.5*((x-[1])/[2])^2)+[3]+x*[4]",xmax-range,xmax+range);
152  prt_gaust->SetNpx(500);
153  prt_gaust->SetParNames("const","mean","sigma");
154  prt_gaust->SetLineColor(2);
155  double integral = h->Integral(h->GetXaxis()->FindBin(xmax-range),h->GetXaxis()->FindBin(xmax+range));
156  double xxmin, xxmax, sigma1(0), mean1(0), mean2(0);
157  xxmax = xmax;
158  xxmin = xxmax;
159  int nfound(1);
160  if(integral>threshold){
161 
162  if(peakSearch == 1){
163  prt_gaust->SetParameter(1,xmax);
164  prt_gaust->SetParameter(2,0.005);
165  prt_gaust->SetParLimits(2,0.003,limit);
166  h->Fit("prt_gaust",opt,"",xxmin-range, xxmax+range);
167  }
168 
169  if(peakSearch>1){
170  nfound = prt_spect->Search(h,4,"goff",0.1);
171  std::cout<<"nfound "<<nfound <<std::endl;
172  if(nfound==1){
173  prt_gaust =new TF1("prt_gaust","gaus(0)",xmax-range,xmax+range);
174  prt_gaust->SetNpx(500);
175  prt_gaust->SetParameter(1,prt_spect->GetPositionX()[0]);
176  }else if(nfound>=2){
177  double p1 = prt_spect->GetPositionX()[0];
178  double p2 = prt_spect->GetPositionX()[1];
179  if(p1>p2) {
180  xxmax = p1;
181  xxmin = p2;
182  }else {
183  xxmax = p2;
184  xxmin = p1;
185  }
186  if(peakSearch==20){
187  xxmax=xxmin;
188  prt_gaust =new TF1("prt_gaust","gaus(0)",xxmin-range,xxmin+range);
189  prt_gaust->SetNpx(500);
190  prt_gaust->SetParameter(1,prt_spect->GetPositionX()[0]);
191  }else{
192  prt_gaust =new TF1("prt_gaust","gaus(0)+gaus(3)",xmax-range,xmax+range);
193  prt_gaust->SetNpx(500);
194  prt_gaust->SetParameter(0,1000);
195  prt_gaust->SetParameter(3,1000);
196 
197  prt_gaust->FixParameter(1,xxmin);
198  prt_gaust->FixParameter(4,xxmax);
199  prt_gaust->SetParameter(2,0.1);
200  prt_gaust->SetParameter(5,0.1);
201  h->Fit("prt_gaust",opt,"",xxmin-range, xxmax+range);
202  prt_gaust->ReleaseParameter(1);
203  prt_gaust->ReleaseParameter(4);
204  }
205  }
206 
207  prt_gaust->SetParameter(2,0.2);
208  prt_gaust->SetParameter(5,0.2);
209  }
210 
211  h->Fit("prt_gaust",opt,"",xxmin-range, xxmax+range);
212  mean1 = prt_gaust->GetParameter(1);
213  sigma1 = prt_gaust->GetParameter(2);
214  if(sigma1>10) sigma1=10;
215 
216  if(peakSearch == 2){
217  mean2 = (nfound==1) ? prt_gaust->GetParameter(1) : prt_gaust->GetParameter(4);
218  // sigma2 = (nfound==1) ? prt_gaust->GetParameter(2) : prt_gaust->GetParameter(5);
219  }
220  }
221  delete prt_gaust;
222  return TVector3(mean1,sigma1,mean2);
223 }
224 
225 TGraph *prt_fitslices(TH2F *hh,double minrange=0, double maxrange=0, double fitrange=1,int rebin=1,int ret=0){
226  TH2F *h =(TH2F*) hh->Clone("h");
227  h->RebinY(rebin);
228  int point(0);
229  TGraph *gres = new TGraph();
230  for (int i=1;i<h->GetNbinsY();i++){
231  double x = h->GetYaxis()->GetBinCenter(i);
232  TH1D* hp;
233  if(minrange!=maxrange){
234  TCutG *cut = new TCutG("prt_onepeakcut",5);
235  cut->SetVarX("y");
236  cut->SetVarY("x");
237  cut->SetPoint(0,minrange,-1E6);
238  cut->SetPoint(1,minrange, 1E6);
239  cut->SetPoint(2,maxrange, 1E6);
240  cut->SetPoint(3,maxrange,-1E6);
241  cut->SetPoint(4,minrange,-1E6);
242 
243  hp = h->ProjectionX(Form("bin%d",i),i,i,"[prt_onepeakcut]");
244  }else{
245  hp = h->ProjectionX(Form("bin%d",i),i,i);
246  }
247 
248  TVector3 res = prt_fit((TH1F*)hp,fitrange,100,2,1,1);
249  double y=0;
250  if(ret==0) y = res.X();
251  if(ret==1) y = res.Y();
252  if(ret==2) y = res.X() + 0.5*res.Y();
253  if(ret==3) y = res.X() - 0.5*res.Y();
254  if(y==0 || y<minrange || y>maxrange) continue;
255 
256  gres->SetPoint(point,y,x);
257  gres->SetLineWidth(2);
258  gres->SetLineColor(kRed);
259  point++;
260  }
261  return gres;
262 }
263 
264 void prt_createMap(int setupid = 2019){
265  prt_geometry = setupid;
266  if(prt_geometry==2015) prt_ntdc = prt_ntdc_jun2015;
267  if(prt_geometry==2016) prt_ntdc = prt_ntdc_oct2016;
268  if(prt_geometry==2017) prt_ntdc = prt_ntdc_jul2018;
269  if(prt_geometry==2018) prt_ntdc = prt_ntdc_jul2018;
270  if(prt_geometry==2019) prt_ntdc = prt_ntdc_jul2019;
271 
272  for(int i=0; i<prt_ntdc; i++){
273  if(prt_geometry==2015) prt_tdcsid[i] = prt_tdcsid_jun2015[i];
274  if(prt_geometry==2016) prt_tdcsid[i] = prt_tdcsid_oct2016[i];
275  if(prt_geometry==2017) prt_tdcsid[i] = prt_tdcsid_jul2018[i];
276  if(prt_geometry==2018) prt_tdcsid[i] = prt_tdcsid_jul2018[i];
277  if(prt_geometry==2019) prt_tdcsid[i] = prt_tdcsid_jul2019[i];
278  }
279 
280  TGaxis::SetMaxDigits(4);
281  for(int i=0; i<prt_maxnametdc; i++) map_tdc[i]=-1;
282  for(int i=0; i<prt_ntdc; i++){
283  int dec = TString::BaseConvert(prt_tdcsid[i],16,10).Atoi();
284  map_tdc[dec]=i;
285  }
286 
287  // for(int ch=0; ch<prt_maxdircch; ch++){
288  for(int ch=0; ch<prt_maxch; ch++){
289  int mcp = ch/prt_npix;
290  int pix = ch%prt_npix;
291  int col = pix/2 - 8*(pix/16);
292  int row = pix%2 + 2*(pix/16);
293  // pix = col+sqrt(prt_npix)*row;
294 
295  map_mpc[mcp][pix]=ch;
296  map_mcp[ch] = mcp;
297  map_pix[ch] = pix;
298  map_row[ch] = row;
299  map_col[ch] = col;
300  }
301 
302  for(int i=0; i<5; i++){
303  prt_particleArray[prt_pdg[i]]=i;
304  }
305  prt_particleArray[212]=2;
306 }
307 
308 int prt_getChannelNumber(int tdc, int tdcChannel){
309  int ch = -1;
310  if(prt_geometry==2018){
311  ch=0;
312  for(int i=0; i<=tdc; i++){
313  if(i==tdc && (i==2 || i==6 || i==10 || i==14)) ch += tdcChannel-32;
314  else if(i==tdc) ch += tdcChannel;
315  else if(i==1 || i==2 || i==5 || i==6 || i==9 || i==10 || i==13 || i==14) ch += 16;
316  else ch += 48;
317  }
318  }else if(prt_geometry==2023){
319  ch = 32*tdc+tdcChannel;
320  }else{
321  ch = 48*tdc+tdcChannel;
322  }
323  return ch;
324 }
325 
326 int prt_getTdcId(int ch){
327  int tch=0, tdcid=0;
328  if(prt_geometry==2018){
329  for(int i=0; i<=prt_ntdc; i++){
330  tdcid=i;
331  if(i==2 || i==6 || i==10 || i==14) tch += 16;
332  else if(i==1 || i==2 || i==5 || i==6 || i==9 || i==10 || i==13 || i==14) tch += 16;
333  else tch += 48;
334  if(tch>ch) break;
335  }
336  }else if(prt_geometry==2023){
337  tdcid=ch/32;
338  }else{
339  tdcid=ch/48;
340  }
341  return tdcid;
342 }
343 
344 TString prt_getTdcName(int ch){
345  return prt_tdcsid[prt_getTdcId(ch)];
346 }
347 
348 int prt_getTdcChannel(int ch){
349  int tch=0,tdcc=0;
350  if(prt_geometry==2018){
351  for(int i=0; i<=prt_ntdc; i++){
352  tdcc=ch-tch+1;
353  if(i==2 || i==6 || i==10 || i==14 || i==1 || i==2 || i==5 || i==6 || i==9 || i==10 || i==13 || i==14) tch += 16;
354  else tch += 48;
355  if(tch>ch){
356  break;
357  };
358  }
359  }else if(prt_geometry==2023){
360  tdcc=ch%32+1;
361  }else{
362  tdcc=ch%48+1;
363  }
364  return tdcc;
365 }
366 
367 int prt_removeRefChannels(int ch, int tdcSeqId){
368  return ch - tdcSeqId;
369 }
370 
371 int prt_addRefChannels(int ch,int tdcSeqId){
372  return ch + tdcSeqId;
373 }
374 
375 Bool_t prt_isBadChannel(int ch){
376  if(ch<0 || ch>=prt_maxdircch) return true;
377 
378  // // bad pixels july15
379 
380  // if(ch==202) return true;
381  // if(ch==204) return true;
382  // if(ch==206) return true;
383  // if(ch==830) return true;
384  // if(ch==831) return true;
385  // if(ch==828) return true;
386  // if(ch==815) return true;
387  // if(ch>383 && ch<400) return true; //dead chain
388 
389  return false;
390 }
391 
392 TString prt_randstr(int len = 10){
393  TString str = "";
394  static const char alphanum[] =
395  "0123456789"
396  "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
397  "abcdefghijklmnopqrstuvwxyz";
398 
399  for (int i = 0; i < len; ++i) {
400  str += alphanum[rand() % (sizeof(alphanum) - 1)];
401  }
402  return str;
403 }
404 
405 // layoutId == 5 - 5 row's design for the PANDA Barrel DIRC
406 // layoutId == 2015 - cern 2015
407 // layoutId == 2016 - cern 2016
408 // layoutId == 2017 - cern 2017
409 // layoutId == 2018 - cern 2018
410 // layoutId == 2021 - new 3.6 row's design for the PANDA Barrel DIRC
411 // layoutId == 2023 - new 2x4 layout for the PANDA Barrel DIRC
412 // layoutId == 2031 - EIC DIRC beam test
413 // layoutId == 2030 - EIC DIRC prism
414 // layoutId == 2032 - EIC DIRC focusing prism
415 TH1 * prt_cdigi_th;
416 TCanvas *prt_drawDigi(int layoutId = 0, double maxz = 0, double minz = 0, TCanvas *cdigi = NULL){
417 
418  prt_last_layoutId=layoutId;
419  prt_last_maxz=maxz;
420  prt_last_minz=minz;
421 
422  if(prt_geometry==2021) layoutId=2021;
423  // if(prt_geometry==2019) layoutId=2018;
424 
425  TString sid = prt_randstr(3);
426  if(cdigi) cdigi->cd();
427  else cdigi = new TCanvas("hp="+sid,"hp_"+sid,800,400);
428 
429  TPad* prt_hpads[prt_nmcp];
430  TPad* prt_hpglobal;
431 
432  if(layoutId==2015 || layoutId==5) prt_hpglobal = new TPad(sid,"T",0.04,0.04,0.88,0.96);
433  else if(layoutId==2021) prt_hpglobal = new TPad(sid,"T",0.12,0.02,0.78,0.98);
434  else if(layoutId==2016) prt_hpglobal = new TPad(sid,"T",0.2,0.02,0.75,0.98);
435  else if(layoutId==2017) prt_hpglobal = new TPad(sid,"T",0.15,0.02,0.80,0.98);
436  else if(layoutId==2018) prt_hpglobal = new TPad(sid,"T",0.05,0.07,0.9,0.93);
437  else if(layoutId==2023) prt_hpglobal = new TPad(sid,"T",0.073,0.02,0.877,0.98);
438  else if(layoutId==2030) prt_hpglobal = new TPad(sid,"T",0.10,0.025,0.82,0.975);
439  else if(layoutId==2032) prt_hpglobal = new TPad(sid,"T",0.04,0.025,0.91,0.975);
440  else if(layoutId==2031) prt_hpglobal = new TPad(sid,"T",0.12,0.01,0.80,0.99);
441  else prt_hpglobal = new TPad(sid,"T",0.04,0.04,0.96,0.96);
442 
443  prt_hpglobal->SetFillStyle(0);
444  prt_hpglobal->Draw();
445  prt_hpglobal->cd();
446 
447  int nrow = 3, ncol = 5;
448  if(layoutId ==2016) ncol=3;
449  if(layoutId ==2017) ncol=4;
450  if(layoutId ==2018 || layoutId ==2023) {nrow=2; ncol=4;}
451  if(layoutId ==2021) ncol=4;
452  if(layoutId ==2030) {nrow=4; ncol=6;}
453  if(layoutId ==2032) {nrow=4; ncol=7;}
454  if(layoutId ==2031) {nrow=3; ncol=4;}
455 
456  if(layoutId > 1){
457  float tbw(0.02), tbh(0.01), shift(0),shiftw(0.02),shifth(0),margin(0.01);
458  int padi(0);
459 
460  for(int i=0; i<ncol; i++){
461  for(int j=0; j<nrow; j++){
462  if(j==1) shift = -0.028;
463  else shift = 0;
464  shifth=0;
465  if(layoutId == 5) {shift =0; shiftw=0.001; tbw=0.001; tbh=0.001;}
466  if(layoutId == 2021) {
467  if(i==0 && j == nrow-1) continue;
468  shift =0; shiftw=0.001; tbw=0.001; tbh=0.001;
469  if(i==0) shifth=0.167;
470  }
471  if(layoutId == 2016) {
472  shift = -0.01; shiftw=0.01; tbw=0.03; tbh=0.006;
473  if(j==1) shift += 0.015;
474  }
475  if(layoutId == 2017) {
476  margin= 0.1;
477  shift = 0; shiftw=0.01; tbw=0.005; tbh=0.006;
478  //if(j==1) shift += 0.015;
479  }
480  if(layoutId == 2018) {
481  margin= 0.1;
482  shift = 0; shiftw=0.01; tbw=0.005; tbh=0.006;
483  }
484  if(layoutId == 2023) {
485  margin= 0.1;
486  shift = 0; shiftw=0.01; tbw=0.0015; tbh=0.042;
487  }
488  if(layoutId == 2030) {
489  margin= 0.1;
490  shift = 0; shiftw=0.01; tbw=0.001; tbh=0.001;
491  padi=j*ncol+i;
492  }
493  if(layoutId == 2032) {
494  margin= 0.1;
495  shift = 0; shiftw=0.01; tbw=0.001; tbh=0.001;
496  padi=j*ncol+i;
497  }
498  if(layoutId == 2031) {
499  margin= 0.1;
500  shift = 0; shiftw=0.01; tbw=0.001; tbh=0.001;
501  padi=i*nrow+j;
502  }
503 
504  prt_hpads[padi] = new TPad(sid+Form("P%d",i*10+j),"T",
505  i/(ncol+2*margin)+tbw+shift+shiftw,
506  j/(double)nrow+tbh+shifth,
507  (i+1)/(ncol+2*margin)-tbw+shift+shiftw,
508  (1+j)/(double)nrow-tbh+shifth, 21);
509  prt_hpads[padi]->SetFillColor(kCyan-8);
510  prt_hpads[padi]->SetMargin(0.055,0.055,0.055,0.055);
511  prt_hpads[padi]->Draw();
512  padi++;
513  }
514  }
515  }else{
516  float tbw(0.02), tbh(0.01), shift(0),shiftw(-0.02);
517  int padi(0);
518  for(int ii=0; ii<ncol; ii++){
519  for(int j=0; j<nrow; j++){
520  if(j==1) shift = 0.04;
521  else shift = 0;
522  prt_hpads[padi] = new TPad(Form("P%d",ii*10+j),"T", ii/(double)ncol+tbw+shift+shiftw , j/(double)nrow+tbh, (ii+1)/(double)ncol-tbw+shift+shiftw, (1+j)/(double)nrow-tbh, 21);
523  prt_hpads[padi]->SetFillColor(kCyan-8);
524  prt_hpads[padi]->SetMargin(0.04,0.04,0.04,0.04);
525  prt_hpads[padi]->Draw();
526  padi++;
527  }
528  }
529  }
530 
531  int np;
532  double max=0;
533 
534  {
535  double tmax;
536  if(maxz==0){
537  for(int p=0; p<nrow*ncol;p++){
538  tmax = prt_hdigi[p]->GetBinContent(prt_hdigi[p]->GetMaximumBin());
539  if(max<tmax) max = tmax;
540  }
541  }else{
542  max = maxz;
543  }
544 
545  if(maxz==-2 || minz==-2){ // optimize range
546  for(int p=0; p<nrow*ncol;p++){
547  tmax = prt_hdigi[p]->GetMaximum();
548  if(max<tmax) max = tmax;
549  }
550  int tbins = 2000;
551  TH1F *h = new TH1F("","",tbins,0,max);
552  for(int p=0; p<nrow*ncol;p++){
553  for(int i=0; i<8; i++){
554  for(int j=0; j<8; j++){
555  double val = prt_hdigi[p]->GetBinContent(i+1,j+1);
556  if(val!=0) h->Fill(val);
557 
558  }
559  }
560  }
561  double integral;
562  for(int i=0; i<tbins; i++){
563  integral = h->Integral(0,i);
564  if(integral>0) {
565  if(minz==-2) minz = h->GetBinCenter(i);
566  break;
567  }
568  }
569 
570  for(int i=tbins; i>0; i--){
571  integral = h->Integral(i,tbins);
572  if(integral>10) {
573  if(maxz==-2) max = h->GetBinCenter(i);
574  break;
575  }
576  }
577  }
578  }
579 
580  prt_last_max = max;
581  int nnmax(0);
582  for(int p=0; p<nrow*ncol;p++){
583  if(layoutId == 1 || layoutId == 4) np =p%nrow*ncol + p/3;
584  else np = p;
585 
586  if(layoutId == 6 && p>10) continue;
587 
588  prt_hpads[p]->cd();
589  prt_hdigi[np]->Draw("col");//"col+text"
590  if(maxz==-1) max = prt_hdigi[np]->GetBinContent(prt_hdigi[np]->GetMaximumBin());
591  if(nnmax<prt_hdigi[np]->GetEntries()) nnmax=np;
592  prt_hdigi[np]->SetMaximum(max);
593  prt_hdigi[np]->SetMinimum(minz);
594  }
595 
596  // prt_cdigi_palette = (TPaletteAxis*)prt_hdigi[nnmax]->GetListOfFunctions()->FindObject("palette");
597  // prt_cdigi_palette->SetX1NDC(0.89);
598  // prt_cdigi_palette->SetY1NDC(0.1);
599  // prt_cdigi_palette->SetX2NDC(0.93);
600  // prt_cdigi_palette->SetY2NDC(0.90);
601 
602  cdigi->cd();
603  TPaletteAxis* prt_cdigi_palette;
604  if(layoutId==2018 || layoutId==2023) prt_cdigi_palette = new TPaletteAxis(0.89,0.1,0.93,0.90,(TH1 *) prt_hdigi[nnmax]);
605  else if(layoutId==2032) prt_cdigi_palette = new TPaletteAxis(0.91,0.1,0.94,0.90,(TH1 *) prt_hdigi[nnmax]);
606  else prt_cdigi_palette = new TPaletteAxis(0.82,0.1,0.86,0.90,(TH1 *) prt_hdigi[nnmax]);
607  prt_cdigi_palette->SetName("prt_palette");
608  prt_cdigi_palette->Draw();
609 
610  cdigi->Modified();
611  cdigi->Update();
612 
613  return cdigi;
614 }
615 
616 TString prt_getPixData(TString s="m,p,v\n", int layoutId=1){
617  int nrow = 3, ncol = 5, np, nmax=0, npix=8;
618  if(layoutId ==2016) ncol=3;
619  if(layoutId ==2017) ncol=4;
620  if(layoutId ==2018 || layoutId ==2023) {nrow=2; ncol=4;}
621  if(layoutId ==2021) ncol=4;
622  if(layoutId ==2030) {nrow=4; ncol=6; npix=16;}
623  if(layoutId ==2032) {nrow=4; ncol=7; npix=16;}
624  if(layoutId ==2031) {nrow=3; ncol=4; npix=16;}
625 
626  int nnmax(0);
627  for(int p=0; p<nrow*ncol;p++){
628  if(layoutId == 1 || layoutId == 4) np =p%nrow*ncol + p/3;
629  else if(layoutId == 2030) np =p%ncol*nrow + p/ncol;
630  else np = p;
631 
632  if(prt_last_maxz==-1) prt_last_max = prt_hdigi[p]->GetBinContent(prt_hdigi[p]->GetMaximumBin());
633  if(nnmax<prt_hdigi[p]->GetEntries()) nnmax=p;
634  prt_hdigi[p]->SetMaximum(prt_last_max);
635  prt_hdigi[p]->SetMinimum(prt_last_minz);
636  for(int i=1; i<=npix; i++){
637  for(int j=1; j<=npix; j++){
638  double weight = (double)(prt_hdigi[p]->GetBinContent(i,j))/(double)prt_last_max*255;
639  if(weight>255) weight=255;
640  if(weight > 0) s += Form("%d,%d,%d\n", np, (i-1)*npix+j-1, (int)weight);
641  }
642  }
643  }
644  return s;
645 }
646 
647 void prt_initDigi(int type=1){
648  if(type == 1){
649  for(int m=0; m<prt_nmcp;m++){
650  if(prt_hdigi[m]) prt_hdigi[m]->Reset("M");
651  else{
652  prt_hdigi[m] = new TH2F( Form("mcp%d", m),Form("mcp%d", m),8,0.,8.,8,0.,8.);
653  prt_hdigi[m]->SetStats(0);
654  prt_hdigi[m]->SetTitle(0);
655  prt_hdigi[m]->GetXaxis()->SetNdivisions(10);
656  prt_hdigi[m]->GetYaxis()->SetNdivisions(10);
657  prt_hdigi[m]->GetXaxis()->SetLabelOffset(100);
658  prt_hdigi[m]->GetYaxis()->SetLabelOffset(100);
659  prt_hdigi[m]->GetXaxis()->SetTickLength(1);
660  prt_hdigi[m]->GetYaxis()->SetTickLength(1);
661  prt_hdigi[m]->GetXaxis()->SetAxisColor(15);
662  prt_hdigi[m]->GetYaxis()->SetAxisColor(15);
663  }
664  }
665  }
666  if(type == 2){ //eic
667  for(int m=0; m<prt_nmcp;m++){
668  if(prt_hdigi[m]) prt_hdigi[m]->Reset("M");
669  else{
670  prt_hdigi[m] = new TH2F( Form("mcp%d", m),Form("mcp%d", m),16,0,16,16,0,16);
671  prt_hdigi[m]->SetStats(0);
672  prt_hdigi[m]->SetTitle(0);
673  prt_hdigi[m]->GetXaxis()->SetNdivisions(20);
674  prt_hdigi[m]->GetYaxis()->SetNdivisions(20);
675  prt_hdigi[m]->GetXaxis()->SetLabelOffset(100);
676  prt_hdigi[m]->GetYaxis()->SetLabelOffset(100);
677  prt_hdigi[m]->GetXaxis()->SetTickLength(1);
678  prt_hdigi[m]->GetYaxis()->SetTickLength(1);
679  prt_hdigi[m]->GetXaxis()->SetAxisColor(15);
680  prt_hdigi[m]->GetYaxis()->SetAxisColor(15);
681  }
682  }
683  }
684 }
685 
686 void prt_resetDigi(){
687  for(int m=0; m<prt_nmcp;m++){
688  prt_hdigi[m]->Reset("M");
689  }
690 }
691 
692 void prt_axisHits800x500(TH2 * hist){
693  hist->SetStats(0);
694  hist->SetTitle(Form("%d hits",(int)hist->GetEntries()));
695  hist->GetXaxis()->SetTitle("z, [cm]");
696  hist->GetXaxis()->SetTitleSize(0.05);
697  hist->GetXaxis()->SetTitleOffset(0.8);
698  hist->GetYaxis()->SetTitle("y, [cm]");
699  hist->GetYaxis()->SetTitleSize(0.05);
700  hist->GetYaxis()->SetTitleOffset(0.7);
701 }
702 
703 void prt_axisAngle800x500(TH2 * hist){
704  hist->SetStats(0);
705  hist->SetTitle(Form("%d hits",(int)hist->GetEntries()));
706  hist->GetXaxis()->SetTitle("#theta, [degree]");
707  hist->GetXaxis()->SetTitleSize(0.05);
708  hist->GetXaxis()->SetTitleOffset(0.8);
709  hist->GetYaxis()->SetTitle("photons per track, [#]");
710  hist->GetYaxis()->SetTitleSize(0.05);
711  hist->GetYaxis()->SetTitleOffset(0.7);
712 }
713 
714 void prt_axisAngle800x500(TH1 * hist){
715  hist->SetStats(0);
716  hist->SetTitle(Form("%d hits",(int)hist->GetEntries()));
717  hist->GetXaxis()->SetTitle("#theta, [degree]");
718  hist->GetXaxis()->SetTitleSize(0.05);
719  hist->GetXaxis()->SetTitleOffset(0.8);
720  hist->GetYaxis()->SetTitle("photons per track, [#]");
721  hist->GetYaxis()->SetTitleSize(0.05);
722  hist->GetYaxis()->SetTitleOffset(0.7);
723 }
724 
725 void prt_axisTime800x500(TH2 * hist){
726  hist->GetXaxis()->SetTitle("time, [ns]");
727  hist->GetXaxis()->SetTitleSize(0.05);
728  hist->GetXaxis()->SetTitleOffset(0.8);
729  hist->GetYaxis()->SetTitle("entries, #");
730  hist->GetYaxis()->SetTitleSize(0.05);
731  hist->GetYaxis()->SetTitleOffset(0.7);
732  hist->SetLineColor(1);
733 }
734 
735 void prt_axisTime800x500(TH1 * hist, TString xtitle = "time [ns]"){
736  TGaxis::SetMaxDigits(3);
737  hist->GetXaxis()->SetTitle(xtitle);
738  hist->GetXaxis()->SetTitleSize(0.06);
739  hist->GetXaxis()->SetTitleOffset(0.8);
740  hist->GetXaxis()->SetLabelSize(0.05);
741  hist->GetYaxis()->SetTitle("entries [#]");
742  hist->GetYaxis()->SetTitleSize(0.06);
743  hist->GetYaxis()->SetTitleOffset(0.7);
744  hist->GetYaxis()->SetLabelSize(0.05);
745  hist->SetLineColor(1);
746 }
747 
748 void prt_setPrettyStyle(){
749  // Canvas printing details: white bg, no borders.
750  gStyle->SetCanvasColor(0);
751  gStyle->SetCanvasBorderMode(0);
752  gStyle->SetCanvasBorderSize(0);
753 
754  // Canvas frame printing details: white bg, no borders.
755  gStyle->SetFrameFillColor(0);
756  gStyle->SetFrameBorderMode(0);
757  gStyle->SetFrameBorderSize(0);
758 
759  // Plot title details: centered, no bg, no border, nice font.
760  gStyle->SetTitleX(0.1);
761  gStyle->SetTitleW(0.8);
762  gStyle->SetTitleBorderSize(0);
763  gStyle->SetTitleFillColor(0);
764 
765  // Font details for titles and labels.
766  gStyle->SetTitleFont(42, "xyz");
767  gStyle->SetTitleFont(42, "pad");
768  gStyle->SetLabelFont(42, "xyz");
769  gStyle->SetLabelFont(42, "pad");
770 
771  // Details for stat box.
772  gStyle->SetStatColor(0);
773  gStyle->SetStatFont(42);
774  gStyle->SetStatBorderSize(1);
775  gStyle->SetStatX(0.975);
776  gStyle->SetStatY(0.9);
777 
778  // gStyle->SetOptStat(0);
779 }
780 
781 void prt_setGStyle(TGraph *g, int id){
782  int prt_coll[]={kBlack,kRed+1,kGreen, kBlue, 4,kCyan-6,kOrange, 7,8,9,10,1,1,1};
783  int prt_colm[]={kBlack,kRed+1,kGreen+2,kBlue+1,4,kCyan-6,kOrange+1,7,8,9,10,1,1,1};
784 
785  int cl=(id<10)? prt_coll[id]:id;
786  int cm=(id<10)? prt_colm[id]:id;
787  g->SetLineColor(cl);
788  g->SetMarkerColor(cm);
789  g->SetMarkerStyle(20);
790  g->SetMarkerSize(0.8);
791  g->SetName(Form("gr_%d",id));
792 }
793 
794 void prt_setRootPalette(int pal = 0){
795 
796  // pal = 1: rainbow\n"
797  // pal = 2: reverse-rainbow\n"
798  // pal = 3: amber\n"
799  // pal = 4: reverse-amber\n"
800  // pal = 5: blue/white\n"
801  // pal = 6: white/blue\n"
802  // pal = 7: red temperature\n"
803  // pal = 8: reverse-red temperature\n"
804  // pal = 9: green/white\n"
805  // pal = 10: white/green\n"
806  // pal = 11: orange/blue\n"
807  // pal = 12: blue/orange\n"
808  // pal = 13: white/black\n"
809  // pal = 14: black/white\n"
810 
811  const int NRGBs = 5;
812  const int NCont = 255;
813  gStyle->SetNumberContours(NCont);
814 
815  if (pal < 1 && pal> 15) return;
816  else pal--;
817 
818  double stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
819  double red[15][NRGBs] = {{ 0.00, 0.00, 0.87, 1.00, 0.51 },
820  { 0.51, 1.00, 0.87, 0.00, 0.00 },
821  { 0.17, 0.39, 0.62, 0.79, 1.00 },
822  { 1.00, 0.79, 0.62, 0.39, 0.17 },
823  { 0.00, 0.00, 0.00, 0.38, 1.00 },
824  { 1.00, 0.38, 0.00, 0.00, 0.00 },
825  { 0.00, 0.50, 0.89, 0.95, 1.00 },
826  { 1.00, 0.95, 0.89, 0.50, 0.00 },
827  { 0.00, 0.00, 0.38, 0.75, 1.00 },
828  { 0.00, 0.34, 0.61, 0.84, 1.00 },
829  { 0.75, 1.00, 0.24, 0.00, 0.00 },
830  { 0.00, 0.00, 0.24, 1.00, 0.75 },
831  { 0.00, 0.34, 0.61, 0.84, 1.00 },
832  { 1.00, 0.84, 0.61, 0.34, 0.00 },
833  { 0.00, 0.00, 0.80, 1.00, 0.80 }
834  };
835  double green[15][NRGBs] = {{ 0.00, 0.81, 1.00, 0.20, 0.00 },
836  { 0.00, 0.20, 1.00, 0.81, 0.00 },
837  { 0.01, 0.02, 0.39, 0.68, 1.00 },
838  { 1.00, 0.68, 0.39, 0.02, 0.01 },
839  { 0.00, 0.00, 0.38, 0.76, 1.00 },
840  { 1.00, 0.76, 0.38, 0.00, 0.00 },
841  { 0.00, 0.00, 0.27, 0.71, 1.00 },
842  { 1.00, 0.71, 0.27, 0.00, 0.00 },
843  { 0.00, 0.35, 0.62, 0.85, 1.00 },
844  { 1.00, 0.75, 0.38, 0.00, 0.00 },
845  { 0.24, 1.00, 0.75, 0.18, 0.00 },
846  { 0.00, 0.18, 0.75, 1.00, 0.24 },
847  { 0.00, 0.34, 0.61, 0.84, 1.00 },
848  { 1.00, 0.84, 0.61, 0.34, 0.00 },
849  { 0.00, 0.85, 1.00, 0.30, 0.00 }
850  };
851  double blue[15][NRGBs] = {{ 0.51, 1.00, 0.12, 0.00, 0.00 },
852  { 0.00, 0.00, 0.12, 1.00, 0.51 },
853  { 0.00, 0.09, 0.18, 0.09, 0.00 },
854  { 0.00, 0.09, 0.18, 0.09, 0.00 },
855  { 0.00, 0.47, 0.83, 1.00, 1.00 },
856  { 1.00, 1.00, 0.83, 0.47, 0.00 },
857  { 0.00, 0.00, 0.00, 0.40, 1.00 },
858  { 1.00, 0.40, 0.00, 0.00, 0.00 },
859  { 0.00, 0.00, 0.00, 0.47, 1.00 },
860  { 1.00, 0.47, 0.00, 0.00, 0.00 },
861  { 0.00, 0.62, 1.00, 0.68, 0.12 },
862  { 0.12, 0.68, 1.00, 0.62, 0.00 },
863  { 0.00, 0.34, 0.61, 0.84, 1.00 },
864  { 1.00, 0.84, 0.61, 0.34, 0.00 },
865  { 0.60, 1.00, 0.10, 0.00, 0.00 }
866  };
867 
868 
869  TColor::CreateGradientColorTable(NRGBs, stops, red[pal], green[pal], blue[pal], NCont);
870 
871 }
872 
873 #if defined(prt__sim) || defined(eic__sim)
874 bool prt_init(TString inFile="../build/hits.root", int bdigi=0, TString savepath="", int setupid = 2019){
875 
876  if(inFile=="" || gSystem->AccessPathName(inFile)) return false;
877  if(savepath!="") prt_savepath=savepath;
878  prt_setRootPalette(1);
879  prt_createMap(setupid);
880  delete prt_ch;
881 
882  prt_ch = new TChain("data");
883 
884  prt_ch->Add(inFile);
885  prt_ch->SetBranchAddress("PrtEvent", &prt_event);
886 
887  prt_entries = prt_ch->GetEntries();
888  std::cout<<"Entries in chain: "<<prt_entries <<std::endl;
889  if(bdigi) prt_initDigi(bdigi);
890  return true;
891 }
892 
893 void prt_nextEvent(int ievent, int printstep){
894  prt_ch->GetEntry(ievent);
895  if(ievent%printstep==0 && ievent!=0) std::cout<<"Event # "<<ievent<< " # hits "<<prt_event->GetHitSize()<<std::endl;
896  if(ievent == 0){
897  if(gROOT->GetApplication()){
898  TIter next(gROOT->GetApplication()->InputFiles());
899  TObjString *os=0;
900  while((os = (TObjString*)next())){
901  prt_info += os->GetString()+" ";
902  }
903  prt_info += "\n";
904  }
905  prt_info += prt_event->PrintInfo();
906  prt_mom = prt_event->GetMomentum().Mag() +0.01;
907  prt_theta = prt_event->GetAngle();
908  prt_phi = prt_event->GetPhi();
909  prt_geometry= prt_event->GetGeometry();
910  prt_beamx= prt_event->GetBeamX();
911  prt_beamz= prt_event->GetBeamZ();
912  prt_test1 = prt_event->GetTest1();
913  prt_test2 = prt_event->GetTest2();
914  }
915  prt_particle = prt_event->GetParticle();
916  if(prt_event->GetParticle()<3000 && prt_event->GetParticle()>0){
917  prt_pid=prt_particleArray[prt_event->GetParticle()];
918  }
919 }
920 #endif
921 
922 #ifdef prt__beam
923 bool prt_init(TString inFile="../build/hits.root", int bdigi=0, TString savepath="", int setupid=2019){
924 
925  if(inFile=="") return false;
926  if(savepath!="") prt_savepath=savepath;
927 
928  prt_createMap(setupid);
929  prt_setRootPalette(1);
930  delete prt_ch;
931 
932  prt_ch = new TChain("data");
933 
934  prt_ch->Add(inFile);
935  prt_ch->SetBranchAddress("PrtEvent", &prt_event);
936 
937  // prt_ch->SetBranchStatus("fHitArray.fLocalPos", 0);
938  // prt_ch->SetBranchStatus("fHitArray.fGlobalPos", 0);
939  // prt_ch->SetBranchStatus("fHitArray.fDigiPos", 0);
940  // prt_ch->SetBranchStatus("fHitArray.fMomentum", 0);
941  // prt_ch->SetBranchStatus("fHitArray.fPosition", 0);
942 
943  prt_ch->SetBranchStatus("fHitArray.fParentParticleId", 0);
944  prt_ch->SetBranchStatus("fHitArray.fNreflectionsInPrizm", 0);
945  prt_ch->SetBranchStatus("fHitArray.fPathInPrizm", 0);
946  prt_ch->SetBranchStatus("fHitArray.fCherenkovMC", 0);
947 
948  prt_ch->SetBranchStatus("fPosition", 0);
949 
950  prt_entries = prt_ch->GetEntries();
951  std::cout<<"Entries in chain: "<<prt_entries <<std::endl;
952  if(bdigi) prt_initDigi(bdigi);
953  return true;
954 }
955 
956 void prt_nextEvent(int ievent, int printstep){
957  prt_ch->GetEntry(ievent);
958  if(ievent%printstep==0 && ievent!=0) cout<<"Event # "<<ievent<< " # hits "<<prt_event->GetHitSize()<<endl;
959  if(ievent == 0){
960  if(gROOT->GetApplication()){
961  prt_info += "beam test";
962  TIter next(gROOT->GetApplication()->InputFiles());
963  TObjString *os=0;
964  while((os = (TObjString*)next())){
965  prt_info += os->GetString()+" ";
966  }
967  prt_info += "\n";
968  }
969  prt_info += prt_event->PrintInfo();
970  prt_mom = prt_event->GetMomentum().Mag() +0.01;
971  prt_theta = prt_event->GetAngle();
972  prt_phi = prt_event->GetPhi();
973  prt_geometry= prt_event->GetGeometry();
974  prt_beamx= prt_event->GetBeamX();
975  prt_beamz= prt_event->GetBeamZ();
976  prt_test1 = prt_event->GetTest1();
977  prt_test2 = prt_event->GetTest2();
978  }
979  prt_particle = prt_event->GetParticle();
980  if(prt_event->GetParticle()<3000 && prt_event->GetParticle()>0){
981  prt_pid=prt_particleArray[prt_event->GetParticle()];
982  }
983 }
984 #endif
985 
986 int prt_getColorId(int ind, int style =0){
987  int cid = 1;
988  if(style==0) {
989  cid=ind+1;
990  if(cid==5) cid =8;
991  if(cid==3) cid =kOrange+2;
992  }
993  if(style==1) cid=ind+300;
994  return cid;
995 }
996 
997 int prt_shiftHist(TH1 *hist, double double_shift){
998  int bins=hist->GetXaxis()->GetNbins();
999  double xmin=hist->GetXaxis()->GetBinLowEdge(1);
1000  double xmax=hist->GetXaxis()->GetBinUpEdge(bins);
1001  double_shift=double_shift*(bins/(xmax-xmin));
1002  int shift=0;
1003  if(double_shift<0) shift=TMath::FloorNint(double_shift);
1004  if(double_shift>0) shift=TMath::CeilNint(double_shift);
1005  if(shift==0) return 0;
1006  if(shift>0){
1007  for(int i=1; i<=bins; i++){
1008  if(i+shift<=bins) hist->SetBinContent(i,hist->GetBinContent(i+shift));
1009  if(i+shift>bins) hist->SetBinContent(i,0);
1010  }
1011  return 0;
1012  }
1013  if(shift<0){
1014  for(int i=bins; i>0; i--){
1015  if(i+shift>0) hist->SetBinContent(i,hist->GetBinContent(i+shift));
1016  if(i+shift<=0) hist->SetBinContent(i,0);
1017  }
1018  return 0;
1019  }
1020  return 1;
1021 }
1022 
1023 void prt_addInfo(TString str){
1024  prt_info += str+"\n";
1025 }
1026 
1027 void prt_writeInfo(TString filename){
1028  std::ofstream myfile;
1029  myfile.open (filename);
1030  myfile << prt_info+"\n";
1031  myfile.close();
1032 }
1033 
1034 void prt_writeString(TString filename, TString str){
1035  std::ofstream myfile;
1036  myfile.open (filename);
1037  myfile << str+"\n";
1038  myfile.close();
1039  std::cout<<"output: "<<filename<<std::endl;
1040 }
1041 
1042 TString prt_createDir(TString inpath=""){
1043  if(inpath != "") prt_savepath = inpath;
1044  TString finalpath = prt_savepath;
1045 
1046  if(finalpath =="") return "";
1047 
1048  if(prt_savepath.EndsWith("auto")) {
1049  TString dir = prt_savepath.ReplaceAll("auto","data");
1050  gSystem->mkdir(dir);
1051  TDatime *time = new TDatime();
1052  TString path(""), stime = Form("%d.%d.%d", time->GetDay(),time->GetMonth(),time->GetYear());
1053  gSystem->mkdir(dir+"/"+stime);
1054  for(int i=0; i<1000; i++){
1055  path = stime+"/"+Form("arid-%d",i);
1056  if(gSystem->mkdir(dir+"/"+path)==0) break;
1057  }
1058  gSystem->Unlink(dir+"/last");
1059  gSystem->Symlink(path, dir+"/last");
1060  finalpath = dir+"/"+path;
1061  prt_savepath=finalpath;
1062  }else{
1063  gSystem->mkdir(prt_savepath,kTRUE);
1064  }
1065 
1066  prt_writeInfo(finalpath+"/readme");
1067  return finalpath;
1068 }
1069 
1070 void prt_canvasPrint(TPad *c, TString name="", TString path="", int what=0){
1071  c->Modified();
1072  c->Update();
1073  c->Print(path+"/"+name+".png");
1074  if(what>0) c->Print(path+"/"+name+".C");
1075  if(what>1) c->Print(path+"/"+name+".pdf");
1076  if(what>2) c->Print(path+"/"+name+".eps");
1077 }
1078 
1079 void prt_set_style(TCanvas *c){
1080  // prt_setRootPalette(1);
1081  if(fabs(c->GetBottomMargin()-0.1)<0.001) c->SetBottomMargin(0.12);
1082  TIter next(c->GetListOfPrimitives());
1083  TObject *obj;
1084 
1085  while((obj = next())){
1086  if(obj->InheritsFrom("TH1")){
1087  TH1F *hh = (TH1F*)obj;
1088  hh->GetXaxis()->SetTitleSize(0.06);
1089  hh->GetYaxis()->SetTitleSize(0.06);
1090 
1091  hh->GetXaxis()->SetLabelSize(0.05);
1092  hh->GetYaxis()->SetLabelSize(0.05);
1093 
1094  hh->GetXaxis()->SetTitleOffset(0.85);
1095  hh->GetYaxis()->SetTitleOffset(0.76);
1096 
1097  if(c->GetWindowHeight()>700){
1098  hh->GetXaxis()->SetTitleSize(0.045);
1099  hh->GetYaxis()->SetTitleSize(0.045);
1100 
1101  hh->GetXaxis()->SetLabelSize(0.035);
1102  hh->GetYaxis()->SetLabelSize(0.035);
1103 
1104  hh->GetXaxis()->SetTitleOffset(0.85);
1105  hh->GetYaxis()->SetTitleOffset(0.98);
1106  c->SetRightMargin(0.12);
1107  }
1108 
1109  if(c->GetWindowWidth()<=600){
1110  hh->GetXaxis()->SetTitleSize(0.05);
1111  hh->GetYaxis()->SetTitleSize(0.05);
1112 
1113  hh->GetXaxis()->SetLabelSize(0.04);
1114  hh->GetYaxis()->SetLabelSize(0.04);
1115 
1116  hh->GetXaxis()->SetTitleOffset(0.85);
1117  hh->GetYaxis()->SetTitleOffset(0.92);
1118  c->SetRightMargin(0.12);
1119  }
1120 
1121 
1122  if(fabs(c->GetBottomMargin()-0.12)<0.001){
1123  TPaletteAxis *palette = (TPaletteAxis*)hh->GetListOfFunctions()->FindObject("palette");
1124  if(palette) {
1125  palette->SetY1NDC(0.12);
1126  c->Modified();
1127  }
1128  }
1129  c->Modified();
1130  c->Update();
1131  }
1132 
1133  if(obj->InheritsFrom("TGraph")){
1134  TGraph *gg = (TGraph*)obj;
1135  gg->GetXaxis()->SetLabelSize(0.05);
1136  gg->GetXaxis()->SetTitleSize(0.06);
1137  gg->GetXaxis()->SetTitleOffset(0.84);
1138 
1139  gg->GetYaxis()->SetLabelSize(0.05);
1140  gg->GetYaxis()->SetTitleSize(0.06);
1141  gg->GetYaxis()->SetTitleOffset(0.8);
1142  }
1143 
1144  if(obj->InheritsFrom("TMultiGraph")){
1145  TMultiGraph *gg = (TMultiGraph*)obj;
1146  gg->GetXaxis()->SetLabelSize(0.05);
1147  gg->GetXaxis()->SetTitleSize(0.06);
1148  gg->GetXaxis()->SetTitleOffset(0.84);
1149 
1150  gg->GetYaxis()->SetLabelSize(0.05);
1151  gg->GetYaxis()->SetTitleSize(0.06);
1152  gg->GetYaxis()->SetTitleOffset(0.8);
1153  }
1154 
1155  if(obj->InheritsFrom("TF1")){
1156  TF1 *f = (TF1*)obj;
1157  f->SetNpx(500);
1158  }
1159  }
1160 }
1161 
1162 void prt_save(TPad *c= NULL,TString path="", int what=0, int style=0){
1163  TString name = c->GetName();
1164  Bool_t batch = gROOT->IsBatch();
1165  gROOT->SetBatch(1);
1166 
1167  if(c && path != "") {
1168  int w = 800, h = 400;
1169  if(style != -1){
1170  if(style == 1) {w = 800; h = 500;}
1171  if(style == 2) {w = 800; h = 600;}
1172  if(style == 3) {w = 800; h = 400;}
1173  if(style == 5) {w = 800; h = 900;}
1174  if(style == 0){
1175  w = ((TCanvas*)c)->GetWindowWidth();
1176  h = ((TCanvas*)c)->GetWindowHeight();
1177  }
1178 
1179  TCanvas *cc;
1180  if(TString(c->GetName()).Contains("hp") || TString(c->GetName()).Contains("cdigi")) {
1181  cc = prt_drawDigi(prt_last_layoutId,prt_last_maxz,prt_last_minz);
1182  cc->SetCanvasSize(800,400);
1183  if(name.Contains("=")) name = name.Tokenize('=')->First()->GetName();
1184  }else{
1185  cc = new TCanvas(TString(c->GetName())+"exp","cExport",0,0,w,h);
1186  cc = (TCanvas*) c->DrawClone();
1187  cc->SetCanvasSize(w,h);
1188  }
1189 
1190  if(style == 0) prt_set_style(cc);
1191 
1192  prt_canvasPrint(cc,name,path,what);
1193  }else{
1194  c->SetCanvasSize(w,h);
1195  prt_canvasPrint(c,name,path,what);
1196  }
1197  }
1198 
1199  gROOT->SetBatch(batch);
1200 }
1201 
1202 TString prt_createSubDir(TString dir="dir"){
1203  gSystem->mkdir(dir);
1204  return dir;
1205 }
1206 
1207 TList *prt_canvaslist;
1208 
1209 TCanvas *prt_canvasGet(TString name="c"){
1210  TIter next(prt_canvaslist);
1211  TCanvas *c = 0;
1212  while((c = (TCanvas*) next())){
1213  if(c->GetName()==name || name=="*") break;
1214  }
1215  return c;
1216 }
1217 
1218 void prt_canvasAdd(TString name="c",int w=800, int h=400){
1219  if(!prt_canvasGet(name)){
1220  if(!prt_canvaslist) prt_canvaslist = new TList();
1221  TCanvas *c = new TCanvas(name,name,0,0,w,h);
1222  prt_canvaslist->Add(c);
1223  }
1224 }
1225 
1226 void prt_canvasAdd(TCanvas *c){
1227  if(!prt_canvaslist) prt_canvaslist = new TList();
1228  c->cd();
1229  prt_canvaslist->Add(c);
1230 }
1231 
1232 void prt_canvasCd(){
1233 
1234 }
1235 
1236 void prt_canvasDel(TString name="c"){
1237  TIter next(prt_canvaslist);
1238  TCanvas *c=0;
1239  while((c = (TCanvas*) next())){
1240  if(c->GetName()==name || name=="*") prt_canvaslist->Remove(c);
1241  c->Delete();
1242  }
1243 }
1244 
1245 // style = 0 - for web blog
1246 // style = 1 - for talk
1247 // what = 0 - save in png
1248 // what = 1 - save in png, C
1249 // what = 2 - save in png, C, pdf
1250 // what = 3 - save in png, C, pdf, eps
1251 void prt_canvasSave(int what=1, int style=0, Bool_t rm=false){
1252  TIter next(prt_canvaslist);
1253  TCanvas *c=0;
1254  TString path = prt_createDir();
1255  while((c = (TCanvas*) next())){
1256  prt_set_style(c);
1257  prt_save(c, path, what,style);
1258  if(rm){
1259  prt_canvaslist->Remove(c);
1260  c->Close();
1261  }
1262  }
1263 }
1264 
1265 // path - folder for saving
1266 void prt_canvasSave(TString path, int what=1, int style=0, Bool_t rm=false){
1267  prt_savepath = path;
1268  prt_canvasSave(what,style,rm);
1269 }
1270 
1271 void prt_set_style(){
1272  TIter next(prt_canvaslist);
1273  TCanvas *c=0;
1274  TString path = prt_createDir();
1275  while((c = (TCanvas*) next())){
1276  prt_set_style(c);
1277  }
1278 }
1279 
1280 void prt_waitPrimitive(TString name, TString prim=""){
1281  TIter next(prt_canvaslist);
1282  TCanvas *c=0;
1283  while((c = (TCanvas*) next())){
1284  if(TString(c->GetName())==name){
1285  c->Modified();
1286  c->Update();
1287  c->WaitPrimitive(prim);
1288  }
1289  }
1290 }
1291 
1292 double prt_integral(TH1F *h,double xmin, double xmax){
1293  TAxis *axis = h->GetXaxis();
1294  int bmin = axis->FindBin(xmin);
1295  int bmax = axis->FindBin(xmax);
1296  double integral = h->Integral(bmin,bmax);
1297  integral -= h->GetBinContent(bmin)*(xmin-axis->GetBinLowEdge(bmin))/axis->GetBinWidth(bmin);
1298  integral -= h->GetBinContent(bmax)*(axis->GetBinUpEdge(bmax)-xmax)/axis->GetBinWidth(bmax);
1299  return integral;
1300 }
1301 
1302 void prt_normalize(TH1F* hists[],int size){
1303  // for(int i=0; i<size; i++){
1304  // hists[i]->Scale(1/hists[i]->Integral(), "width");
1305  // }
1306 
1307  double max = 0;
1308  double min = 0;
1309  for(int i=0; i<size; i++){
1310  double tmax = hists[i]->GetBinContent(hists[i]->GetMaximumBin());
1311  double tmin = hists[i]->GetMinimum();
1312  if(tmax>max) max = tmax;
1313  if(tmin<min) min = tmin;
1314  }
1315  max += 0.05*max;
1316  for(int i=0; i<size; i++){
1317  hists[i]->GetYaxis()->SetRangeUser(min,max);
1318  }
1319 }
1320 
1321 void prt_normalizeto(TH1F* hists[],int size, double max=1){
1322  for(int i=0; i<size; i++){
1323  double tmax = hists[i]->GetBinContent(hists[i]->GetMaximumBin());
1324  if(tmax>0)hists[i]->Scale(max/tmax);
1325  }
1326 }
1327 
1328 void prt_normalize(TH1F* h1,TH1F* h2){
1329  double max = (h1->GetMaximum()>h2->GetMaximum())? h1->GetMaximum() : h2->GetMaximum();
1330  max += max*0.1;
1331  h1->GetYaxis()->SetRangeUser(0,max);
1332  h2->GetYaxis()->SetRangeUser(0,max);
1333 }
1334 
1335 // just x for now
1336 TGraph* prt_smooth(TGraph* g,int smoothness=1){
1337  double x, y;
1338  int n = g->GetN();
1339  TH1F *h = new TH1F("h","h",g->GetN(),0,n);
1340  TGraph *gr = new TGraph();
1341  gr->SetName(g->GetName());
1342  for(auto i=0; i<n; i++){
1343  g->GetPoint(i,x,y);
1344  h->Fill(i,x);
1345  }
1346 
1347  h->Smooth(smoothness);
1348 
1349  for(auto i=0;i<n;i++){
1350  g->GetPoint(i,x,y);
1351  gr->SetPoint(i,h->GetBinContent(i),y);
1352  }
1353  return gr;
1354 }
1355 
1356 int prt_get_pid(int pdg){
1357  int pid=0;
1358  if(pdg==11) pid=0; //e
1359  if(pdg==13) pid=1; //mu
1360  if(pdg==211) pid=2; //pi
1361  if(pdg==321) pid=3; //K
1362  if(pdg==2212) pid=4; //p
1363  return pid;
1364 }
1365 
1366 double prt_get_momentum_from_tof(double dist,double dtof){
1367  double s = dtof*0.299792458/dist;
1368  double x = s*s;
1369  double a = prt_mass[2]*prt_mass[2]; //pi
1370  double b = prt_mass[4]*prt_mass[4]; //p
1371  double p = sqrt((a - 2*sqrt(a*a+a*b*x-2*a*b+b*b)/s + b)/(x - 4));
1372 
1373  return p;
1374 }
1375 
1376 // return TOF difference [ns] for mominum p [GeV] and flight path l [m]
1377 double prt_get_tof_diff(int pid1=211, int pid2=321, double p=1, double l=2){
1378  double c = 299792458;
1379  double m1 = prt_mass[prt_get_pid(pid1)];
1380  double m2 = prt_mass[prt_get_pid(pid2)];
1381  double td = l*(sqrt(p*p+m1*m1)-sqrt(p*p+m2*m2))/(p*c)*1E9;
1382 
1383  // relativistic
1384  // td = l*(m1*m1-m2*m2)/(2*p*p*c)*1E9;
1385 
1386  return td;
1387 }
1388 
1389 bool prt_ispath(TString path){
1390  Long_t *id(0),*size(0),*flags(0),*modtime(0);
1391  return !gSystem->GetPathInfo(path,id,size,flags,modtime);
1392 }
1393 
1394 int prt_get3digit(TString str){
1395  TPRegexp e("[0-9][0-9][0-9]");
1396  return ((TString)str(e)).Atoi();
1397 }
1398 
1399 #endif