loon: MadContour.cxx Source File

00001 #include "MadContour.h"
00002 
00003 MadContour::MadContour(){
00004   fdraw = false;
00005   finit = false;
00006   fMinChi2 = 1000000.;
00007   fMinPar1 = 0;
00008   fMinPar2 = 0;
00009   fContCount = 0;
00010   fCloseLoop = true;
00011 }
00012 
00013 MadContour::MadContour(TH2F *Chi2Hist){
00014   fdraw = false;
00015   finit = false;
00016   fMinChi2 = 1000000.;
00017   fMinPar1 = 0;
00018   fMinPar2 = 0;
00019   fContCount = 0;
00020   Init(SetChiMap(Chi2Hist));
00021   fCloseLoop = true;
00022 }
00023 
00024 MadContour::~MadContour(){
00025   if(fdraw){
00026     delete fcen;
00027     delete fgr1;
00028     delete fgr2;
00029     delete fgr3;
00030     delete femap;
00031   }
00032 }
00033 
00034 int MadContour::SetChiMap(TH2F *Chi2Hist){
00035 
00036   finit = false;
00037     
00038   if(!Chi2Hist){
00039     std::cerr << "No TH2F" << std::endl;
00040     return 0;
00041   }
00042 
00043   if(Chi2Hist->Integral()==0){
00044     std::cerr << "TH2F has no entries" << std::endl;
00045     return 0;
00046   }
00047   
00048   fChi2Hist = Chi2Hist;
00049   return 1;
00050 
00051 }
00052 
00053 TH2F *MadContour::GetEmptyMap(){
00054   
00055   TH2F *hist = new TH2F("emptymap","Confidence Limits",
00056                         fChi2Hist->GetNbinsX(),
00057                         fChi2Hist->GetBinLowEdge(1),
00058                         fChi2Hist->GetBinCenter(fChi2Hist->GetNbinsX())
00059                         +fChi2Hist->GetBinWidth(1),
00060                         fChi2Hist->GetNbinsY(),
00061                         fChi2Hist->GetYaxis()->GetBinLowEdge(1),
00062                         fChi2Hist->GetYaxis()
00063                         ->GetBinCenter(fChi2Hist->GetNbinsY())
00064                         +fChi2Hist->GetYaxis()->GetBinWidth(1));
00065   
00066   return hist;
00067 
00068 }
00069 
00070 void MadContour::Init(int use){
00071 
00072   if(use==0) return;
00073 
00074   int minbinpar1 = 0;
00075   int minbinpar2 = 0;
00076 
00077   float binwidthpar1 = fChi2Hist->GetXaxis()->GetBinWidth(1);
00078   float binwidthpar2 = fChi2Hist->GetYaxis()->GetBinWidth(1);
00079   
00080   for(int i=1;i<=fChi2Hist->GetNbinsX();i++){
00081     for(int j=1;j<=fChi2Hist->GetNbinsY();j++){
00082       if(fChi2Hist->GetBinContent(i,j)<fMinChi2) {
00083         fMinChi2 = fChi2Hist->GetBinContent(i,j);
00084         fMinPar1 = fChi2Hist->GetXaxis()->GetBinCenter(i);
00085         minbinpar1 = i;
00086         fMinPar2 = fChi2Hist->GetYaxis()->GetBinCenter(j);
00087         minbinpar2 = j;
00088       }
00089     }
00090   }
00091   
00092   if(minbinpar1!=1&&minbinpar1!=fChi2Hist->GetNbinsX()) 
00093     fMinPar1 =  
00094       (fChi2Hist->GetBinContent(minbinpar1-1,minbinpar2)-fMinChi2)*binwidthpar1
00095       /(fChi2Hist->GetBinContent(minbinpar1+1,minbinpar2)
00096         + fChi2Hist->GetBinContent(minbinpar1-1,minbinpar2) - 2.*fMinChi2) 
00097       - binwidthpar1/2. + fMinPar1;
00098   
00099   if(minbinpar2!=1&&minbinpar2!=fChi2Hist->GetNbinsY()) 
00100     fMinPar2 =
00101       (fChi2Hist->GetBinContent(minbinpar1,minbinpar2-1)-fMinChi2)*binwidthpar2
00102       /(fChi2Hist->GetBinContent(minbinpar1,minbinpar2+1)
00103         + fChi2Hist->GetBinContent(minbinpar1,minbinpar2-1) - 2.*fMinChi2) 
00104       - binwidthpar2/2. + fMinPar2;
00105   
00106   finit = true;
00107 
00108 }
00109 
00110 void MadContour::CloseLoop(bool closeloop)
00111 {
00112   fCloseLoop = closeloop;
00113 }
00114 
00115 TGraph *MadContour::GetCentre(){
00116   
00117   float plotpar1[1];
00118   float plotpar2[1];
00119   plotpar1[0] = fMinPar1;
00120   plotpar2[0] = fMinPar2;
00121   TGraph *centre = new TGraph(1,plotpar1,plotpar2);
00122   centre->SetName("centre");
00123   return centre;
00124   
00125 }
00126 
00127 TGraph *MadContour::GetContour(float limit){
00128 
00129   if(limit==67||limit==68.27)         limit=2.30;
00130   else if(limit==90)                  limit=4.61;
00131   else if(limit==95||limit==95.45)    limit=6.18;
00132   else if(limit==99)                  limit=9.21;
00133   else if(limit==99.73)               limit=11.83;
00134   else                                limit=4.61;
00135   
00136   fContCount+=1;
00137 
00138   float *plotpar1 = new float[4*fChi2Hist->GetNbinsX()];
00139   float *plotpar2 = new float[4*fChi2Hist->GetNbinsY()];
00140   int count = 0;
00141   
00142   char histname[256];
00143   for(int i = fChi2Hist->GetNbinsX();i>=1;i--){
00144     
00145     sprintf(histname,"temp%i",i);
00146     TH1F *temphist = (TH1F*) fChi2Hist->ProjectionY(histname,i,i);
00147     
00148     if(temphist->GetMinimum()<fMinChi2+limit){
00149       
00150       int minbin = temphist->GetMinimumBin();
00151       int theBin = 1;
00152       
00153       while(minbin-theBin>1){
00154         
00155         if(temphist->GetBinContent(minbin-theBin)>=fMinChi2+limit){
00156           
00157           while(temphist->GetBinContent(minbin-theBin)>=fMinChi2+limit
00158                 && minbin-theBin>1) theBin+=1;
00159           
00160           if(temphist->GetBinContent(minbin-theBin)<fMinChi2+limit){
00161             plotpar1[count] = fChi2Hist->GetXaxis()->GetBinCenter(i);
00162             if((temphist->GetBinContent(minbin-theBin) -
00163                 temphist->GetBinContent(minbin-theBin+1)) != 0) {
00164               plotpar2[count] = temphist->GetBinCenter(minbin-theBin+1) + 
00165                 (temphist->GetBinCenter(minbin-theBin) - 
00166                  temphist->GetBinCenter(minbin-theBin+1)) *
00167                 ((fMinChi2+limit) - temphist->GetBinContent(minbin-theBin+1)) /
00168                 (temphist->GetBinContent(minbin-theBin) -
00169                  temphist->GetBinContent(minbin-theBin+1));
00170             }
00171             else plotpar2[count] = (minbin-theBin+1);
00172             count += 1;
00173           }
00174           else break;
00175         }
00176         
00177         while(temphist->GetBinContent(minbin-theBin)<fMinChi2+limit
00178               && minbin-theBin>1) theBin+=1;
00179 
00180         plotpar1[count] = fChi2Hist->GetXaxis()->GetBinCenter(i);
00181         if((temphist->GetBinContent(minbin-theBin) -
00182             temphist->GetBinContent(minbin-theBin+1)) != 0){
00183           plotpar2[count] = temphist->GetBinCenter(minbin-theBin+1) + 
00184             (temphist->GetBinCenter(minbin-theBin) - 
00185              temphist->GetBinCenter(minbin-theBin+1)) *
00186             ((fMinChi2+limit) - temphist->GetBinContent(minbin-theBin+1)) /
00187             (temphist->GetBinContent(minbin-theBin) -
00188              temphist->GetBinContent(minbin-theBin+1));
00189         }
00190         else plotpar2[count] = temphist->GetBinCenter(minbin-theBin+1);
00191         count += 1;
00192       }
00193     }
00194     delete temphist;
00195   }
00196   
00197   for(int i = 1;i<=fChi2Hist->GetNbinsX();i++){
00198 
00199     sprintf(histname,"temp%i",i);
00200     TH1F *temphist = (TH1F*) fChi2Hist->ProjectionY(histname,i,i);
00201 
00202     if(temphist->GetMinimum()<fMinChi2+limit){
00203 
00204       int minbin = temphist->GetMinimumBin();
00205       int theBin = 1;
00206 
00207       while(minbin+theBin<temphist->GetNbinsX()){
00208         
00209         if(temphist->GetBinContent(minbin+theBin)>=fMinChi2+limit){
00210           
00211           while(temphist->GetBinContent(minbin+theBin)>=fMinChi2+limit
00212                 && minbin+theBin<temphist->GetNbinsX()) theBin+=1;
00213           
00214           if(temphist->GetBinContent(minbin+theBin)<fMinChi2+limit){
00215             plotpar1[count] = fChi2Hist->GetXaxis()->GetBinCenter(i);
00216             if((temphist->GetBinContent(minbin+theBin) -
00217                 temphist->GetBinContent(minbin+theBin-1)) != 0) {
00218               plotpar2[count] = temphist->GetBinCenter(minbin+theBin-1) + 
00219                 (temphist->GetBinCenter(minbin+theBin) - 
00220                  temphist->GetBinCenter(minbin+theBin-1)) *
00221                 ((fMinChi2+limit) - temphist->GetBinContent(minbin+theBin-1)) /
00222                 (temphist->GetBinContent(minbin+theBin) -
00223                  temphist->GetBinContent(minbin+theBin-1));
00224             }
00225             else plotpar2[count] = temphist->GetBinCenter(minbin+theBin-1);
00226             count += 1;
00227           }
00228           else break;
00229         }
00230         
00231         while(temphist->GetBinContent(minbin+theBin)<fMinChi2+limit 
00232               && minbin+theBin<temphist->GetNbinsX()) theBin+=1;
00233 
00234         plotpar1[count] = fChi2Hist->GetXaxis()->GetBinCenter(i);
00235         if((temphist->GetBinContent(minbin+theBin) - 
00236             temphist->GetBinContent(minbin+theBin-1)) != 0) {
00237         plotpar2[count] = temphist->GetBinCenter(minbin+theBin-1) + 
00238           (temphist->GetBinCenter(minbin+theBin) - 
00239            temphist->GetBinCenter(minbin+theBin-1)) *
00240           ((fMinChi2+limit) - temphist->GetBinContent(minbin+theBin-1)) /
00241           (temphist->GetBinContent(minbin+theBin) -
00242            temphist->GetBinContent(minbin+theBin-1));
00243         }
00244         else plotpar2[count] = temphist->GetBinCenter(minbin+theBin-1);
00245         count+=1;
00246       }
00247     }
00248     delete temphist;
00249   }
00250 
00251   //close loop
00252   if(fCloseLoop){
00253     plotpar1[count] = plotpar1[0];
00254     plotpar2[count] = plotpar2[0];
00255     count+=1;
00256   }
00257 
00258   TGraph *gr = new TGraph(count,plotpar1,plotpar2);
00259   char ContName[256];
00260   sprintf(ContName,"contour%i",fContCount);
00261   gr->SetName(ContName);
00262   gr->SetLineWidth(3);
00263   gr->SetLineColor(fContCount+1);
00264 
00265   delete plotpar1;
00266   delete plotpar2;
00267   
00268   return gr;
00269 
00270 }
00271 
00272 float MadContour::GetErrorOnBestPar(float limit,int par,int Dim)
00273 {
00274 
00275   if(Dim==2){ //considering size of contour in 2D
00276     if(limit==67||limit==68.27)         limit=2.30;
00277     else if(limit==90)                  limit=4.61;
00278     else if(limit==95)                  limit=5.99;
00279     else if(limit==95.45)               limit=6.18;
00280     else if(limit==99)                  limit=9.21;
00281     else if(limit==99.73)               limit=11.83;
00282     else                                limit=4.61;
00283   }
00284   else if(Dim==1){ //considering 1-D errors i.e. deltachi2=1 <=> 1 sigma error
00285     if(limit==67||limit==68.27)         limit=1.00;
00286     else if(limit==90)                  limit=2.71;
00287     else if(limit==95)                  limit=3.84;
00288     else if(limit==95.45)               limit=4.00;
00289     else if(limit==99)                  limit=6.63;
00290     else if(limit==99.73)               limit=9.00;
00291     else                                limit=2.71;    
00292   }
00293   else {
00294     limit = 4.61; //assume 90%, 2D
00295   }
00296 
00297   int bin = 1;
00298   TH1D *temphist;
00299   if(par==1){
00300     bin = fChi2Hist->GetYaxis()->FindBin(fMinPar2);
00301     temphist = fChi2Hist->ProjectionX("temphist",bin,bin);
00302   }
00303   else if(par==2){
00304     bin = fChi2Hist->GetXaxis()->FindBin(fMinPar1);
00305     temphist = fChi2Hist->ProjectionY("temphist",bin,bin);
00306   }
00307   else {
00308     //Invalid parameter selection
00309     return -1;
00310   }
00311 
00312   int minbin = temphist->GetMinimumBin();
00313 
00314   int theBin = 1;
00315   while(temphist->GetBinContent(minbin-theBin)<fMinChi2+limit
00316         && minbin-theBin>1) theBin+=1;
00317   
00318   float LowerErr = temphist->GetBinCenter(minbin-theBin+1) + 
00319     (temphist->GetBinCenter(minbin-theBin) - 
00320      temphist->GetBinCenter(minbin-theBin+1)) *
00321     ((fMinChi2+limit) - temphist->GetBinContent(minbin-theBin+1)) /
00322     (temphist->GetBinContent(minbin-theBin) -
00323      temphist->GetBinContent(minbin-theBin+1));
00324 
00325   theBin = 1;
00326   while(temphist->GetBinContent(minbin+theBin)<fMinChi2+limit 
00327         && minbin+theBin<temphist->GetNbinsX()) theBin+=1;
00328 
00329   float UpperErr = temphist->GetBinCenter(minbin+theBin-1) + 
00330     (temphist->GetBinCenter(minbin+theBin) - 
00331      temphist->GetBinCenter(minbin+theBin-1)) *
00332     ((fMinChi2+limit) - temphist->GetBinContent(minbin+theBin-1)) /
00333     (temphist->GetBinContent(minbin+theBin) -
00334      temphist->GetBinContent(minbin+theBin-1));
00335   
00336   delete temphist;
00337 
00338   return UpperErr-LowerErr;
00339 
00340 }
00341 
00342 void MadContour::Draw()
00343 {
00344 
00345   if(!finit) {
00346     std::cerr << "ChiMap cannot be used" << std::endl;
00347     return;
00348   }
00349 
00350   fdraw = true;
00351   femap = GetEmptyMap();
00352   fcen = GetCentre();
00353   fgr1 = GetContour(90);
00354   fgr2 = GetContour(95);
00355   fgr3 = GetContour(99);
00356 
00357   femap->Draw();
00358   fcen->Draw("*");
00359   fgr1->Draw("L");
00360   fgr2->Draw("L");
00361   fgr3->Draw("L");
00362 
00363 }
00364 
00365 void MadContour::ContDraw()
00366 {
00367   double lev[4] = {2.30,4.61,6.18,9.21};
00368   lev[0]+=fChi2Hist->GetMinimum();
00369   lev[1]+=fChi2Hist->GetMinimum();
00370   lev[2]+=fChi2Hist->GetMinimum();
00371   lev[3]+=fChi2Hist->GetMinimum();
00372   fChi2Hist->SetContour(4,lev);
00373   fChi2Hist->Draw("CONT1");
00374   fcen = GetCentre();
00375   fcen->Draw("*");  
00376 }