dijetcor.C

const double kJPi = TMath::Pi();
const double kJTwoPi = kJPi*2;
void dijetcor(){

  TH1D *hCor = new TH1D ( "hCor", "hCor", 300, -1/3.0, ( 2-1/3.0 ) ); 
  for( int ievt=0;ievt<1000;ievt++ ){
    double phi1 = gRandom->Uniform( kJTwoPi );
    //double phi2 = phi1 + kJPi ;
    double phi2 = gRandom->Gaus( phi1 + kJPi, 0.3 ) ;
    double dphi1 = 0.1*kJPi;
    double dphi2 = 0.1*kJPi;

    int n1 = 5;
    int n2 = 5;

    vector<double> t1;
    for( int i=0;i<n1;i++ ){
      //t1.push_back( gRandom->Gaus(phi1, dphi1) );
      t1.push_back( gRandom->Uniform(dphi1)-dphi1/2+phi1 );
    }
    for( int j=0;j<n2;j++ ){
      //t1.push_back( gRandom->Gaus(phi2, dphi2) );
      t1.push_back( gRandom->Uniform(dphi2)-dphi2/2+phi2 );
    }

    for( unsigned int ii=0;ii<t1.size();ii++ ){
      for( unsigned int jj=0;jj<t1.size();jj++ ){
        if( ii == jj ) continue;
        hCor->Fill( DeltaPhi( t1[ii], t1[jj] )/kJPi );
      }
    }
  }
  hCor->Draw();

}

double DeltaPhi(double phi1, double phi2) {
  double res =  atan2(sin(phi1-phi2), cos(phi1-phi2));
  // if (res<0) res *= -1.0;
  return res>-kJPi/3.0 ? res : kJTwoPi+res ; 
}