StRoot: StHighptPool/macros/commonmacro/histutil.h Source File

00001 //#ifndef HISTUTIL_H
00002 //#define HISTUTIL_H
00003 
00004 #include <cstring>
00005 
00006 // 3d to 2d
00007 TH2D* HistSlice(TH3* h3,const char* basename, const char* basetitle,
00008                const char* slicename, Axis_t min, Axis_t max,
00009                const char* slicetype,
00010                Option_t* opt="");
00011 
00012 
00013 // 3d to 2d but keep varying width histos
00014 TH2D* HistSlice(TH3* h3,const char* basename, const char* basetitle,
00015                const char* slicename, const char* slicetype);
00016 
00017 // 2d to 1d
00018 TH1D* HistSlice(TH2* h2,const char* basename, const char* basetitle,
00019                const char* slicename,Axis_t min, Axis_t max,
00020                const char* projaxis,
00021                Option_t* opt="");
00022 
00023 // 3d to 1d
00024 // if project to x then slicename1="y",slicename2="z"
00025 // if project to y then slicename2="z",slicename2="x"
00026 // if project to z then slicename1="x",slicename2="y"
00027 
00028 TH1D* HistSlice(TH3* h3,const char* basename, const char* basetitle,
00029                const char* slicename1,Axis_t min1, Axis_t max1,
00030                const char* slicename2,Axis_t min2, Axis_t max2,
00031                const char* projecttype,
00032                Option_t* opt="");
00033 
00034 TProfile* Profile(TH3* h3,const char* basename, const char* basetitle,
00035                   const char* slicename, Axis_t min, Axis_t max,
00036                   const char* slicetype,
00037                   const char* projecttype,
00038                   Option_t* opt="");
00039 
00040 TProfile* Profile(TH2* h2,const char* basename,const char* basetitle,
00041                   const char* slicename, Axis_t min, Axis_t max,
00042                   const char* profaxis,
00043                   Option_t* opt="");
00044 
00045 TH1D* Rms(TH3* h3,const char* basename, const char* basetitle,
00046                   const char* slicename, Axis_t min, Axis_t max,
00047                   const char* slicetype,
00048                   const char* projecttype);
00049 TH1D* Rms(TH2* h2,const char* basename,
00050               const char* projecttype);
00051 // rebin an axis of a 3D histogram (with varying bins);
00052 
00053 void Rebin(TH3* h3,int axisType,int nBin,Axis_t* ary);
00054 
00055 Int_t FindMinBin(const TAxis* axis,Stat_t val);
00056 Int_t FindMaxBin(const TAxis* axis,Stat_t val);
00057 void  SetRange(const TAxis* axis,Stat_t min,Stat_t max);
00058 Stat_t FindMax(TH1* ha,TH1* hb);
00059 Stat_t FindMin(TH1* ha,TH1* hb);
00060 Stat_t FindMax(TH1* ha,TH1* hb,TH1* hc);
00061 void Divide(TCanvas* c1,int nx,int ny,const char* title=0,const char* file=0);
00062 void Print(TCanvas* c1, const char* outDir, TString basename);
00063 
00064 void SetMinMax(TH1* h);
00065 void SetMinMax(TH1* h,float min, float max);
00066 void GetMeanRms(TH1* ha,float& cmean,float& crms);
00067 void GetMeanRmsYield(TH1* ha,float& cmean,float& crms);
00068 void PrintMeanRms(TH1* ha,float x, float y,float size=0.04);
00069 void PrintMeanRmsYield(TH1* ha,float x, float y,float size=0.04);
00070 void ReplaceTitle(TH1* h,char* current, char* replace);
00071 void DrawLine(TH1* h,float y=0);
00072 
00073 void showHistogramValues(TH1D* h1, char* more=0);
00074 void showTGraphValues(TGraphAsymmErrors* graph,char* more=0);
00075 void dump(TH1* h1);
00076 void dump(TGraph* g);
00077 void dump(TGraphAsymmErrors* g);
00078 void scale2D(TH2*);
00079 //-----------------------------------------------
00080 
00081 void DrawLine(TH1* h,float y)
00082 {
00083   TLine* line=new TLine; TAxis* axis=h->GetXaxis();
00084   line->DrawLine(axis->GetBinLowEdge(axis->GetFirst()),y,
00085                  axis->GetBinUpEdge(axis->GetLast()),y);
00086 
00087 }
00088 
00089 void ReplaceTitle(TH1* h,char* current, char* replace)
00090 {
00091   TString s=h->GetTitle();
00092   s.ReplaceAll(current,replace);
00093   h->SetTitle(s.Data());
00094 }
00095 
00096 void SetMinMax(TH1* h,float min, float max)
00097 {
00098   h->SetMinimum(min);
00099   h->SetMaximum(max);
00100 }
00101 
00102 void SetMinMax(TH1* h)
00103 {
00104   h->SetMinimum(h->GetMinimum()*0.8);
00105   h->SetMaximum(h->GetMaximum()*1.2);
00106 }
00107 
00108 void GetMeanRms(TH1* ha,float& mean,float& rms)
00109 {
00110   TAxis* axis=ha->GetXaxis();
00111   mean=0;rms=0;
00112   int n=0;
00113 
00114   for(int i=axis->GetFirst(); i<=axis->GetLast(); i++){
00115     if(ha->GetBinContent(i)>0){
00116       float val = (ha->GetBinCenter(i)*ha->GetBinContent(i));
00117       mean += val;
00118       rms += val*val;
00119 
00120     }
00121   }
00122   n=ha->Integral(axis->GetFirst(),axis->GetLast());
00123   //n=ha->Integral(1,axis->GetNbins());
00124   if(n>0){
00125     mean /= n; rms /= n; rms=sqrt(rms - mean*mean);
00126   }
00127 }
00128 
00129 void GetMeanRmsYield(TH1* ha,float& mean,float& rms)
00130 {
00131   TAxis* axis=ha->GetXaxis();
00132   mean=0;rms=0;
00133   int n=0;
00134 
00135   //  cout << ha->GetName() << endl;
00136 
00137   for(int i=axis->GetFirst(); i<=axis->GetLast(); i++){
00138     if(ha->GetBinContent(i)>0){
00139       float val = ha->GetBinContent(i);
00140       mean += val;
00141       rms += val*val;
00142       //  cout << "bin=" << i << ",val=" << val << endl;
00143       
00144     }
00145   }
00146   n=axis->GetLast()-axis->GetFirst()+1;
00147   //  cout << "n=" << n << endl;
00148   //n=ha->Integral(1,axis->GetNbins());
00149   if(n>0){
00150     mean /= n; rms /= n; rms=sqrt(rms - mean*mean);
00151   }
00152   //  cout << mean << ", " << rms << endl;
00153 }
00154 
00155 void PrintMeanRms(TH1* ha,float x,float y,float size)
00156 {
00157   TText* text=new TText; float mean(0),rms(0);
00158   text->SetTextSize(size);
00159   char buf[100];
00160   mean = ha->GetMean(); rms = ha->GetRMS();
00161   //  if(compute)GetMeanRms(ha,mean,rms); 
00162   sprintf(buf,"mean: %.3f",mean);
00163   text->DrawTextNDC(x,y,buf);
00164   sprintf(buf,"rms: %.3f",rms);
00165   text->DrawTextNDC(x,y-0.05,buf);
00166 }
00167 
00168 void PrintMeanRmsYield(TH1* ha,float x,float y,float size)
00169 {
00170   TText* text=new TText; float mean(0),rms(0);
00171   text->SetTextSize(size);
00172   char buf[100];
00173   //mean = ha->GetMean(); rms = ha->GetRMS();
00174   GetMeanRmsYield(ha,mean,rms); 
00175   sprintf(buf,"mean: %.3f",mean);
00176   text->DrawTextNDC(x,y,buf);
00177   sprintf(buf,"rms: %.3f",rms);
00178   text->DrawTextNDC(x,y-0.05,buf);
00179 }
00180 
00181 // divide a canvas and print a title
00182 
00183 void Divide(TCanvas *c1,int nx,int ny,const char* title,const char* file)
00184 {
00185   c1->Clear();  
00186   TString cTitle = (title) ? title : c1->GetTitle();
00187   if(title) c1->SetTitle(title);
00188   
00189   TPaveLabel *mP = new TPaveLabel(.02,.95,.5,.98,title);
00190   mP->SetFillColor(0);
00191   mP->SetBorderSize(0);
00192   mP->Draw();
00193   
00194   if(file){
00195     TText tFile;
00196     tFile.SetTextSize(0.02);
00197     tFile.DrawTextNDC(0.1,0.0,file);
00198   }
00199 
00200   char buf[10];
00201   
00202   Double_t ymax = .94;
00203   Double_t dy = ymax/ny;
00204   Double_t dx = 1./nx;
00205   Double_t ymargin = 0.02;
00206   double   xmargin = 0.02;
00207   int n=0;
00208   for (Int_t iy=0;iy<ny;iy++) {
00209     double yhigh = ymax - iy*dy - ymargin;
00210     double ylow = yhigh - dy + 3*ymargin;
00211     if (ylow < 0) ylow = 0;
00212     if (ylow > yhigh) continue;
00213     for (Int_t ix=0;ix<nx;ix++) {
00214       n++;
00215       double xlow = ix*dx + xmargin;
00216       double xhigh = xlow +dx -2*xmargin;
00217       if (xlow > xhigh) continue;
00218       
00219       TPad* p = new TPad(buf,buf,xlow,ylow,xhigh,yhigh);
00220       p->SetNumber(n);
00221       p->Draw();
00222     }
00223   }
00224   
00225 }
00226                 
00227 // given a maximum value, find it's bin
00228 Int_t 
00229 FindMinBin(const TAxis* axis,Stat_t val)
00230 {
00231   Int_t bin = axis->FindBin(val);
00232   if(bin<0) return 1;
00233   return (axis->GetBinCenter(bin)<val) ? ++bin : bin;
00234   
00235 }
00236  
00237 // given a minimum value, find it's bin
00238 Int_t FindMaxBin(const TAxis* axis,Stat_t val)
00239 {
00240   Int_t bin = axis->FindBin(val);
00241   if(bin>axis->GetNbins()) return axis->GetNbins();
00242   return (axis->GetBinCenter(bin)>val) ? --bin : bin;
00243 }
00244 
00245 void SetRange(const TAxis* axis, Stat_t min, Stat_t max)
00246 {
00247   axis->SetRange(FindMinBin(axis,min),FindMaxBin(axis,max));
00248 }
00249 Stat_t
00250 FindMax(TH1* ha, TH1* hb)
00251 {
00252   Stat_t maxa = ha->GetMaximum();
00253   Stat_t maxb = hb->GetMaximum();
00254   return (maxa > maxb) ? maxa : maxb;
00255 
00256 }
00257 
00258 Stat_t
00259 FindMin(TH1* ha, TH1* hb)
00260 {
00261   Stat_t mina = ha->GetMinimum();
00262   Stat_t minb = hb->GetMinimum();
00263   return (mina > minb) ? mina : minb;
00264 
00265 }
00266 
00267 Stat_t 
00268 FindMax(TH1* ha,TH1* hb,TH1* hc)
00269 {
00270   Stat_t max1 = FindMax(ha,hb);
00271   Stat_t max2 = FindMax(hb,hc);
00272   return (max1>max2) ? max1 : max2;
00273 }
00274 
00275 
00276 // 3d to 2d
00277 
00278 TH2D* HistSlice(TH3* h3,
00279                const char* basename,const char* basetitle,
00280                const char* slicename, Axis_t min, Axis_t max,
00281                const char* slicetype,
00282                Option_t* opt)
00283 {
00284   char name[200], title[200],option[20],buf[100];
00285   Int_t minBin, maxBin;
00286   Axis_t mmin,mmax;
00287   TH2D* h2;
00288 
00289   if(!h3){
00290     cout << "Null pointer for h3 : " << basename << endl;
00291     return 0;
00292   }
00293 
00294   // set the range
00295   // x slice
00296 
00297   TAxis* axis = 0; TAxis* yAxis=0,*xAxis=0;
00298   if(strcmp(slicetype,"yz")==0 || strcmp(slicetype,"zy")==0 ){
00299     axis = h3->GetXaxis();
00300     if(strcmp(slicetype,"yz")==0){ yAxis=h3->GetYaxis(); xAxis=h3->GetZaxis();}
00301     else { yAxis=h3->GetZaxis(); xAxis=h3->GetYaxis(); }
00302   }
00303   // y slice
00304   else if (strcmp(slicetype,"xz")==0 || strcmp(slicetype,"zx")==0){
00305     axis = h3->GetYaxis();
00306     if(strcmp(slicetype,"xz")==0){ yAxis=h3->GetXaxis(); xAxis=h3->GetZaxis();}
00307     else { yAxis=h3->GetZaxis(); xAxis=h3->GetXaxis(); }
00308   }
00309   // z slice
00310   else if (strcmp(slicetype,"xy")==0 || strcmp(slicetype,"yx")==0){
00311     axis = h3->GetZaxis();
00312     if(strcmp(slicetype,"xy")==0){ yAxis=h3->GetXaxis(); xAxis=h3->GetYaxis();}
00313     else { yAxis=h3->GetYaxis(); xAxis=h3->GetXaxis(); }
00314   }
00315   else {
00316     cerr << "Wrong value for slice axis for " << basename << endl;
00317     return 0;
00318   }    
00319   
00320   if(min>=max){
00321     minBin = 1; maxBin = axis->GetNbins();
00322   }
00323   else{
00324     minBin = FindMinBin(axis,min);
00325     maxBin = FindMaxBin(axis,max);
00326   }
00327   sprintf(option,"%s%s",slicetype,opt);
00328 
00329   if(!slicename){
00330     strcpy(buf,axis->GetTitle()); slicename=buf;
00331   }
00332 
00333 
00334   axis->SetRange(minBin,maxBin);
00335 
00336   h2 = (TH2D*) h3->Project3D(option);
00337   h2->SetYTitle(yAxis->GetTitle());
00338   h2->SetXTitle(xAxis->GetTitle());
00339 
00340 
00341   mmin = axis->GetBinLowEdge(minBin);
00342   mmax = axis->GetBinUpEdge(maxBin);
00343   
00344   // set the name and title
00345   sprintf(name,"%s(%.2f<%s<%.2f)",basename,mmin,slicename,mmax);
00346   sprintf(title,"%s (%.2f<%s<%.2f)",basetitle,mmin,slicename,mmax);
00347   
00348   h2->SetName(name); h2->SetTitle(title);
00349 
00350   axis->SetRange(0,9999999);
00351 
00352   return h2;
00353 }
00354 
00355 // 3d to 2d var bins
00356 
00357 TH2D* HistSlice(TH3* h3,const char* basename, const char* basetitle,
00358                const char* slicename, const char* slicetype)
00359 {
00360   TString opt = slicetype;
00361 
00362   int projX=-1,projY=-1,projZ=-1;
00363 
00364   if(opt.Contains("xy")){ projX=1; projY=0; projZ=2; }
00365   else if(opt.Contains("yx")){ projX=0; projY=1; projZ=2;}
00366   else if(opt.Contains("zy")){ projX=1; projY=2; projZ=0; }
00367   else if(opt.Contains("yz")){ projX=2; projY=1; projZ=0; }
00368   else if(opt.Contains("xz")){ projX=2; projY=0; projZ=1; }
00369   else if(opt.Contains("zx")){ projX=0; projY=2; projZ=1; }
00370 
00371   TAxis* xAxis=0, *yAxis=0, *zAxis=0;
00372   switch(projX){
00373   case 0: xAxis = h3->GetXaxis(); break;
00374   case 1: xAxis = h3->GetYaxis(); break;
00375   case 2: xAxis = h3->GetZaxis(); break;
00376   }
00377   switch(projY){
00378   case 0: yAxis = h3->GetXaxis(); break;
00379   case 1: yAxis = h3->GetYaxis(); break;
00380   case 2: yAxis = h3->GetZaxis(); break;
00381   }
00382   switch(projZ){
00383   case 0: zAxis = h3->GetXaxis(); break;
00384   case 1: zAxis = h3->GetYaxis(); break;
00385   case 2: zAxis = h3->GetZaxis(); break;
00386   }
00387   
00388   // create the 2d histogram
00389   int ixmin = xAxis->GetFirst();
00390   int ixmax = xAxis->GetLast();
00391   int iymin = yAxis->GetFirst();
00392   int iymax = yAxis->GetLast();
00393   int izmin = zAxis->GetFirst();
00394   int izmax = zAxis->GetLast();
00395   int nx = ixmax-ixmin+1;
00396   int ny = iymax-iymin+1;
00397   int nz = izmax-izmin+1;
00398 
00399   char* name = new char[strlen(basename)+200];
00400   char* title = new char[strlen(basename)+200];
00401   
00402   sprintf(name,"%sSbin%dto%d",basename,izmin,izmax);
00403   sprintf(title,"%s(%.2f<%s<%.2f)",basetitle,
00404           zAxis->GetBinLowEdge(izmin),slicename,zAxis->GetBinUpEdge(izmax));
00405 
00406   /*
00407   cout << "xmin=" << xAxis->GetBinLowEdge(ixmin) 
00408        << ",xmax=" << xAxis->GetBinUpEdge(ixmax)
00409        << endl
00410        << "ymin=" << yAxis->GetBinLowEdge(iymin)
00411        << ",ymax=" << yAxis->GetBinUpEdge(iymax)
00412        << endl
00413        << "zmin=" << zAxis->GetBinLowEdge(izmin)
00414        << ",zmax=" << zAxis->GetBinUpEdge(izmax)
00415        << endl;
00416   */
00417 
00418   double* xAry=0; xAry=xAxis->GetXbins()->GetArray();
00419   double* yAry=0; yAry=yAxis->GetXbins()->GetArray();
00420 
00421   // normal binning
00422   TH2D* h2=0;
00423   if(!xAry){
00424     h2=(TH2D*) h3->Project3D(slicetype);
00425     h2->SetName(name); h2->SetTitle(title);
00426     return h2;
00427   }
00428 
00429   // varying binning
00430   h2=new TH2D(name,title,
00431               nx,&xAry[ixmin-1],ny,&yAry[iymin-1]);
00432 
00433   // fill the new histogram
00434 
00435   int iBinXreal=0, iBinYreal=0, iBinZreal=0;
00436   for(int ix=ixmin; ix<=ixmax; ix++){
00437     switch(projX){
00438     case 0: iBinXreal = ix; break;
00439     case 1: iBinYreal = ix; break;
00440     case 2: iBinZreal = ix; break;
00441     }
00442     for(int iy=iymin; iy<=iymax; iy++){
00443       switch(projY){
00444       case 0: iBinXreal = iy; break;
00445       case 1: iBinYreal = iy; break;
00446       case 2: iBinZreal = iy; break;
00447       }
00448       double content=0;
00449       for(int iz=izmin; iz<=izmax; iz++){
00450         switch(projZ){
00451         case 0: iBinXreal = iz; break;
00452         case 1: iBinYreal = iz; break;
00453         case 2: iBinZreal = iz; break;
00454         }
00455         int bin = h3->GetBin(iBinXreal,iBinYreal,iBinZreal);
00456         h2->Fill(xAxis->GetBinCenter(ix),
00457                  yAxis->GetBinCenter(iy),h3->GetBinContent(bin));
00458       }
00459     }
00460   }
00461 
00462   return h2;
00463 }
00464 
00465 
00466 // 2d to 1d
00467 
00468 TH1D* HistSlice(TH2* h2,const char* basename, const char* basetitle,
00469                const char* slicename,Axis_t min, Axis_t max,
00470                const char* projaxis,
00471                Option_t* opt)
00472 {
00473   char name[200], title[200];
00474   Axis_t mmin(0),mmax(0);
00475   Int_t minBin(0), maxBin(0);
00476   TH1D* h1;
00477   TAxis* axis = 0;
00478   char buf[100];
00479 
00480   if(strcmp(projaxis,"x")==0){
00481     axis = h2->GetYaxis();
00482   }
00483   else if(strcmp(projaxis,"y")==0){
00484     axis = h2->GetXaxis();
00485   }
00486   else{
00487     cout << "The argument must be either 'x' or 'y' " << endl;
00488     return 0;
00489   }
00490   
00491   if(min>=max){
00492     minBin = 1; maxBin = axis->GetNbins();
00493  
00494   }
00495   else{
00496     minBin = FindMinBin(axis,min);
00497     maxBin = FindMaxBin(axis,max);
00498   }
00499 
00500   if(strcmp(projaxis,"x")==0){
00501     h1= h2->ProjectionX("dummy",minBin,maxBin,opt);
00502     if(!slicename){
00503       strcpy(buf,axis->GetTitle()); slicename=buf;
00504     }    
00505     h1->SetXTitle(h2->GetXaxis()->GetTitle());
00506   }
00507   else if(strcmp(projaxis,"y")==0){
00508     h1= h2->ProjectionY("dummy",minBin,maxBin,opt);
00509     if(!slicename){
00510       strcpy(buf,axis->GetTitle()); slicename=buf;
00511     }    
00512     h1->SetXTitle(h2->GetYaxis()->GetTitle());
00513   }
00514 
00515   mmin = axis->GetBinLowEdge(minBin);
00516   mmax = axis->GetBinUpEdge(maxBin);
00517   
00518 
00519   //
00520   // change the name and title
00521   //
00522   sprintf(name,"%s(%.2f<%s<%.2f)",basename,mmin,slicename,mmax);
00523   sprintf(title,"%s (%.2f<%s<%.2f)",basetitle,mmin,slicename,mmax);
00524   h1->SetTitle(title);
00525   h1->SetName(name);
00526   
00527   return h1;
00528 }
00529 
00530 TH1D* HistSlice(TH3* h3,const char* basename, const char* basetitle,
00531                const char* slicename1,Axis_t min1, Axis_t max1,
00532                const char* slicename2,Axis_t min2, Axis_t max2,
00533                const char* projecttype,
00534                Option_t* opt)
00535 {
00536   char slicetype[10];
00537   char buf1[100],buf2[100];
00538   TAxis* axis1, *axis2;
00539   if(strcmp(projecttype,"x")==0){ // slice out y
00540     strcpy(slicetype,"zx");
00541     axis1=h3->GetYaxis(); axis2=h3->GetZaxis();
00542   }
00543   else if(strcmp(projecttype,"y")==0){ // slice out z
00544     strcpy(slicetype,"xy");
00545     axis1=h3->GetZaxis(); axis2=h3->GetXaxis();
00546   }
00547   else if(strcmp(projecttype,"z")==0){ // slice out x
00548     strcpy(slicetype,"yz");
00549     axis1=h3->GetXaxis(); axis2=h3->GetYaxis();
00550   }
00551   else{
00552     cout << "error: wrong project type " << basename << endl;
00553     return;
00554   }
00555   if(!slicename1){
00556     strcpy(buf1,axis1->GetTitle()); slicename1=buf1;
00557   }
00558   if(!slicename2){
00559     strcpy(buf2,axis2->GetTitle()); slicename2=buf2;
00560   }
00561 
00562   TH2D* h2 =
00563     (TH2D*) HistSlice(h3,basename,basetitle,
00564                       slicename1,min1,max1,
00565                       slicetype,opt);
00566 
00567   return HistSlice(h2,h2->GetName(),h2->GetTitle(),
00568                    slicename2,min2,max2,
00569                    "x",opt);
00570                       
00571 }
00572 
00573 TProfile* Profile(TH3* h3,const char* basename, const char* basetitle,
00574                   const char* slicename, Axis_t min, Axis_t max,
00575                   const char* slicetype,
00576                   const char* projecttype,
00577                   Option_t* opt)
00578 {
00579   // first project to 2d
00580   TH2D* h2 = HistSlice(h3,basename,basetitle,slicename,min,max,
00581                       slicetype);
00582   char name[200];
00583   TProfile* p;
00584   // profile
00585   if(strstr(projecttype,"x")){
00586     sprintf(name,"%s_profx",h2->GetName());
00587     p=h2->ProfileX(name,0,9999999,opt);
00588   }
00589   else{
00590     sprintf(name,"%s_profy",h2->GetName());
00591     p=h2->ProfileY(name,0,9999999,opt);
00592   }
00593   p->SetTitle(h2->GetTitle());
00594   return p;
00595 }
00596 
00597 TProfile* Profile(TH2* h2,const char* basename,const char* basetitle,
00598                   const char* slicename, Axis_t min, Axis_t max,
00599                   const char* profaxis,
00600                   Option_t* opt)
00601 {
00602   TProfile* p;
00603   char name[200], title[200];
00604   Axis_t mmin(0),mmax(0);
00605   Int_t minBin(0), maxBin(0);
00606   TAxis* axis = 0;
00607   char buf[100];
00608 
00609   if(strcmp(profaxis,"x")==0){
00610     axis = h2->GetYaxis();
00611   }
00612   else if(strcmp(profaxis,"y")==0){
00613     axis = h2->GetXaxis();
00614   }
00615   else{
00616     cout << "The argument must be either 'x' or 'y' " << endl;
00617     return 0;
00618   }
00619   
00620   if(min>=max){
00621     minBin = 1; maxBin = axis->GetNbins();
00622  
00623   }
00624   else{
00625     minBin = FindMinBin(axis,min);
00626     maxBin = FindMaxBin(axis,max);
00627   }
00628 
00629   if(strcmp(profaxis,"x")==0){
00630     p= h2->ProfileX("dummy",minBin,maxBin,opt);
00631     if(!slicename){
00632       strcpy(buf,axis->GetTitle()); slicename=buf;
00633     }    
00634     p->SetXTitle(h2->GetXaxis()->GetTitle());
00635   }
00636   else if(strcmp(profaxis,"y")==0){
00637     p= h2->ProfileY("dummy",minBin,maxBin,opt);
00638     if(!slicename){
00639       strcpy(buf,axis->GetTitle()); slicename=buf;
00640     }    
00641     p->SetXTitle(h2->GetYaxis()->GetTitle());
00642   }
00643 
00644   mmin = axis->GetBinLowEdge(minBin);
00645   mmax = axis->GetBinUpEdge(maxBin);
00646   
00647 
00648   //
00649   // change the name and title
00650   //
00651   sprintf(name,"%s(%.2f<%s<%.2f)",basename,mmin,slicename,mmax);
00652   sprintf(title,"%s (%.2f<%s<%.2f)",basetitle,mmin,slicename,mmax);
00653   p->SetTitle(title);
00654   p->SetName(name);
00655   
00656   return p;
00657 
00658 }
00659 
00660 TProfile* Rms(TH2* h2,const char* basename,
00661               const char* projecttype)
00662 {
00663   char name[100], xtitle[100];
00664   TProfile* p;
00665   if(strstr(projecttype,"x")){
00666     p=h2->ProfileX("dummy",1,h2->GetNbinsY(),"s");
00667     strcpy(xtitle,h2->GetXaxis()->GetTitle());
00668   }
00669   else if(strstr(projecttype,"y")){
00670     p=h2->ProfileY("dummy",1,h2->GetNbinsX(),"s");
00671     strcpy(xtitle,h2->GetYaxis()->GetTitle());
00672   }
00673   else{
00674     cout << "wrong type " << projecttype << endl; exit(1);
00675   }
00676   sprintf(name,"%s%s_%s",h2->GetName(),basename,projecttype);
00677   p->SetName(name);
00678   sprintf(name,"%s_%s_rms",h2->GetName(),basename,projecttype);
00679   Int_t nBin= p->GetNbinsX();
00680   TH1* h1 = new TH1D(name,name,nBin,
00681                      p->GetXaxis()->GetBinLowEdge(1),
00682                      p->GetXaxis()->GetBinUpEdge(nBin));
00683   // now loop over the tprofile and fill the histogram
00684   for(int i=1;i<=nBin;i++){
00685     h1->SetBinContent(i,p->GetBinError(i));
00686     // set the error as 1% of the content
00687     h1->SetBinError(i,p->GetBinContent(i)*.001);
00688   }
00689   h1->SetXTitle(xtitle);
00690   return h1;
00691 }
00692 
00693 TH1D* Rms(TH3* h3,const char* basename, const char* basetitle,
00694               const char* slicename, Axis_t min, Axis_t max,
00695               const char* slicetype,
00696               const char* projecttype)
00697 {
00698   char name[200]; sprintf(name,"%s_prms",basename);
00699   TProfile* p = Profile(h3,name,basetitle,slicename,
00700                         min,max,slicetype,projecttype,"s");
00701   
00702   // make a 1 d histogram
00703   
00704   sprintf(name,"%s_h1",p->GetName());
00705   Int_t nBin= p->GetNbinsX();
00706   TH1D* h1 = new TH1D(name,p->GetTitle(),nBin,
00707                      p->GetXaxis()->GetBinLowEdge(1),
00708                      p->GetXaxis()->GetBinUpEdge(nBin));
00709 
00710   // now loop over the tprofile and fill the histogram
00711   for(int i=1;i<=nBin;i++){
00712     h1->SetBinContent(i,p->GetBinError(i));
00713     // set the error as 1% of the content
00714     h1->SetBinError(i,p->GetBinContent(i)*.001);
00715   }
00716   h1->SetXTitle(p->GetXaxis()->GetTitle());
00717   return h1;
00718 
00719 }
00720 
00721 void
00722 Printps(TCanvas* c1, const char* outDir, TString basename)
00723 {
00724 
00725   if(!basename.EndsWith(".ps")) basename += ".ps";
00726   TString s = outDir; s += "/"; s += basename;
00727   c1->Print(s.Data());
00728 
00729 }
00730 
00731 void
00732 Printgif(TCanvas* c1, const char* outDir, TString basename)
00733 {
00734 
00735   if(!basename.EndsWith(".gif")) basename += ".gif";
00736   TString s = outDir; s += "/"; s += basename;
00737   c1->Print(s.Data());
00738 
00739 }
00740 
00741 // errors are not reset...
00742 void
00743 Rebin(TH3* h3,int axisType,int nBin, Axis_t* ary)
00744 {
00745   if(!h3) { cout << "null pointer?" << endl; return; }
00746 
00747 
00748   // clone the histogram
00749   TH3* hClone = (TH3*)h3->Clone();
00750   TString name = h3->GetName(); name += "Clone";
00751 
00752   hClone->SetName(name.Data());
00753 
00754   TAxis *rebinAxis=0, *aAxis=0, *bAxis=0;
00755   TAxis *xAxis=hClone->GetXaxis(), 
00756     *yAxis=hClone->GetYaxis(), *zAxis=hClone->GetZaxis();
00757   TAxis *newAxis,*newAaxis,*newBaxis;
00758 
00759   switch(axisType){
00760   case 0: 
00761     rebinAxis=xAxis; aAxis=yAxis; bAxis=zAxis; 
00762     newAxis=h3->GetXaxis(); newAaxis=h3->GetYaxis(); newBaxis=h3->GetZaxis();
00763     break;
00764   case 1: 
00765     rebinAxis=yAxis; aAxis=xAxis; bAxis=zAxis; 
00766     newAxis=h3->GetYaxis(); newAaxis=h3->GetXaxis(); newBaxis=h3->GetZaxis();
00767     break;
00768   case 2: 
00769     rebinAxis=zAxis; aAxis=xAxis; bAxis=yAxis; 
00770     newAxis=h3->GetZaxis(); newAxis=h3->GetXaxis(); newBaxis=h3->GetYaxis();
00771     break;
00772   default: cout << "Wrong axis type: " << axisType << endl; exit(-1);
00773   }
00774 
00775   
00776   // set new axis of the original histogram
00777   cout << "old axis " << rebinAxis->GetNbins() << endl;
00778   newAxis->Set(nBin,ary);
00779   cout << "new axis " << newAxis->GetNbins() << endl;
00780 
00781   int aNbin=aAxis->GetNbins(), bNbin=bAxis->GetNbins();
00782   double* aAry = new double[aAxis->GetNbins()+1];
00783   double* bAry = new double[bAxis->GetNbins()+1];
00784 
00785   for(int i=0; i<aNbin; i++){
00786     int bin=i+1;
00787     aAry[i]=aAxis->GetBinLowEdge(bin);
00788   }
00789   aAry[aNbin] = aAxis->GetBinUpEdge(aNbin);
00790   newAaxis->Set(aNbin,aAry);
00791 
00792   for(int i=0; i<bNbin; i++){
00793     int bin=i+1;
00794     bAry[i]=bAxis->GetBinLowEdge(bin);
00795   }
00796   bAry[bNbin] = bAxis->GetBinUpEdge(bNbin);
00797   newBaxis->Set(bNbin,bAry);
00798 
00799   delete aAry; delete bAry;
00800 
00801   // now fill it
00802 
00803   double newval=0;
00804   int iBinXreal=0, iBinYreal=0, iBinZreal=0, bin=0;
00805   
00806   for(int aBin=1; aBin<=aAxis->GetNbins(); aBin++){
00807     for(int bBin=1; bBin<=bAxis->GetNbins(); bBin++){
00808       for(int newBin=1; newBin<=nBin; newBin++){
00809         newval=0;
00810         
00811         for(int rebinBin=1; rebinBin<=rebinAxis->GetNbins(); rebinBin++){
00812           if(newAxis->FindBin(rebinAxis->GetBinCenter(rebinBin))==newBin){
00813             switch(axisType){
00814             case 0: iBinXreal=rebinBin; iBinYreal=aBin; iBinZreal=bBin; break;
00815             case 1: iBinXreal=aBin; iBinYreal=rebinBin; iBinZreal=bBin; break;
00816             case 2: iBinXreal=aBin; iBinYreal=bBin; iBinZreal=rebinBin; break;
00817             }
00818             bin = hClone->GetBin(iBinXreal, iBinYreal, iBinZreal);
00819             newval += hClone->GetBinContent(bin);
00820           }
00821         }
00822         // set new bin content
00823         switch(axisType){
00824         case 0: iBinXreal=newBin; iBinYreal=aBin; iBinZreal=bBin; break;
00825         case 1: iBinXreal=aBin; iBinYreal=newBin; iBinZreal=bBin; break;
00826         case 2: iBinXreal=aBin; iBinYreal=bBin; iBinZreal=newBin; break;
00827         }
00828         bin = h3->GetBin(iBinXreal,iBinYreal,iBinZreal);
00829         h3->SetBinContent(bin,newval);
00830         
00831       } // newBin
00832     }
00833   }
00834   // reality check
00835   /*
00836   TH1* h1Clone=hClone->Project3D("x");
00837   TH1* h1New  =h3->Project3D("x");
00838 
00839   for(int i=1; i<=h1New->GetNbinsX(); i++){
00840     cout << "new bin " << i 
00841          << ",low = "<< h1New->GetXaxis()->GetBinLowEdge(i)
00842          << ",high = " << h1New->GetXaxis()->GetBinUpEdge(i)
00843          << ",val = " << h1New->GetBinContent(i) << endl;
00844     double sum=0;
00845     for(int j=1; j<=h1Clone->GetNbinsX(); j++){
00846       if(h1New->GetXaxis()->FindBin(h1Clone->GetXaxis()->GetBinCenter(j))==i){
00847         cout << "\told bin " << j
00848              << ",low = " << h1Clone->GetXaxis()->GetBinLowEdge(j)
00849              << ",high = " << h1Clone->GetXaxis()->GetBinUpEdge(j)
00850              << ",val = " << h1Clone->GetBinContent(j) << endl;
00851         sum+= h1Clone->GetBinContent(j);
00852       }
00853     }
00854     cout << "\tsum = " << sum << endl;
00855   }
00856   */
00857   
00858   // get rid of the clone
00859   delete hClone; delete aAry; delete bAry; 
00860   
00861 }
00862 //#endif
00863 
00864 //____________________
00865 
00866 void
00867 dump(TH1* h1)
00868 {
00869   cout << h1->GetName() << endl;
00870   Int_t nBin = h1->GetNbinsX();
00871 
00872   for(Int_t i=1; i<=nBin; i++){
00873     float errorFrac = (h1->GetBinContent(i)>0)?
00874       h1->GetBinError(i)/h1->GetBinContent(i) : 0;
00875     
00876     printf("\t%d : %.4f -- %.4f \t%.4f \t%.4f \t%.4f\n",
00877            i,
00878            h1->GetXaxis()->GetBinLowEdge(i),
00879            h1->GetXaxis()->GetBinUpEdge(i),
00880            h1->GetBinContent(i),
00881            h1->GetBinError(i),
00882            errorFrac);
00883   }
00884 }
00885 
00886 void 
00887 dump(TGraph* g)
00888 {
00889   cout << g->GetName() << endl;
00890   double *x = g->GetX();
00891   double *y = g->GetY();
00892   
00893   for(int i=0; i<g->GetN(); i++){
00894     printf("\t%d : %.4f \t %.4f\n",i,x[i],y[i]);
00895   }
00896 }
00897 
00898 void 
00899 dump(TGraphAsymmErrors* graph)
00900 {
00901   double* xValues = graph->GetX();
00902   double* yValues = graph->GetY();
00903   double* errXLow = graph->GetEXlow();
00904   double* errXHigh = graph->GetEXhigh();
00905   
00906   double* errYLow  = graph->GetEYlow();
00907   double* errYHigh = graph->GetEYhigh();
00908 
00909   for(int i=0; i<graph->GetN(); i++){
00910     printf("\t%d : %.4f < %.4f <%.4f \t %.3e %.3e\n",
00911            i,xValues[i]-errXLow[i],xValues[i],xValues[i]+errXHigh[i],
00912            yValues[i],errYLow[i]);
00913   }
00914 
00915 
00916 }
00917 
00918 
00919 void
00920 showHistogramValues(TH1D* h1, char* more)
00921 {
00922   *ofVal << "#################################################" << endl;
00923   *ofVal << h1->GetName();
00924   if(more) *ofVal << ":::::"<<more << endl;
00925   else *ofVal << endl;
00926 
00927   Int_t nBin = h1->GetNbinsX();
00928 
00929   Double_t errorFrac;
00930 
00931   // no under and over flows
00932 
00933   for(Int_t i=1; i<=nBin; i++){
00934     errorFrac = (h1->GetBinContent(i)>0)?
00935       h1->GetBinError(i)/h1->GetBinContent(i) : 0;
00936 
00937     *ofVal << "\t" << h1->GetXaxis()->GetBinLowEdge(i) << " < pt < "
00938            << h1->GetXaxis()->GetBinUpEdge(i) << "    "
00939            << "\t" << h1->GetBinContent(i) << "  \t" << h1->GetBinError(i)
00940            << "\t" << errorFrac << endl;
00941   }
00942 
00943 }
00944 
00945 //____________________
00946 
00947 void 
00948 showTGraphValues(TGraphAsymmErrors* graph,char* more)
00949 {
00950   *ofVal << "#################################################" << endl;
00951   *ofVal << graph->GetName() << ":::::" << more << endl;
00952   
00953   const Int_t nBin = graph->GetN();
00954   Double_t* xValues = graph->GetX();
00955   Double_t* yValues = graph->GetY();
00956   Double_t* errYHigh = graph->GetEYhigh();
00957   Double_t* errXLow = graph->GetEXlow();
00958   Double_t* errXHigh = graph->GetEXhigh();
00959 
00960   Double_t lowEdge, upEdge, errorFrac;
00961 
00962   for(Int_t i=0; i<nBin; i++){
00963     lowEdge = xValues[i] - errXLow[i];
00964     upEdge  = xValues[i] + errXHigh[i];
00965     
00966     errorFrac = errYHigh[i]/yValues[i];
00967 
00968     *ofVal << "\t" << lowEdge << "<" << xValues[i] << "<" << upEdge
00969            << "\t" << yValues[i] << "\t" << errYHigh[i] << "\t" << errorFrac 
00970            << endl;
00971   }
00972 }
00973 // for each x slice, divide the y bins by the integral of the x slice
00974 void scale2D(TH2* h2)
00975 {
00976   TAxis *xAxis=h2->GetXaxis();
00977   TAxis *yAxis=h2->GetYaxis();
00978 
00979   for(int ix=1; ix<=xAxis->GetNbins(); ix++){
00980     // first get the sum of this slice
00981     float sum=0;
00982     for(int iy=1; iy<=yAxis->GetNbins(); iy++){
00983       int bin=h2->GetBin(ix,iy);
00984       sum += h2->GetBinContent(bin);
00985     }
00986     // now scale
00987     for(int iy=1; iy<=yAxis->GetNbins(); iy++){
00988       int bin=h2->GetBin(ix,iy);
00989       if(sum)
00990         h2->SetBinContent(bin,h2->GetBinContent(bin)/sum);
00991     }
00992   }
00993 }