Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-30 10:30:56 

0001 // prttools - useful functions 
0002 // original author: Roman Dzhygadlo - GSI Darmstadt 
0003 
0004 #ifndef prttools_h
0005 #define prttools_h 1
0006 
0007 #include "TROOT.h"
0008 #include "TSystem.h"
0009 #include "TStyle.h"
0010 #include "TCanvas.h"
0011 #include "TPad.h"
0012 #include "TH1.h"
0013 #include "TH2.h"
0014 #include "TGraph.h"
0015 #include "TMultiGraph.h"
0016 #include "TSpline.h"
0017 #include "TF1.h"
0018 #include "TFile.h"
0019 #include "TTree.h"
0020 #include "TClonesArray.h"
0021 #include "TVector3.h"
0022 #include "TMath.h"
0023 #include "TChain.h"
0024 #include "TGaxis.h"
0025 #include "TColor.h"
0026 #include "TString.h"
0027 #include "TArrayD.h"
0028 //#include "TSpectrum.h"
0029 //#include "TSpectrum2.h"
0030 #include "Math/Minimizer.h"
0031 #include "Math/Factory.h"
0032 #include "Math/Functor.h"
0033 #include "TRandom2.h"
0034 #include "TError.h"
0035 #include "TPaveStats.h"
0036 #include "TObjString.h"
0037 #include "TApplication.h"
0038 #include <TLegend.h>
0039 #include <TAxis.h>
0040 #include <TPaletteAxis.h>
0041 #include <TRandom.h>
0042 #include <TCutG.h>
0043 #include <TKey.h>
0044 #include "TPRegexp.h"
0045 #include "TFitResult.h"
0046 
0047 #include <iostream>
0048 #include <fstream>
0049 #include <sstream>
0050 
0051 const int prt_nmcp = 24;
0052 const int prt_npix = 16 * 16; 
0053 
0054 const int prt_maxdircch = prt_nmcp*prt_npix;
0055 const int prt_maxnametdc=10000;
0056 
0057 TRandom  prt_rand;
0058 TChain*  prt_ch = 0;
0059 int    prt_entries(0),prt_particle(0),prt_geometry(2023),prt_beamx(0),prt_beamz(0);
0060 int    prt_last_layoutId,prt_last_max,prt_last_maxz,prt_last_minz;
0061 double prt_theta(0), prt_test1(0),prt_test2(0),prt_mom(0),prt_phi(0);
0062 TString  prt_savepath(""),prt_info("");
0063 TH2F*    prt_hdigi[prt_nmcp];
0064 //TSpectrum *prt_spect = new TSpectrum(2);
0065 
0066 const int prt_ntdcm = 80; //41
0067 int prt_ntdc = prt_ntdcm;
0068 int prt_maxch = prt_ntdc*48;
0069 TString prt_tdcsid[prt_ntdcm];
0070 const int prt_maxchm = prt_ntdcm*48;
0071 int map_tdc[prt_maxnametdc];
0072 int map_mpc[prt_maxchm/prt_npix][prt_npix];
0073 int map_mcp[prt_maxchm];
0074 int map_pix[prt_maxchm];
0075 int map_row[prt_maxchm];
0076 int map_col[prt_maxchm];
0077 
0078 int prt_pid(0), prt_pdg[]={-11,13,211,321,2212};
0079 double prt_mass[]={0.000511,0.1056584,0.139570,0.49368,0.9382723};
0080 TString  prt_name[]={"e","muon","pion","kaon","proton"};
0081 TString  prt_lname[]={"e","#mu","#pi","K","p"};
0082 int    prt_color[]={kOrange+6,kCyan+1,kBlue+1,kRed+1,kRed+1};
0083 double prt_particleArray[3000];
0084 
0085 TF1 *prt_gaust;
0086 /*TVector3 prt_fit(TH1 *h, double range = 3, double threshold=20, double limit=2, int peakSearch=1,int bkg = 0, TString opt="MQ"){
0087   int binmax = h->GetMaximumBin();
0088   double xmax = h->GetXaxis()->GetBinCenter(binmax);
0089   if(bkg==0) prt_gaust = new TF1("prt_gaust","[0]*exp(-0.5*((x-[1])/[2])^2)",xmax-range,xmax+range);
0090   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);
0091   prt_gaust->SetNpx(500);
0092   prt_gaust->SetParNames("const","mean","sigma");
0093   prt_gaust->SetLineColor(2);
0094   double integral = h->Integral(h->GetXaxis()->FindBin(xmax-range),h->GetXaxis()->FindBin(xmax+range));
0095   double xxmin, xxmax, sigma1(0), mean1(0), mean2(0);
0096   xxmax = xmax;
0097   xxmin = xxmax;
0098   int nfound(1);
0099   if(integral>threshold){
0100     
0101     if(peakSearch == 1){
0102       prt_gaust->SetParameter(1,xmax);
0103       prt_gaust->SetParameter(2,0.005);
0104       prt_gaust->SetParLimits(2,0.003,limit);
0105       h->Fit("prt_gaust",opt,"",xxmin-range, xxmax+range);
0106     }
0107     
0108     if(peakSearch>1){
0109       nfound = prt_spect->Search(h,4,"goff",0.1);
0110       std::cout<<"nfound  "<<nfound <<std::endl;
0111       if(nfound==1){
0112     prt_gaust =new TF1("prt_gaust","gaus(0)",xmax-range,xmax+range);
0113     prt_gaust->SetNpx(500);
0114     prt_gaust->SetParameter(1,prt_spect->GetPositionX()[0]);
0115       }else if(nfound>=2){
0116     double p1 = prt_spect->GetPositionX()[0];
0117     double p2 = prt_spect->GetPositionX()[1];
0118     if(p1>p2) {
0119       xxmax = p1;
0120       xxmin = p2;
0121     }else {
0122       xxmax = p2;
0123       xxmin = p1;
0124     }
0125     if(peakSearch==20){
0126       xxmax=xxmin;
0127       prt_gaust =new TF1("prt_gaust","gaus(0)",xxmin-range,xxmin+range);
0128       prt_gaust->SetNpx(500);
0129       prt_gaust->SetParameter(1,prt_spect->GetPositionX()[0]);
0130     }else{
0131       prt_gaust =new TF1("prt_gaust","gaus(0)+gaus(3)",xmax-range,xmax+range);
0132       prt_gaust->SetNpx(500);
0133       prt_gaust->SetParameter(0,1000);
0134       prt_gaust->SetParameter(3,1000);
0135     
0136       prt_gaust->FixParameter(1,xxmin);
0137       prt_gaust->FixParameter(4,xxmax);
0138       prt_gaust->SetParameter(2,0.1);
0139       prt_gaust->SetParameter(5,0.1);
0140       h->Fit("prt_gaust",opt,"",xxmin-range, xxmax+range);
0141       prt_gaust->ReleaseParameter(1);
0142       prt_gaust->ReleaseParameter(4);
0143     }
0144       }
0145     
0146       prt_gaust->SetParameter(2,0.2);
0147       prt_gaust->SetParameter(5,0.2);
0148     }
0149 
0150     h->Fit("prt_gaust",opt,"",xxmin-range, xxmax+range);
0151     mean1 = prt_gaust->GetParameter(1);
0152     sigma1 = prt_gaust->GetParameter(2);
0153     if(sigma1>10) sigma1=10;
0154     
0155     if(peakSearch == 2){
0156       mean2 = (nfound==1) ? prt_gaust->GetParameter(1) : prt_gaust->GetParameter(4);
0157       // sigma2 = (nfound==1) ? prt_gaust->GetParameter(2) : prt_gaust->GetParameter(5);
0158     }
0159   }
0160   delete prt_gaust;
0161   return TVector3(mean1,sigma1,mean2);
0162   }
0163 
0164 TGraph *prt_fitslices(TH2F *hh,double minrange=0, double maxrange=0, double fitrange=1,int rebin=1,int ret=0){
0165   TH2F *h =(TH2F*) hh->Clone("h");
0166   h->RebinY(rebin);
0167   int point(0);
0168   TGraph *gres = new TGraph();
0169   for (int i=1;i<h->GetNbinsY();i++){
0170     double x = h->GetYaxis()->GetBinCenter(i);
0171     TH1D* hp;
0172     if(minrange!=maxrange){
0173       TCutG *cut = new TCutG("prt_onepeakcut",5);
0174       cut->SetVarX("y");
0175       cut->SetVarY("x");
0176       cut->SetPoint(0,minrange,-1E6);
0177       cut->SetPoint(1,minrange, 1E6);
0178       cut->SetPoint(2,maxrange, 1E6);
0179       cut->SetPoint(3,maxrange,-1E6);
0180       cut->SetPoint(4,minrange,-1E6);
0181     
0182       hp = h->ProjectionX(Form("bin%d",i),i,i,"[prt_onepeakcut]");
0183     }else{
0184       hp = h->ProjectionX(Form("bin%d",i),i,i);
0185     }
0186 
0187     TVector3 res = prt_fit((TH1F*)hp,fitrange,100,2,1,1);
0188     double y=0;
0189     if(ret==0) y = res.X();
0190     if(ret==1) y = res.Y();
0191     if(ret==2) y = res.X() + 0.5*res.Y();
0192     if(ret==3) y = res.X() - 0.5*res.Y();
0193     if(y==0 || y<minrange || y>maxrange) continue;
0194     
0195     gres->SetPoint(point,y,x);
0196     gres->SetLineWidth(2);
0197     gres->SetLineColor(kRed);
0198     point++;
0199   }
0200   return gres;
0201   }*/
0202 
0203 void prt_createMap(int setupid = 2019){
0204   prt_geometry = setupid;
0205 
0206   TGaxis::SetMaxDigits(4);
0207   for(int i=0; i<prt_maxnametdc; i++) map_tdc[i]=-1;
0208   for(int i=0; i<prt_ntdc; i++){
0209     int dec = TString::BaseConvert(prt_tdcsid[i],16,10).Atoi();
0210     map_tdc[dec]=i;
0211   }
0212   
0213   //  for(int ch=0; ch<prt_maxdircch; ch++){
0214   for(int ch=0; ch<prt_maxch; ch++){
0215     int mcp = ch/prt_npix;
0216     int pix = ch%prt_npix;
0217     int col = pix/2 - 8*(pix/16);
0218     int row = pix%2 + 2*(pix/16);
0219     // pix = col+sqrt(prt_npix)*row;
0220 
0221     map_mpc[mcp][pix]=ch;
0222     map_mcp[ch] = mcp;
0223     map_pix[ch] = pix;
0224     map_row[ch] = row;
0225     map_col[ch] = col;
0226   } 
0227 
0228   for(int i=0; i<5; i++){
0229     prt_particleArray[prt_pdg[i]]=i;
0230   }
0231   prt_particleArray[212]=2;
0232 }
0233 
0234 TString prt_randstr(int len = 10){
0235   TString str = ""; 
0236   static const char alphanum[] =
0237     "0123456789"
0238     "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
0239     "abcdefghijklmnopqrstuvwxyz";
0240 
0241   for (int i = 0; i < len; ++i) {
0242     str += alphanum[rand() % (sizeof(alphanum) - 1)];
0243   }
0244   return str;
0245 }
0246 
0247 // layoutId == 5    - 5 row's design for the PANDA Barrel DIRC
0248 // layoutId == 2015 - cern 2015
0249 // layoutId == 2016 - cern 2016
0250 // layoutId == 2017 - cern 2017
0251 // layoutId == 2018 - cern 2018
0252 // layoutId == 2021 - new 3.6 row's design for the PANDA Barrel DIRC
0253 // layoutId == 2023 - new 2x4 layout for the PANDA Barrel DIRC
0254 // layoutId == 2031 - EIC DIRC beam test
0255 // layoutId == 2030 - EIC DIRC prism
0256 // layoutId == 2032 - EIC DIRC focusing prism 
0257 TH1 * prt_cdigi_th;
0258 TCanvas *prt_drawDigi(int layoutId = 0, double maxz = 0, double minz = 0, TCanvas *cdigi = NULL){
0259 
0260   prt_last_layoutId=layoutId;
0261   prt_last_maxz=maxz;
0262   prt_last_minz=minz;
0263 
0264   if(prt_geometry==2021) layoutId=2021;
0265   // if(prt_geometry==2019) layoutId=2018;
0266 
0267   TString sid = prt_randstr(3);
0268   if(cdigi) cdigi->cd();
0269   else cdigi = new TCanvas("hp="+sid,"hp_"+sid,800,400);
0270 
0271   TPad* prt_hpads[prt_nmcp];
0272   TPad* prt_hpglobal;
0273   
0274   if(layoutId==2015 ||  layoutId==5) prt_hpglobal = new TPad(sid,"T",0.04,0.04,0.88,0.96);
0275   else if(layoutId==2021) prt_hpglobal = new TPad(sid,"T",0.12,0.02,0.78,0.98);
0276   else if(layoutId==2016) prt_hpglobal = new TPad(sid,"T",0.2,0.02,0.75,0.98);
0277   else if(layoutId==2017) prt_hpglobal = new TPad(sid,"T",0.15,0.02,0.80,0.98);
0278   else if(layoutId==2018) prt_hpglobal = new TPad(sid,"T",0.05,0.07,0.9,0.93);
0279   else if(layoutId==2023) prt_hpglobal = new TPad(sid,"T",0.073,0.02,0.877,0.98);
0280   else if(layoutId==2030) prt_hpglobal = new TPad(sid,"T",0.10,0.025,0.82,0.975);
0281   else if(layoutId==2032) prt_hpglobal = new TPad(sid,"T",0.04,0.025,0.91,0.975);
0282   else if(layoutId==2031) prt_hpglobal = new TPad(sid,"T",0.12,0.01,0.80,0.99);
0283   else  prt_hpglobal = new TPad(sid,"T",0.04,0.04,0.96,0.96);
0284 
0285   prt_hpglobal->SetFillStyle(0);
0286   prt_hpglobal->Draw();  
0287   prt_hpglobal->cd();
0288   
0289   int nrow = 3, ncol = 5;
0290   if(layoutId ==2016) ncol=3;
0291   if(layoutId ==2017) ncol=4;
0292   if(layoutId ==2018 || layoutId ==2023) {nrow=2; ncol=4;}
0293   if(layoutId ==2021) ncol=4;
0294   if(layoutId ==2030) {nrow=4; ncol=6;}
0295   if(layoutId ==2032) {nrow=4; ncol=7;}
0296   if(layoutId ==2031) {nrow=3; ncol=4;}
0297   
0298   if(layoutId > 1){
0299     float tbw(0.02), tbh(0.01), shift(0),shiftw(0.02),shifth(0),margin(0.01);
0300     int padi(0);
0301       
0302     for(int i=0; i<ncol; i++){
0303       for(int j=0; j<nrow; j++){
0304     if(j==1) shift = -0.028;
0305     else shift = 0;
0306     shifth=0;
0307     if(layoutId == 5) {shift =0; shiftw=0.001; tbw=0.001; tbh=0.001;}
0308     if(layoutId == 2021) {
0309       if(i==0 && j == nrow-1) continue;
0310       shift =0; shiftw=0.001; tbw=0.001; tbh=0.001;
0311       if(i==0) shifth=0.167;
0312     }
0313     if(layoutId == 2016) {
0314       shift = -0.01; shiftw=0.01; tbw=0.03; tbh=0.006;
0315       if(j==1) shift += 0.015;
0316     }
0317     if(layoutId == 2017) {
0318       margin= 0.1;
0319       shift = 0; shiftw=0.01; tbw=0.005; tbh=0.006;
0320       //if(j==1) shift += 0.015;
0321     }
0322     if(layoutId == 2018) {
0323       margin= 0.1;
0324       shift = 0; shiftw=0.01; tbw=0.005; tbh=0.006;
0325     }
0326     if(layoutId == 2023) {
0327       margin= 0.1;
0328       shift = 0; shiftw=0.01; tbw=0.0015; tbh=0.042;
0329     }
0330     if(layoutId == 2030) {
0331       margin= 0.1;
0332       shift = 0; shiftw=0.01; tbw=0.001; tbh=0.001;
0333       padi=j*ncol+i;
0334     }
0335     if(layoutId == 2032) {
0336       margin= 0.1;
0337       shift = 0; shiftw=0.01; tbw=0.001; tbh=0.001;
0338       padi=j*ncol+i;
0339     }
0340     if(layoutId == 2031) {
0341       margin= 0.1;
0342       shift = 0; shiftw=0.01; tbw=0.001; tbh=0.001;
0343       padi=i*nrow+j;
0344     }
0345 
0346     prt_hpads[padi] =  new TPad(sid+Form("P%d",i*10+j),"T",
0347                     i/(ncol+2*margin)+tbw+shift+shiftw,
0348                     j/(double)nrow+tbh+shifth,
0349                     (i+1)/(ncol+2*margin)-tbw+shift+shiftw,
0350                     (1+j)/(double)nrow-tbh+shifth, 21);
0351     prt_hpads[padi]->SetFillColor(kCyan-8);
0352     prt_hpads[padi]->SetMargin(0.055,0.055,0.055,0.055);
0353     prt_hpads[padi]->Draw();
0354     padi++;
0355       }
0356     }
0357   }else{
0358     float tbw(0.02), tbh(0.01), shift(0),shiftw(-0.02);
0359     int padi(0);
0360     for(int ii=0; ii<ncol; ii++){
0361       for(int j=0; j<nrow; j++){
0362     if(j==1) shift = 0.04;
0363     else shift = 0;
0364     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);
0365     prt_hpads[padi]->SetFillColor(kCyan-8);
0366     prt_hpads[padi]->SetMargin(0.04,0.04,0.04,0.04);
0367     prt_hpads[padi]->Draw(); 
0368     padi++;
0369       }
0370     }
0371   }
0372 
0373   int np;
0374   double max=0;
0375 
0376   {
0377     double tmax;
0378     if(maxz==0){
0379       for(int p=0; p<nrow*ncol;p++){
0380     tmax = prt_hdigi[p]->GetBinContent(prt_hdigi[p]->GetMaximumBin());
0381     if(max<tmax) max = tmax;
0382       }
0383     }else{
0384       max = maxz;
0385     }
0386   
0387     if(maxz==-2 || minz==-2){ // optimize range
0388       for(int p=0; p<nrow*ncol;p++){
0389     tmax = prt_hdigi[p]->GetMaximum();
0390     if(max<tmax) max = tmax;
0391       }
0392       int tbins = 2000;
0393       TH1F *h = new TH1F("","",tbins,0,max);
0394       for(int p=0; p<nrow*ncol;p++){
0395     for(int i=0; i<8; i++){
0396       for(int j=0; j<8; j++){
0397         double val = prt_hdigi[p]->GetBinContent(i+1,j+1);
0398         if(val!=0) h->Fill(val);
0399 
0400       }
0401     }
0402       }
0403       double integral;
0404       for(int i=0; i<tbins; i++){
0405     integral = h->Integral(0,i);
0406     if(integral>0) {
0407       if(minz==-2) minz = h->GetBinCenter(i);
0408       break;
0409     } 
0410       }
0411 
0412       for(int i=tbins; i>0; i--){
0413     integral = h->Integral(i,tbins);
0414     if(integral>10) {
0415       if(maxz==-2) max = h->GetBinCenter(i);
0416       break;
0417     } 
0418       }
0419     }
0420   }
0421 
0422   prt_last_max = max;
0423   int nnmax(0);
0424   for(int p=0; p<nrow*ncol;p++){
0425     if(layoutId == 1 || layoutId == 4)  np =p%nrow*ncol + p/3;
0426     else np = p;
0427 
0428     if(layoutId == 6 && p>10) continue;
0429     
0430     prt_hpads[p]->cd();
0431     prt_hdigi[np]->Draw("col");//"col+text"
0432     if(maxz==-1)  max = prt_hdigi[np]->GetBinContent(prt_hdigi[np]->GetMaximumBin());
0433     if(nnmax<prt_hdigi[np]->GetEntries()) nnmax=np;
0434     prt_hdigi[np]->SetMaximum(max);
0435     prt_hdigi[np]->SetMinimum(minz);
0436   }
0437   
0438   // prt_cdigi_palette = (TPaletteAxis*)prt_hdigi[nnmax]->GetListOfFunctions()->FindObject("palette");
0439   // prt_cdigi_palette->SetX1NDC(0.89);
0440   // prt_cdigi_palette->SetY1NDC(0.1);
0441   // prt_cdigi_palette->SetX2NDC(0.93);
0442   // prt_cdigi_palette->SetY2NDC(0.90);
0443   
0444   cdigi->cd();
0445   TPaletteAxis* prt_cdigi_palette;  
0446   if(layoutId==2018 || layoutId==2023)  prt_cdigi_palette = new TPaletteAxis(0.89,0.1,0.93,0.90,(TH1 *) prt_hdigi[nnmax]);
0447   else if(layoutId==2032) prt_cdigi_palette = new TPaletteAxis(0.91,0.1,0.94,0.90,(TH1 *) prt_hdigi[nnmax]);
0448   else prt_cdigi_palette = new TPaletteAxis(0.82,0.1,0.86,0.90,(TH1 *) prt_hdigi[nnmax]);
0449   prt_cdigi_palette->SetName("prt_palette");  
0450   prt_cdigi_palette->Draw();
0451 
0452   cdigi->Modified();
0453   cdigi->Update();
0454   
0455   return cdigi;
0456 }
0457 
0458 
0459 TString prt_getPixData(TString s="m,p,v\n", int layoutId=1){
0460   int nrow = 3, ncol = 5, np, nmax=0, npix=8;
0461   if(layoutId ==2016) ncol=3;
0462   if(layoutId ==2017) ncol=4;
0463   if(layoutId ==2018 || layoutId ==2023) {nrow=2; ncol=4;}
0464   if(layoutId ==2021) ncol=4;
0465   if(layoutId ==2030) {nrow=4; ncol=6; npix=16;}
0466   if(layoutId ==2032) {nrow=4; ncol=7; npix=16;}
0467   if(layoutId ==2031) {nrow=3; ncol=4; npix=16;}
0468   
0469   int nnmax(0);  
0470   for(int p=0; p<nrow*ncol;p++){
0471     if(layoutId == 1 || layoutId == 4)  np =p%nrow*ncol + p/3;
0472     else if(layoutId == 2030)  np =p%ncol*nrow + p/ncol;
0473     else np = p;
0474 
0475     if(prt_last_maxz==-1) prt_last_max = prt_hdigi[p]->GetBinContent(prt_hdigi[p]->GetMaximumBin());
0476     if(nnmax<prt_hdigi[p]->GetEntries()) nnmax=p;
0477     prt_hdigi[p]->SetMaximum(prt_last_max);
0478     prt_hdigi[p]->SetMinimum(prt_last_minz);
0479     for(int i=1; i<=npix; i++){
0480       for(int j=1; j<=npix; j++){
0481     double weight = (double)(prt_hdigi[p]->GetBinContent(i,j))/(double)prt_last_max*255;
0482     if(weight>255) weight=255;
0483     if(weight > 0) s += Form("%d,%d,%d\n", np, (i-1)*npix+j-1, (int)weight);
0484       }
0485     }
0486   }
0487   return s;
0488 }
0489 
0490 void prt_initDigi(int type=1){  
0491   if(type == 1){
0492     for(int m=0; m<prt_nmcp;m++){
0493       if(prt_hdigi[m]) prt_hdigi[m]->Reset("M");
0494       else{
0495     prt_hdigi[m] = new TH2F( Form("mcp%d", m),Form("mcp%d", m),8,0.,8.,8,0.,8.);
0496     prt_hdigi[m]->SetStats(0);
0497     prt_hdigi[m]->SetTitle(0);
0498     prt_hdigi[m]->GetXaxis()->SetNdivisions(10);
0499     prt_hdigi[m]->GetYaxis()->SetNdivisions(10);
0500     prt_hdigi[m]->GetXaxis()->SetLabelOffset(100);
0501     prt_hdigi[m]->GetYaxis()->SetLabelOffset(100);
0502     prt_hdigi[m]->GetXaxis()->SetTickLength(1);
0503     prt_hdigi[m]->GetYaxis()->SetTickLength(1);
0504     prt_hdigi[m]->GetXaxis()->SetAxisColor(15);
0505     prt_hdigi[m]->GetYaxis()->SetAxisColor(15);
0506       }
0507     }
0508   }
0509   if(type == 2){ //eic
0510     for(int m=0; m<prt_nmcp;m++){
0511       if(prt_hdigi[m]) prt_hdigi[m]->Reset("M");
0512       else{
0513     prt_hdigi[m] = new TH2F( Form("mcp%d", m),Form("mcp%d", m),16,0,16,16,0,16);
0514     prt_hdigi[m]->SetStats(0);
0515     prt_hdigi[m]->SetTitle(0);
0516     prt_hdigi[m]->GetXaxis()->SetNdivisions(20);
0517     prt_hdigi[m]->GetYaxis()->SetNdivisions(20);
0518     prt_hdigi[m]->GetXaxis()->SetLabelOffset(100);
0519     prt_hdigi[m]->GetYaxis()->SetLabelOffset(100);
0520     prt_hdigi[m]->GetXaxis()->SetTickLength(1);
0521     prt_hdigi[m]->GetYaxis()->SetTickLength(1);
0522     prt_hdigi[m]->GetXaxis()->SetAxisColor(15);
0523     prt_hdigi[m]->GetYaxis()->SetAxisColor(15);
0524       }
0525     }
0526   }
0527 }
0528 
0529 void prt_resetDigi(){
0530   for(int m=0; m<prt_nmcp;m++){ 
0531     prt_hdigi[m]->Reset("M");
0532   }
0533 }
0534 
0535 void prt_axisHits800x500(TH2 * hist){
0536   hist->SetStats(0);
0537   hist->SetTitle(Form("%d hits",(int)hist->GetEntries()));
0538   hist->GetXaxis()->SetTitle("z, [cm]");
0539   hist->GetXaxis()->SetTitleSize(0.05);
0540   hist->GetXaxis()->SetTitleOffset(0.8);
0541   hist->GetYaxis()->SetTitle("y, [cm]");
0542   hist->GetYaxis()->SetTitleSize(0.05);
0543   hist->GetYaxis()->SetTitleOffset(0.7);
0544 }
0545 
0546 void prt_axisAngle800x500(TH2 * hist){
0547   hist->SetStats(0);
0548   hist->SetTitle(Form("%d hits",(int)hist->GetEntries()));
0549   hist->GetXaxis()->SetTitle("#theta, [degree]");
0550   hist->GetXaxis()->SetTitleSize(0.05);
0551   hist->GetXaxis()->SetTitleOffset(0.8);
0552   hist->GetYaxis()->SetTitle("photons per track, [#]");
0553   hist->GetYaxis()->SetTitleSize(0.05);
0554   hist->GetYaxis()->SetTitleOffset(0.7);
0555 }
0556 
0557 void prt_axisAngle800x500(TH1 * hist){
0558   hist->SetStats(0);
0559   hist->SetTitle(Form("%d hits",(int)hist->GetEntries()));
0560   hist->GetXaxis()->SetTitle("#theta, [degree]");
0561   hist->GetXaxis()->SetTitleSize(0.05);
0562   hist->GetXaxis()->SetTitleOffset(0.8);
0563   hist->GetYaxis()->SetTitle("photons per track, [#]");
0564   hist->GetYaxis()->SetTitleSize(0.05);
0565   hist->GetYaxis()->SetTitleOffset(0.7);
0566 }
0567 
0568 void prt_axisTime800x500(TH2 * hist){
0569   hist->GetXaxis()->SetTitle("time, [ns]");
0570   hist->GetXaxis()->SetTitleSize(0.05);
0571   hist->GetXaxis()->SetTitleOffset(0.8);
0572   hist->GetYaxis()->SetTitle("entries, #");
0573   hist->GetYaxis()->SetTitleSize(0.05);
0574   hist->GetYaxis()->SetTitleOffset(0.7);
0575   hist->SetLineColor(1);
0576 }
0577 
0578 void prt_axisTime800x500(TH1 * hist, TString xtitle = "time [ns]"){
0579   TGaxis::SetMaxDigits(3);
0580   hist->GetXaxis()->SetTitle(xtitle);
0581   hist->GetXaxis()->SetTitleSize(0.06);
0582   hist->GetXaxis()->SetTitleOffset(0.8);
0583   hist->GetXaxis()->SetLabelSize(0.05);
0584   hist->GetYaxis()->SetTitle("entries [#]");
0585   hist->GetYaxis()->SetTitleSize(0.06);
0586   hist->GetYaxis()->SetTitleOffset(0.7);
0587   hist->GetYaxis()->SetLabelSize(0.05);
0588   hist->SetLineColor(1);
0589 }
0590 
0591 void prt_setPrettyStyle(){
0592   // Canvas printing details: white bg, no borders.
0593   gStyle->SetCanvasColor(0);
0594   gStyle->SetCanvasBorderMode(0);
0595   gStyle->SetCanvasBorderSize(0);
0596 
0597   // Canvas frame printing details: white bg, no borders.
0598   gStyle->SetFrameFillColor(0);
0599   gStyle->SetFrameBorderMode(0);
0600   gStyle->SetFrameBorderSize(0);
0601 
0602   // Plot title details: centered, no bg, no border, nice font.
0603   gStyle->SetTitleX(0.1);
0604   gStyle->SetTitleW(0.8);
0605   gStyle->SetTitleBorderSize(0);
0606   gStyle->SetTitleFillColor(0);
0607 
0608   // Font details for titles and labels.
0609   gStyle->SetTitleFont(42, "xyz");
0610   gStyle->SetTitleFont(42, "pad");
0611   gStyle->SetLabelFont(42, "xyz");
0612   gStyle->SetLabelFont(42, "pad");
0613 
0614   // Details for stat box.
0615   gStyle->SetStatColor(0);
0616   gStyle->SetStatFont(42);
0617   gStyle->SetStatBorderSize(1);
0618   gStyle->SetStatX(0.975);
0619   gStyle->SetStatY(0.9);
0620 
0621   // gStyle->SetOptStat(0);
0622 }
0623 
0624 void prt_setGStyle(TGraph *g, int id){
0625   int prt_coll[]={kBlack,kRed+1,kGreen,  kBlue,  4,kCyan-6,kOrange,  7,8,9,10,1,1,1};
0626   int prt_colm[]={kBlack,kRed+1,kGreen+2,kBlue+1,4,kCyan-6,kOrange+1,7,8,9,10,1,1,1};
0627 
0628   int cl=(id<10)? prt_coll[id]:id;
0629   int cm=(id<10)? prt_colm[id]:id;
0630   g->SetLineColor(cl);
0631   g->SetMarkerColor(cm);
0632   g->SetMarkerStyle(20);
0633   g->SetMarkerSize(0.8);
0634   g->SetName(Form("gr_%d",id));
0635 }
0636 
0637 void prt_setRootPalette(int pal = 0){
0638 
0639   // pal =  1: rainbow\n"
0640   // pal =  2: reverse-rainbow\n"
0641   // pal =  3: amber\n"
0642   // pal =  4: reverse-amber\n"
0643   // pal =  5: blue/white\n"
0644   // pal =  6: white/blue\n"
0645   // pal =  7: red temperature\n"
0646   // pal =  8: reverse-red temperature\n"
0647   // pal =  9: green/white\n"
0648   // pal = 10: white/green\n"
0649   // pal = 11: orange/blue\n"
0650   // pal = 12: blue/orange\n"
0651   // pal = 13: white/black\n"
0652   // pal = 14: black/white\n"
0653 
0654   const int NRGBs = 5;
0655   const int NCont = 255;
0656   gStyle->SetNumberContours(NCont);
0657 
0658   if (pal < 1 && pal> 15) return;
0659   else pal--;
0660 
0661   double stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
0662   double red[15][NRGBs]   = {{ 0.00, 0.00, 0.87, 1.00, 0.51 },
0663                    { 0.51, 1.00, 0.87, 0.00, 0.00 },
0664                    { 0.17, 0.39, 0.62, 0.79, 1.00 },
0665                    { 1.00, 0.79, 0.62, 0.39, 0.17 },
0666                    { 0.00, 0.00, 0.00, 0.38, 1.00 },
0667                    { 1.00, 0.38, 0.00, 0.00, 0.00 },
0668                    { 0.00, 0.50, 0.89, 0.95, 1.00 },
0669                    { 1.00, 0.95, 0.89, 0.50, 0.00 },
0670                    { 0.00, 0.00, 0.38, 0.75, 1.00 },
0671                    { 0.00, 0.34, 0.61, 0.84, 1.00 },
0672                    { 0.75, 1.00, 0.24, 0.00, 0.00 },
0673                    { 0.00, 0.00, 0.24, 1.00, 0.75 },
0674                    { 0.00, 0.34, 0.61, 0.84, 1.00 },
0675                    { 1.00, 0.84, 0.61, 0.34, 0.00 },
0676                    { 0.00, 0.00, 0.80, 1.00, 0.80 }
0677   };
0678   double green[15][NRGBs] = {{ 0.00, 0.81, 1.00, 0.20, 0.00 },          
0679                    { 0.00, 0.20, 1.00, 0.81, 0.00 },
0680                    { 0.01, 0.02, 0.39, 0.68, 1.00 },
0681                    { 1.00, 0.68, 0.39, 0.02, 0.01 },
0682                    { 0.00, 0.00, 0.38, 0.76, 1.00 },
0683                    { 1.00, 0.76, 0.38, 0.00, 0.00 },
0684                    { 0.00, 0.00, 0.27, 0.71, 1.00 },
0685                    { 1.00, 0.71, 0.27, 0.00, 0.00 },
0686                    { 0.00, 0.35, 0.62, 0.85, 1.00 },
0687                    { 1.00, 0.75, 0.38, 0.00, 0.00 },
0688                    { 0.24, 1.00, 0.75, 0.18, 0.00 },
0689                    { 0.00, 0.18, 0.75, 1.00, 0.24 },
0690                    { 0.00, 0.34, 0.61, 0.84, 1.00 },
0691                    { 1.00, 0.84, 0.61, 0.34, 0.00 },
0692                    { 0.00, 0.85, 1.00, 0.30, 0.00 }     
0693   };
0694   double blue[15][NRGBs]  = {{ 0.51, 1.00, 0.12, 0.00, 0.00 },
0695                    { 0.00, 0.00, 0.12, 1.00, 0.51 },
0696                    { 0.00, 0.09, 0.18, 0.09, 0.00 },
0697                    { 0.00, 0.09, 0.18, 0.09, 0.00 },
0698                    { 0.00, 0.47, 0.83, 1.00, 1.00 },
0699                    { 1.00, 1.00, 0.83, 0.47, 0.00 },
0700                    { 0.00, 0.00, 0.00, 0.40, 1.00 },
0701                    { 1.00, 0.40, 0.00, 0.00, 0.00 },
0702                    { 0.00, 0.00, 0.00, 0.47, 1.00 },
0703                    { 1.00, 0.47, 0.00, 0.00, 0.00 },
0704                    { 0.00, 0.62, 1.00, 0.68, 0.12 },
0705                    { 0.12, 0.68, 1.00, 0.62, 0.00 },
0706                    { 0.00, 0.34, 0.61, 0.84, 1.00 },
0707                    { 1.00, 0.84, 0.61, 0.34, 0.00 },
0708                    { 0.60, 1.00, 0.10, 0.00, 0.00 }
0709   };
0710 
0711 
0712   TColor::CreateGradientColorTable(NRGBs, stops, red[pal], green[pal], blue[pal], NCont);
0713  
0714 }
0715 int prt_getColorId(int ind, int style =0){
0716   int cid = 1;
0717   if(style==0) {
0718     cid=ind+1;
0719     if(cid==5) cid =8;
0720     if(cid==3) cid =kOrange+2;
0721   }
0722   if(style==1) cid=ind+300;
0723   return cid;
0724 }
0725 
0726 int prt_shiftHist(TH1 *hist, double double_shift){
0727   int bins=hist->GetXaxis()->GetNbins();
0728   double xmin=hist->GetXaxis()->GetBinLowEdge(1);
0729   double xmax=hist->GetXaxis()->GetBinUpEdge(bins);
0730   double_shift=double_shift*(bins/(xmax-xmin));
0731   int shift=0;
0732   if(double_shift<0) shift=TMath::FloorNint(double_shift);
0733   if(double_shift>0) shift=TMath::CeilNint(double_shift);
0734   if(shift==0) return 0;
0735   if(shift>0){
0736     for(int i=1; i<=bins; i++){
0737       if(i+shift<=bins) hist->SetBinContent(i,hist->GetBinContent(i+shift));
0738       if(i+shift>bins) hist->SetBinContent(i,0);
0739     }
0740     return 0;
0741   }
0742   if(shift<0){
0743     for(int i=bins; i>0; i--){
0744       if(i+shift>0) hist->SetBinContent(i,hist->GetBinContent(i+shift));
0745       if(i+shift<=0) hist->SetBinContent(i,0);
0746     }    
0747     return 0;
0748   }
0749   return 1;
0750 } 
0751 
0752 void prt_addInfo(TString str){
0753   prt_info += str+"\n";
0754 }
0755 
0756 void prt_writeInfo(TString filename){
0757   std::ofstream myfile;
0758   myfile.open (filename);
0759   myfile << prt_info+"\n";
0760   myfile.close();
0761 }
0762 
0763 void prt_writeString(TString filename, TString str){
0764   std::ofstream myfile;
0765   myfile.open (filename);
0766   myfile << str+"\n";
0767   myfile.close();
0768   std::cout<<"output: "<<filename<<std::endl;  
0769 }
0770 
0771 TString prt_createDir(TString inpath=""){
0772   if(inpath != "") prt_savepath = inpath;
0773   TString finalpath = prt_savepath;
0774 
0775   if(finalpath =="") return "";
0776   
0777   if(prt_savepath.EndsWith("auto")) {
0778     TString dir = prt_savepath.ReplaceAll("auto","data");
0779     gSystem->mkdir(dir);
0780     TDatime *time = new TDatime();
0781     TString path(""), stime = Form("%d.%d.%d", time->GetDay(),time->GetMonth(),time->GetYear()); 
0782     gSystem->mkdir(dir+"/"+stime);
0783     for(int i=0; i<1000; i++){
0784       path = stime+"/"+Form("arid-%d",i);
0785       if(gSystem->mkdir(dir+"/"+path)==0) break;
0786     }
0787     gSystem->Unlink(dir+"/last");
0788     gSystem->Symlink(path, dir+"/last");
0789     finalpath = dir+"/"+path;
0790     prt_savepath=finalpath;
0791   }else{
0792     gSystem->mkdir(prt_savepath,kTRUE);
0793   }
0794   
0795   prt_writeInfo(finalpath+"/readme");
0796   return finalpath;
0797 }
0798 
0799 void prt_canvasPrint(TPad *c, TString name="", TString path="", int what=0){
0800   c->Modified();
0801   c->Update();
0802   c->Print(path+"/"+name+".png");
0803   if(what>0) c->Print(path+"/"+name+".C");
0804   if(what>1) c->Print(path+"/"+name+".pdf");
0805   if(what>2) c->Print(path+"/"+name+".eps");    
0806 }
0807 
0808 void prt_set_style(TCanvas *c){
0809   // prt_setRootPalette(1);
0810   if(fabs(c->GetBottomMargin()-0.1)<0.001) c->SetBottomMargin(0.12);
0811   TIter next(c->GetListOfPrimitives());
0812   TObject *obj;
0813     
0814   while((obj = next())){
0815     if(obj->InheritsFrom("TH1")){
0816       TH1F *hh = (TH1F*)obj;
0817       hh->GetXaxis()->SetTitleSize(0.06);
0818       hh->GetYaxis()->SetTitleSize(0.06);
0819 
0820       hh->GetXaxis()->SetLabelSize(0.05);
0821       hh->GetYaxis()->SetLabelSize(0.05);
0822 
0823       hh->GetXaxis()->SetTitleOffset(0.85);
0824       hh->GetYaxis()->SetTitleOffset(0.76);
0825 
0826       if(c->GetWindowHeight()>700){
0827     hh->GetXaxis()->SetTitleSize(0.045);
0828     hh->GetYaxis()->SetTitleSize(0.045);
0829 
0830     hh->GetXaxis()->SetLabelSize(0.035);
0831     hh->GetYaxis()->SetLabelSize(0.035);
0832 
0833     hh->GetXaxis()->SetTitleOffset(0.85);
0834     hh->GetYaxis()->SetTitleOffset(0.98);
0835     c->SetRightMargin(0.12);
0836       }
0837                        
0838       if(c->GetWindowWidth()<=600){
0839         hh->GetXaxis()->SetTitleSize(0.05);
0840         hh->GetYaxis()->SetTitleSize(0.05);
0841 
0842         hh->GetXaxis()->SetLabelSize(0.04);
0843         hh->GetYaxis()->SetLabelSize(0.04);
0844 
0845         hh->GetXaxis()->SetTitleOffset(0.85);
0846         hh->GetYaxis()->SetTitleOffset(0.92);
0847         c->SetRightMargin(0.12);
0848       }
0849 
0850 
0851       if(fabs(c->GetBottomMargin()-0.12)<0.001){
0852     TPaletteAxis *palette = (TPaletteAxis*)hh->GetListOfFunctions()->FindObject("palette");
0853     if(palette) {
0854       palette->SetY1NDC(0.12);
0855       c->Modified();
0856     }
0857       }
0858       c->Modified();
0859       c->Update();
0860     }
0861     
0862     if(obj->InheritsFrom("TGraph")){
0863       TGraph *gg = (TGraph*)obj;
0864       gg->GetXaxis()->SetLabelSize(0.05);
0865       gg->GetXaxis()->SetTitleSize(0.06);
0866       gg->GetXaxis()->SetTitleOffset(0.84);
0867 
0868       gg->GetYaxis()->SetLabelSize(0.05);
0869       gg->GetYaxis()->SetTitleSize(0.06);
0870       gg->GetYaxis()->SetTitleOffset(0.8);
0871     }
0872  
0873     if(obj->InheritsFrom("TMultiGraph")){
0874       TMultiGraph *gg = (TMultiGraph*)obj;
0875       gg->GetXaxis()->SetLabelSize(0.05);
0876       gg->GetXaxis()->SetTitleSize(0.06);
0877       gg->GetXaxis()->SetTitleOffset(0.84);
0878 
0879       gg->GetYaxis()->SetLabelSize(0.05);
0880       gg->GetYaxis()->SetTitleSize(0.06);
0881       gg->GetYaxis()->SetTitleOffset(0.8);
0882     }
0883     
0884     if(obj->InheritsFrom("TF1")){
0885       TF1 *f = (TF1*)obj;
0886       f->SetNpx(500);
0887     }
0888   }
0889 }
0890 
0891 void prt_save(TPad *c= NULL,TString path="", int what=0, int style=0){
0892   TString name = c->GetName();
0893   Bool_t batch = gROOT->IsBatch();
0894   gROOT->SetBatch(1);
0895   
0896   if(c && path != "") {
0897     int w = 800, h = 400;
0898     if(style != -1){
0899       if(style == 1) {w = 800; h = 500;}
0900       if(style == 2) {w = 800; h = 600;}
0901       if(style == 3) {w = 800; h = 400;}
0902       if(style == 5) {w = 800; h = 900;} 
0903       if(style == 0){ 
0904         w = ((TCanvas*)c)->GetWindowWidth();
0905         h = ((TCanvas*)c)->GetWindowHeight();
0906       }
0907       
0908       TCanvas *cc;
0909       if(TString(c->GetName()).Contains("hp") || TString(c->GetName()).Contains("cdigi")) {
0910     cc = prt_drawDigi(prt_last_layoutId,prt_last_maxz,prt_last_minz);
0911     cc->SetCanvasSize(800,400);
0912     if(name.Contains("=")) name =  name.Tokenize('=')->First()->GetName();  
0913       }else{
0914         cc = new TCanvas(TString(c->GetName())+"exp","cExport",0,0,w,h);
0915         cc = (TCanvas*) c->DrawClone();
0916     cc->SetCanvasSize(w,h);
0917       }
0918       
0919       if(style == 0) prt_set_style(cc);
0920       
0921       prt_canvasPrint(cc,name,path,what);
0922     }else{
0923       c->SetCanvasSize(w,h);
0924       prt_canvasPrint(c,name,path,what);
0925     }
0926   }
0927   
0928   gROOT->SetBatch(batch);
0929 }
0930 
0931 TString prt_createSubDir(TString dir="dir"){
0932   gSystem->mkdir(dir);
0933   return dir;
0934 }
0935 
0936 TList *prt_canvaslist;
0937 
0938 TCanvas *prt_canvasGet(TString name="c"){
0939   TIter next(prt_canvaslist);
0940   TCanvas *c = 0;
0941   while((c = (TCanvas*) next())){
0942     if(c->GetName()==name || name=="*") break;
0943   }
0944   return c;
0945 }
0946 
0947 void prt_canvasAdd(TString name="c",int w=800, int h=400){
0948   if(!prt_canvasGet(name)){
0949     if(!prt_canvaslist) prt_canvaslist = new TList();
0950     TCanvas *c = new TCanvas(name,name,0,0,w,h); 
0951     prt_canvaslist->Add(c);
0952   }
0953 }
0954 
0955 void prt_canvasAdd(TCanvas *c){
0956   if(!prt_canvaslist) prt_canvaslist = new TList();
0957   c->cd();
0958   prt_canvaslist->Add(c);
0959 }
0960 
0961 void prt_canvasCd(){
0962   
0963 }
0964 
0965 void prt_canvasDel(TString name="c"){
0966   TIter next(prt_canvaslist);
0967   TCanvas *c=0;
0968   while((c = (TCanvas*) next())){
0969     if(c->GetName()==name || name=="*") prt_canvaslist->Remove(c);
0970     c->Delete();
0971   }
0972 }
0973 
0974 // style = 0 - for web blog
0975 // style = 1 - for talk 
0976 // what = 0 - save in png
0977 // what = 1 - save in png, C
0978 // what = 2 - save in png, C, pdf
0979 // what = 3 - save in png, C, pdf, eps 
0980 void prt_canvasSave(int what=1, int style=0, Bool_t rm=false){
0981   TIter next(prt_canvaslist);
0982   TCanvas *c=0;
0983   TString path = prt_createDir();
0984   while((c = (TCanvas*) next())){
0985     prt_set_style(c);
0986     prt_save(c, path, what,style);
0987     if(rm){
0988       prt_canvaslist->Remove(c);
0989       c->Close();
0990     }
0991   }
0992 }
0993 
0994 // path - folder for saving
0995 void prt_canvasSave(TString path, int what=1, int style=0, Bool_t rm=false){
0996   prt_savepath = path;
0997   prt_canvasSave(what,style,rm);
0998 }
0999 
1000 void prt_set_style(){
1001   TIter next(prt_canvaslist);
1002   TCanvas *c=0;
1003   TString path = prt_createDir();
1004   while((c = (TCanvas*) next())){
1005     prt_set_style(c);
1006   }
1007 }
1008 
1009 void prt_waitPrimitive(TString name, TString prim=""){
1010   TIter next(prt_canvaslist);
1011   TCanvas *c=0;
1012   while((c = (TCanvas*) next())){
1013     if(TString(c->GetName())==name){
1014       c->Modified(); 
1015       c->Update(); 
1016       c->WaitPrimitive(prim);
1017     }
1018   }
1019 }
1020 
1021 double prt_integral(TH1F *h,double xmin, double xmax){
1022   TAxis *axis = h->GetXaxis();
1023   int bmin = axis->FindBin(xmin);
1024   int bmax = axis->FindBin(xmax);
1025   double integral = h->Integral(bmin,bmax);
1026   integral -= h->GetBinContent(bmin)*(xmin-axis->GetBinLowEdge(bmin))/axis->GetBinWidth(bmin);
1027   integral -= h->GetBinContent(bmax)*(axis->GetBinUpEdge(bmax)-xmax)/axis->GetBinWidth(bmax);
1028   return integral;
1029 }
1030 
1031 void prt_normalize(TH1F* hists[],int size){
1032   // for(int i=0; i<size; i++){
1033   //   hists[i]->Scale(1/hists[i]->Integral(), "width"); 
1034   // }
1035   
1036   double max = 0;
1037   double min = 0;
1038   for(int i=0; i<size; i++){
1039     double tmax =  hists[i]->GetBinContent(hists[i]->GetMaximumBin());
1040     double tmin = hists[i]->GetMinimum();
1041     if(tmax>max) max = tmax;
1042     if(tmin<min) min = tmin;
1043   }
1044   max += 0.05*max;
1045   for(int i=0; i<size; i++){
1046     hists[i]->GetYaxis()->SetRangeUser(min,max);
1047   }
1048 }
1049 
1050 void prt_normalizeto(TH1F* hists[],int size, double max=1){
1051   for(int i=0; i<size; i++){
1052     double tmax =  hists[i]->GetBinContent(hists[i]->GetMaximumBin());
1053     if(tmax>0)hists[i]->Scale(max/tmax);
1054   }
1055 }
1056 
1057 void prt_normalize(TH1F* h1,TH1F* h2){
1058   double max = (h1->GetMaximum()>h2->GetMaximum())? h1->GetMaximum() : h2->GetMaximum();
1059   max += max*0.1;
1060   h1->GetYaxis()->SetRangeUser(0,max);
1061   h2->GetYaxis()->SetRangeUser(0,max);
1062 }
1063 
1064 // just x for now
1065 TGraph* prt_smooth(TGraph* g,int smoothness=1){
1066   double x, y;
1067   int n = g->GetN();
1068   TH1F *h = new TH1F("h","h",g->GetN(),0,n);
1069   TGraph *gr = new TGraph();
1070   gr->SetName(g->GetName());
1071   for(auto i=0; i<n; i++){
1072     g->GetPoint(i,x,y);
1073     h->Fill(i,x);
1074   }
1075 
1076   h->Smooth(smoothness);
1077   
1078   for(auto i=0;i<n;i++){
1079     g->GetPoint(i,x,y);
1080     gr->SetPoint(i,h->GetBinContent(i),y);
1081   }
1082   return gr;
1083 }
1084 
1085 int prt_get_pid(int pdg){
1086   int pid=0;
1087   if(pdg==-11)   pid=0; //e
1088   if(pdg==13)   pid=1; //mu
1089   if(pdg==211)  pid=2; //pi
1090   if(pdg==321)  pid=3; //K
1091   if(pdg==2212) pid=4; //p
1092   return pid;
1093 }
1094 
1095 double prt_get_momentum_from_tof(double dist,double dtof){
1096   double s = dtof*0.299792458/dist;
1097   double x = s*s;
1098   double a = prt_mass[2]*prt_mass[2]; //pi
1099   double b = prt_mass[4]*prt_mass[4]; //p
1100   double p = sqrt((a - 2*sqrt(a*a+a*b*x-2*a*b+b*b)/s + b)/(x - 4));
1101 
1102   return p;
1103 }
1104 
1105 // return TOF difference [ns] for mominum p [GeV] and flight path l [m]
1106 double prt_get_tof_diff(int pid1=211, int pid2=321, double p=1, double l=2){
1107   double c = 299792458;
1108   double m1 = prt_mass[prt_get_pid(pid1)];
1109   double m2 = prt_mass[prt_get_pid(pid2)];
1110   double td = l*(sqrt(p*p+m1*m1)-sqrt(p*p+m2*m2))/(p*c)*1E9;
1111 
1112   // relativistic
1113   // td = l*(m1*m1-m2*m2)/(2*p*p*c)*1E9;
1114   
1115   return td;
1116 }
1117 
1118 bool prt_ispath(TString path){
1119   Long_t *id(0),*size(0),*flags(0),*modtime(0);
1120   return !gSystem->GetPathInfo(path,id,size,flags,modtime);
1121 }
1122 
1123 int prt_get3digit(TString str){
1124   TPRegexp e("[0-9][0-9][0-9]");
1125   return ((TString)str(e)).Atoi();
1126 }
1127 
1128 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!