#include "T49HbtMomentumResolutionCorrection.h"
ClassImp(T49HbtMomentumResolutionCorrection)
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//This class produces a 3d histogramms containing the momentum resolution correction. The measred correlation function
//should be multiplied with this histogramm to get the momentum.
//This class fills takes as input particle lists. It combines pairs and fills four histogramms, 
//ideal signal distbritution (pairs are weight according to BE and coulomb), ideal background (weight =1)
//and smeared signal and background. The smearing is either done by default momentum resolution or by
//resolution obtain from Monte Carlo (provided via the T49HbtMomentumResolutionProvider class).
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 T49HbtMomentumResolutionCorrection::T49HbtMomentumResolutionCorrection(Int_t nKt, Double_t minKt, Double_t maxKt, 
								       Int_t nY, Double_t minY, Double_t maxY,
								       Int_t nx, Float_t xmin, Float_t xmax, 
								       Int_t ny, Float_t ymin, Float_t ymax,
								       Int_t nz, Float_t zmin, Float_t zmax)
{
  
  // Init size of TH3D's
  mMinQside = xmin ;
  mMaxQside = xmax ;
  mMinQout  = ymin ;
  mMaxQout  = ymax ;
  mMinQlong = zmin ;
  mMaxQlong = zmax ;
  
  // Init Kt-Y bins
  mNKtBins = nKt ;
  mKtBins = new Double_t[mNKtBins+1] ;
  for(Int_t index=0 ; index<(mNKtBins+1) ; index++)
    {
      mKtBins[index] = ((maxKt-minKt)/mNKtBins) * (Double_t) index + minKt;
    }
  // Y bins
  mNYBins = nY ;
  mYBins  = new Double_t[mNYBins+1] ;
  for(Int_t index=0 ; index<(mNYBins+1) ; index++)
    {
      mYBins[index]= ((maxY-minY)/mNYBins) * (Double_t) index + minY;
    }
  // Max Number of KtY bins
  mMaxKtY = mNKtBins  *  mNYBins ;
  
  //
  // Set default init values
  //
  // HBT parameters
  mLambda = 0.8 ;
  mRside  = 5.0 ;
  mRout   = 6.0 ;
  mRlong  = 7.0 ;
  // Purity
  mPurity = 0.5 ;
  // CoulombCalculator
  mCoulombCalculator = new T49HbtCoulomb() ;
  mCoulombRadius = 10 ; 
  mCoulombCalculator->SetMeanSeparation(mCoulombRadius) ;
  mIncludeCoulombEffect = 1 ;
  
  // Init constants
  mPionMass = 0.139 ;
  mHbarC2 = 0.0389379 ;

  // Init radom generator
  mRand = new TRandom() ;

  // Per default do not use momentum resolution from MC events
  mMomentumResolutionProvider  = NULL ;
  mTotalNumberOfInquries       = NULL ;
  mNumberBadInquiries          = NULL ; 

  // Set default momentum resolution
  mDefaultResolutionQside = 0.01 ;
  mDefaultResolutionQout  = 0.01 ;
  mDefaultResolutionQlong = 0.01 ;

  //
  // Create 5 Histos for signal and background, smeared and unsmeared, and finally the correction
  //
  mNHistos = 5 ;
  // Names for ideal, smeared and correction distributions
  TString lIdealSignal("IdealSignal") ;
  TString lIdealBackground("IdealBackground") ;
  TString lSmearedSignal("SmearedSignal") ;
  TString lSmearedBackground("SmearedBackground") ;
  TString lCorrection("Correction") ;
  
  mHistoArray = new TObjArray(mMaxKtY*mNHistos) ;
  for(Int_t lKtYIndex = 0 ; lKtYIndex < mMaxKtY ; lKtYIndex++)
    {
      for(Int_t lHistoIndex = 0 ; lHistoIndex < mNHistos ; lHistoIndex++)
	{
	  TString lName("Histo_KtYBin") ;
	  lName += lKtYIndex ;
	  
	  switch (lHistoIndex) 
	    {
	    case 0:
	      lName.Append(lIdealSignal) ;
	      break ;
	    case 1:
	      lName.Append(lIdealBackground) ;
	      break ;
	    case 2:
	      lName.Append(lSmearedSignal) ;
	      break ;
	    case 3:
	      lName.Append(lSmearedBackground) ;
	      break ;
	    case 4 :
	      lName.Append(lCorrection) ;
	      break ;
	    default:
	      cerr << "Serious error: should never occur !" << endl ; 
	      break ;
	    }
	  // Create Histogramm
	  TH3D* lHistogram = new TH3D(lName.Data(),lName.Data(),nx,xmin,xmax,ny,ymin,ymax,nz,zmin,zmax);
	  
	  // Store it in TObjArray
	  mHistoArray->AddAt(lHistogram,(lKtYIndex*mNHistos)+lHistoIndex) ;
	}
    }
}
///////////////////////////////////////////////////////////////////////////////////////
 void T49HbtMomentumResolutionCorrection::SetMomentumResolutionProvider(T49HbtMomentumResolutionProvider* aMomentumResolutionProvider)
{
  mMomentumResolutionProvider  = aMomentumResolutionProvider ;
  if(!aMomentumResolutionProvider)
    {
      cerr << "Warning: Object supposed to contain momentum resolution is empty ! " << endl; 
    }
  // Arrange for some statistics
  mTotalNumberOfInquries = new Double_t[mMaxKtY] ;
  mNumberBadInquiries =  new Double_t[mMaxKtY] ;
  
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 void T49HbtMomentumResolutionCorrection::SetHBTRadii(Double_t Lambda, Double_t Rside, Double_t Rout, Double_t Rlong)
{
  mLambda = Lambda ;
  mRside  = Rside ;
  mRout   = Rout ;
  mRlong  = Rlong ;
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 void T49HbtMomentumResolutionCorrection::SetCoulombRadius(Double_t CoulombRadius)
{
  mCoulombRadius = CoulombRadius ; 
  mCoulombCalculator->SetMeanSeparation(mCoulombRadius) ;
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 void T49HbtMomentumResolutionCorrection::SetPurity(Double_t Purity)
{
   mPurity = Purity ;
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 void T49HbtMomentumResolutionCorrection::Fill(TObjArray* ListA, TObjArray* ListB, T49Cut2Track* T2Cut)
{
  // Make sure we get two differen lists ! We need to use the uncorrelated background.
  if(ListA==ListB)
    {
      cerr << "We need two set of tracks from DIFFERENT events !" << endl ;
    }
    
  T49ParticleRoot *Part1;
  T49ParticleRoot *Part2;
  
  Int_t IP1 = 0;
  Int_t IP2 = 0;

  // Loop over ListA
  IP1 = 0;
  TIter NextParticle(ListA);
  while ((Part1 = ((T49ParticleRoot *) NextParticle()))) 
    { 
      IP1++; 
      IP2 = 0;
      TIter NextP2(ListB);
      // Loop over ListB
      while ((Part2 = ((T49ParticleRoot *) NextP2()))) 
	{ 
	  IP2++; 
	  if (Part1 == Part2) 
	    {
	      cerr << "Error: Identical track in two different lists! Should never happen!" << endl  ;
	    }
	  if ((T2Cut == NULL) || (T2Cut->CheckPair(Part1,Part2))) 
	    {
	      if ((Part1->GetPx() == Part2->GetPx()) && 
		  (Part1->GetPy() == Part2->GetPy()) && 		  
		  (Part1->GetPz() == Part2->GetPz()) ) 
		{
		  cerr << "Error: Identical track in two different lists! Should never happen!" << endl ;
		}
	      Fill(Part1,Part2);
	    }
	}
    }
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 Int_t  T49HbtMomentumResolutionCorrection::GetKtYIndex(Double_t lKt, Double_t lYpp)
{
  // kt
  Int_t lKtIndex = 0 ;
  if( lKt < mKtBins[0] || lKt > mKtBins[mNKtBins] )
    { 
       lKtIndex = -1 ;
    } 
  else 
    { 
       while( lKt > mKtBins[lKtIndex+1] )
	 {
	   lKtIndex++;
	 } 
    }
  // Y
  Int_t lYIndex = 0 ;
  if( lYpp < mYBins[0] || lYpp > mYBins[mNYBins] )
    { 
      lYIndex = -1 ;
    } 
  else 
    { 
      while( lYpp > mYBins[lYIndex+1] )
	{
	  lYIndex++;
	} 
    }
  // result
  return (lKtIndex+lYIndex*mNKtBins) ;

}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 void T49HbtMomentumResolutionCorrection::Fill(T49ParticleRoot* t1, T49ParticleRoot* t2)
{
  //
  // Fill pairs into TH3Ds: weighted and unweighted, smeared and unsmeared
  //
  // Calculate Kt-Y bin
  Double_t lKt  = 0.5 * TMath::Sqrt(T49Tool::Px(t1,t2)*T49Tool::Px(t1,t2) + T49Tool::Py(t1,t2)*T49Tool::Py(t1,t2)) ;
  Double_t lYpp = 0.5 * TMath::Log((t1->GetE(mPionMass)+t2->GetE(mPionMass)
				    +t1->GetPz()+t2->GetPz())
				   /(t1->GetE(mPionMass)+t2->GetE(mPionMass)
				     - t1->GetPz()-t2->GetPz()));
  Int_t lKtYIndex = GetKtYIndex(lKt,lYpp) ;
  
  if ((lKtYIndex > -1) && (lKtYIndex < mMaxKtY)) 
    { 
      // Calculate qout,qside and qlong
      Double_t lQside = 0.0 ;
      Double_t lQout  = 0.0 ;
      Double_t lQlong = 0.0 ;
      QCalc(t1,t2, lQside, lQout, lQlong) ;

      if( lQside > mMinQside && lQside < mMaxQside && 
	  lQout  > mMinQout  && lQout  < mMaxQout  &&
	  lQlong > mMinQlong && lQlong < mMaxQlong)
	{
	  //
	  // Ideal Distributions
	  // 

	  // Calculate weight for signal: Weight = purity * [ C_BoseEinstein(qSide,qOut,qLong;Rout,Rside,Rlong) * Ccoulomb(qinv;Rinv ] + (1-p)
	  Double_t mBoseEinsteinWeight = 
	    1.0 + mLambda * exp(((-1)*(lQout*lQout*mRout*mRout) +(-1)*(lQside*lQside*mRside*mRside)+(-1)*(lQlong*lQlong*mRlong*mRlong) )/ mHbarC2 ) ; 
	  Double_t lQinv = T49Tool::Qinv(t1,t2,mPionMass,mPionMass) ;
	  Double_t mCoulombWeight = 1.0 ;
	  if(mIncludeCoulombEffect==1)
	    {
	      mCoulombWeight =  mCoulombCalculator->Weight(lQinv) ;
	    }
	  Double_t mWeight = mPurity *  (mBoseEinsteinWeight * mCoulombWeight) + (1-mPurity) ;
	  	  
	  // Fill IdealSignal
	  ((TH3D*) mHistoArray->At(lKtYIndex*mNHistos))->Fill(lQside,lQout,lQlong,mWeight) ;

	  // Background
	  ((TH3D*) mHistoArray->At(lKtYIndex*mNHistos+1))->Fill(lQside,lQout,lQlong,1.0) ;
     
	  //
	  // Smeared distributions
	  // 
	  // Smear momentum 
	  Double_t MeanQside = 0.0 ;
	  Double_t RMSQside  = 0.0 ;
	  Double_t MeanQout  = 0.0 ;
	  Double_t RMSQout   = 0.0 ;
	  Double_t MeanQlong = 0.0 ;
	  Double_t RMSQlong  = 0.0 ;
	
	  // If defined, get momentum resolution from "database"
	  if(mMomentumResolutionProvider)
	    {
	      mMomentumResolutionProvider->GetMeanAndRMS(lKt,lYpp,lQside,lQout,lQlong,MeanQside,MeanQout,MeanQlong,RMSQside,RMSQout,RMSQlong) ;
	      mTotalNumberOfInquries[lKtYIndex]++;
	      if(RMSQside == 0.0 && RMSQout == 0.0 && RMSQlong == 0.0 && MeanQside== 0.0 && MeanQout == 0.0 && MeanQlong == 0.0)
		{
		  mNumberBadInquiries[lKtYIndex]++;
		}
	    }

	  // Qside
	  Double_t lRandomShiftQside = 0.0 ;
	  if(RMSQside>0.0001 && mMomentumResolutionProvider)
	    {
	      lRandomShiftQside  = mRand->Gaus(MeanQside,RMSQside) ;
	    }
	  else
	    {
	      lRandomShiftQside  = mRand->Gaus(0.0,mDefaultResolutionQside) ;
	    }
	  Double_t lQsideSmeared = lQside + lRandomShiftQside ;
	  if(lQsideSmeared<0.0)
	    {
	      // To ensure symmetry throw these rare cases away
	      return ;
	    }

	  // Qout
	  Double_t lRandomShiftQout = 0.0 ;
	  if(RMSQout>0.0001 && mMomentumResolutionProvider)
	    {
	      lRandomShiftQout  = mRand->Gaus(MeanQout,RMSQout) ;
	    }
	  else
	    {
	      lRandomShiftQout  = mRand->Gaus(0.0,mDefaultResolutionQout) ;
	    }
	  Double_t lQoutSmeared = lQout + lRandomShiftQout ;
	  
	  // Qlong
	  Double_t lRandomShiftQlong = 0.0 ;
	  if(RMSQlong>0.0001 && mMomentumResolutionProvider)
	    {
	      lRandomShiftQlong  = mRand->Gaus(MeanQlong,RMSQlong) ;
	    }
	  else
	    {
	      lRandomShiftQlong  = mRand->Gaus(0.0,mDefaultResolutionQlong) ;
	    }
	  Double_t lQlongSmeared = lQlong + lRandomShiftQlong ;

	  // Fill smeared signal
	  ((TH3D*) mHistoArray->At(lKtYIndex*mNHistos+2))->Fill(lQsideSmeared,lQoutSmeared,lQlongSmeared,mWeight) ;
	  
	  // Fill smeared background
	  ((TH3D*) mHistoArray->At(lKtYIndex*mNHistos+3))->Fill(lQsideSmeared,lQoutSmeared,lQlongSmeared,1.0) ;
	}
    }
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void  T49HbtMomentumResolutionCorrection::
QCalc(T49ParticleRoot* t1,T49ParticleRoot* t2, Double_t& lQside, Double_t& lQout, Double_t& lQlong)
{
  
  // LCMS: calculate pair beta  
  Double_t lBeta = (Double_t) ((t1->GetPz() + t2->GetPz()) /(t1->GetE(mPionMass) + t2->GetE(mPionMass)));
  // Boost pairs to LCMS
  TLorentzVector P1 = T49Tool::zLorentz(t1,mPionMass,lBeta);
  TLorentzVector P2 = T49Tool::zLorentz(t2,mPionMass,lBeta);
  
  // Pair kt,qSide,qOut and qLong
  Double_t ktPAIR = TMath::Sqrt((P1.X()+P2.X())*(P1.X()+P2.X())
				+(P1.Y()+P2.Y())*(P1.Y()+P2.Y()));
  lQside  = ((P1.X()-P2.X())*(P1.Y()+P2.Y()) 
	     -(P1.Y()-P2.Y())*(P1.X()+P2.X())) / ktPAIR; // Qside              
  lQout  = (P1.X()*P1.X()-P2.X()*P2.X() 
	    +P1.Y()*P1.Y()-P2.Y()*P2.Y())     / ktPAIR; // Qout
  lQlong =  P1.Z() - P2.Z();                            // Qlong
  
  // Allow only one orientation
  if ( lQside < 0) 
    {
      lQside = -lQside; 
      lQout  = -lQout; 
      lQlong = -lQlong;
    }      

}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 TObjArray* T49HbtMomentumResolutionCorrection::Get3dCorrectionHistograms() 
{
  // Create TObjArrays with correction histos
  TObjArray* mCorrectionHistogramArray = new TObjArray(mMaxKtY) ;
  
  for(Int_t lKtYIndex = 0 ; lKtYIndex < mMaxKtY ; lKtYIndex++)
    {
      TH3D* mCorrectionHisto = (TH3D*) mHistoArray->At(lKtYIndex*mNHistos+4) ;
      mCorrectionHistogramArray->AddAt(mCorrectionHisto, lKtYIndex) ;

      // Loop over Histo and fill it
      for( Int_t mBinIndex = 0 ; 
	   mBinIndex < (mCorrectionHisto->GetNbinsX()+2)*(mCorrectionHisto->GetNbinsY()+2)*(mCorrectionHisto->GetNbinsZ()+2);
	   mBinIndex++ )
	{
	  // Per default set every bin to 1.0
	  mCorrectionHisto->SetBinContent(mBinIndex,1.0) ;
	  
	  // Now calculate the correction factor
	  Double_t mISignal = ((TH3D*) mHistoArray->At(lKtYIndex*mNHistos+0))->GetBinContent(mBinIndex) ;
	  Double_t mIBack   = ((TH3D*) mHistoArray->At(lKtYIndex*mNHistos+1))->GetBinContent(mBinIndex) ;
	  Double_t mSSignal = ((TH3D*) mHistoArray->At(lKtYIndex*mNHistos+2))->GetBinContent(mBinIndex) ;
	  Double_t mSBack   = ((TH3D*) mHistoArray->At(lKtYIndex*mNHistos+3))->GetBinContent(mBinIndex) ;
	
	  if( mISignal != 0.0 &&  mIBack != 0.0 && mSSignal != 0.0 && mSBack!=0.0 )
	    {
	      Double_t mCorrectionFactor = (mISignal/mIBack) / (mSSignal/mSBack) ;
	      if( mCorrectionFactor > 0.1 && mCorrectionFactor < 10.0 )
		{
		  mCorrectionHisto->SetBinContent(mBinIndex,mCorrectionFactor) ;
		}
	      else
		{
		  cerr << "Warning: correction factor out of bounds for bin " << mBinIndex << endl ;
		}
	    }
	 
	} // Loop over histogram bins
            
      // If we used the database print out some stats
      if(mMomentumResolutionProvider)
	{
	  cout << "Number of inquieries : " << mTotalNumberOfInquries[lKtYIndex]<< "t" ;
	  cout << "t bad inquiries : " << mNumberBadInquiries[lKtYIndex] << "  (" 
	       << (mNumberBadInquiries[lKtYIndex]/mTotalNumberOfInquries[lKtYIndex]) * 100.0<< "%)" << endl ;
	}
    } // Loop over Kt-Y bins
  
  // Return resulting TObjArray
  return mCorrectionHistogramArray ;
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


ROOT page - Class index - Class Hierarchy - Top of the page

This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.