///////////////////////////////////////////////////////////////////
/////// Convolution Function                             //////////
///////                                                  //////////
/////// Written by Kim Alwyn  (09/03/2006)               //////////
///////                                                  //////////
/////// Constructor requires: name of convolution plot,  //////////
///////                       title of convolution plot, //////////
///////                       name of h1 (Breit-Wigner), //////////
///////                       name of h2 (mass spectrum).//////////
/////// This function draws a convolution of the two     //////////
/////// inputted functions.                              //////////
///////                                                  //////////
///////////////////////////////////////////////////////////////////

int Convolution(const char *name, const char *title, const TH1 *h1, const 
TH1 *h2)
{
  
  //Check that both histograms exist.
  
  if (!h1){  
    Error("Convolute", "Breit-Wigner Histogram does not exist");
  }
  if (!h2){
    Error("Convolute", "Mass Spectrum Histogram does not exist");
  }
  
  //Find the number of bins and the range of the Breit-Wigner histogram:
  
  int binNum1 = h1->GetNbinsX();
  float lowerEdge1 = h1->GetBinLowEdge(1);
  float upperEdge1 = h1->GetBinLowEdge(binNum1)+h1->GetBinWidth(binNum1);

  //Find the number of bins and the range of the mass spectrum histogram:
  
  int binNum2 = h2->GetNbinsX();
  float lowerEdge2 = h2->GetBinLowEdge(1);
  float upperEdge2 = h2->GetBinLowEdge(binNum2)+h2->GetBinWidth(binNum2);

  //Find new upper and lower limits: These are the same as for h2 but 
  //binNum1/2 extra bins added onto each end. 
  
  float lowerEdge3 = lowerEdge2 - (0.5*binNum1*h2->GetBinWidth(binNum2));
  float upperEdge3 = upperEdge2 + (0.5*binNum1*h2->GetBinWidth(binNum2));

  //Book the final histogram: using the new limits found above.  
  
  TH1F *h4 = new TH1F( name, title, 
		       binNum2+binNum1, lowerEdge3, upperEdge3);
    
  //i = bin number of h1. There are binNum1 iterations of this loop.
  //j = bin number of h2. There are binNum2 iterations of this loop.
  
  //This loop iterates over each bin in the Breit-Wigner histogram, (h1).
  //It transposes the mass spectrum to centre on the bin and scales it to be 
  //the same height. It is then added cumulatively to the final histogram,(h4).
  
  for (int i=1; i<binNum1+1; i++) {

    //Find the bin corresponding to 0.0 on the x-axis for each histogram.
    //(i.e.the point the graphs are centred on)
    
    //    int midBin1 = h1->FindBin(0.);
    //    int midBin2 = h2->FindBin(0.); 
   
    //Book a histogram identical to the final (h4).
    
    TH1F *h3 = new TH1F(  name, title,
			  binNum2+binNum1, lowerEdge3, upperEdge3);
    
    //Transpose the mass spectrum along by 'i-midBin2+midBin1' in x-direction.
    //This is done by looping over each bin in the Mass Spectrum histogram(h2)
    //and moving it onto h3 in the transposed position.
    //("midBin2" allows for the fact that h1 and h2 may have different scales.
    //"midBin1" compensates for the fact that h2 and h3 have different 
    //lowerEdge values.)
    
    for (int j=1; j<binNum2+1; j++) {
      Double_t tmp = h2->GetBinContent(j);
      h3->SetBinContent(j+i/*-(0.5*midBin2)+midBin1*/,tmp);
      
    } //end of for loop over j
    
    //Find the height of the ith bin of the Breit-Wigner function(h1):
    
    Double_t binH = h1->GetBinContent(i);  
    
    //Find the scale factor needed to scale the highest bin in the Mass 
    //Spectrum (h2) up to the height of the ith bin of the Breit-Wigner (h1).  
 
    Double_t scale = 0.0;
    if (h3->GetBinContent(h3->GetMaximumBin())!=0){
      scale = binH / h3->GetBinContent(h3->GetMaximumBin());
    }
      
    //Add h3 to the final histogram (h4), scaling it by "scale". Then delete h3
    //ready for it to be refilled in the next iteration of the loop.
    
    h4->Add(h3, scale);
    h3->Delete();
       
    }//end of for loop over i
  
  //Draw the final histogram.
  h4->Draw();
  
}