#include "T49HbtMomentumResolutionCalculation.h"
ClassImp(T49HbtMomentumResolutionCalculator)
ClassImp(T49HbtMomentumResolutionProvider)

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// The T49HbtMomentumResolutionCalculator class is intended to produce for each qSide,qLong,qOut bin in each kt-Y bin 
// the momentum resolution in qSide,qOut,qLong and a possible shift in qSide,qOut,qLong. To obtain these values 3 vectors 
// for each cell are filled holding the difference between a Monte Carlo pair (MC) and the reconstructed pair (RC) 
// in qSide, qOut and qLong. A gaus fit to each vector gives the RMS (=resolution) and the shift, which will is needed  
// by the smear procedure in T49HbtMomentumResolutionCorrection.
//
// The constructor wants to know the dimension of the qside,qout,qlong 3d object and the kt-Y binning.
// The mixing macro calls the Use function, enertering lists of MC and the according RC tracks.
// The result can be stored in a file which will be accessed by T49HbtMomentumResolutionCorrection 
// via T49HbtMomentumResolutionProvider
//
// Structure of T49HbtMomentumResolutionCalculator:
//
//             each kt-Y cell -> ListofDisplacedMomenta[#qside * #qout * #qlong] holding TVector3s of MomDiffs
//                            -> TH3D* ResolutionMeanQside (out,long)  results
//                            -> TH3D* ResolutionRMSQside  (out,long)  results
//
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

T49HbtMomentumResolutionCalculator::T49HbtMomentumResolutionCalculator(TH3D* MomentumDifferenceOfMCTracks,
Int_t nKtBins, Double_t minKt, Double_t maxKt, Int_t nYBins, Double_t minY, Double_t maxY)
{
  // Init master 3d histo
  mMomentumDifferenceOfMCTracks = MomentumDifferenceOfMCTracks ;
  mMomentumDifferenceOfMCTracks->SetName("Master3d") ;

  // Number of cells in the master 3d histogramm 
  mNumberOfCells = (mMomentumDifferenceOfMCTracks->GetNbinsX()+2) * 
                   (mMomentumDifferenceOfMCTracks->GetNbinsY()+2) *
                   (mMomentumDifferenceOfMCTracks->GetNbinsZ()+2) ;

  // Create TObjArrays holding MC and Reconstructed tracks
  mSelectedMCTracks = new TObjArray(1000) ;
  mSelectedRETracks = new TObjArray(1000) ;

  // Fix pion mass
  mPionMass = 0.139 ;
 
  // Calculate number of kt-Y bins
  mNKtBins = nKtBins ;
  mNYBins = nYBins ;
  mMaxKtY = mNKtBins * mNYBins ;
  
  // Calculate kt-Y bin intervalls
  // Kt
  mKtBins = new Double_t[mNKtBins+1] ;
  for(Int_t index = 0 ; index < (mNKtBins+1) ; index++)
    {
      mKtBins[index] = minKt + ( ((Double_t)(maxKt-minKt)/(Double_t)mNKtBins ) * (Double_t) index) ;
    }
  // Y
  mYBins  = new Double_t[mNYBins+1] ;
  for(Int_t index = 0 ; index < (mNYBins+1) ; index++)
    {
      mYBins[index] = minY + ( (maxY-minY)/mNYBins * index) ;
    }

  // Create set of objects for each kt-Y bin
  mListOfListOfDisplacedMomenta = new TObjArray(mMaxKtY) ;
  // Number of histor RMS in qSide,qOut,qLong and shift for of them = 6
  mNumberOfResolutionHistos = 6 ; 
  mResolution = new TObjArray(mMaxKtY*mNumberOfResolutionHistos) ;
  mResolution->SetName("resultArray") ;
  
  // Loop over kt-Y bins
  for(Int_t indexKtY = 0 ; indexKtY < mMaxKtY ; indexKtY++)
    {
      // Create list of displaced momenta for each kt-Y bin
      TObjArray* lListOfListOfDisplacedMomenta = new TObjArray(mNumberOfCells) ;
      mListOfListOfDisplacedMomenta->AddAt(lListOfListOfDisplacedMomenta,indexKtY) ;

      // Create resolution 3d histos for each kt-Y bin
      for(Int_t index_MRMS_osl = 0 ; index_MRMS_osl < mNumberOfResolutionHistos ; index_MRMS_osl++)
	{
	  // Compile name of the 6 histos (Mean and RMS for qSide, qOut and qLong) for each ktY bin
	  TString name;
	  switch (index_MRMS_osl)
	    {
	    case 0:
	      name.Append("ResolutionMeanQSide") ; break ;
	    case 1:
	      name.Append("ResolutionMeanQOut") ; break ;
	    case 2:
	      name.Append("ResolutionMeanQLong") ; break ;
	    case 3:
	      name.Append("ResolutionRMSQSide") ; break ;
	    case 4:
	      name.Append("ResolutionRMSQOut") ; break ;
	    case 5:
	      name.Append("ResolutionRMSQLong") ; break ;
	    default: 
	      cerr << "Error in T49HbtMomentumResolutionCalculator Constructor, should never happen ! "<< endl ; break;
	    };
	  
	  // Add kt-Y index to histo's name 
	  name += indexKtY ;
	  
	  // Create histo as a clone from the master
	  TH3D* lResolution =  (TH3D*) mMomentumDifferenceOfMCTracks->Clone(name.Data()) ;
	  lResolution->Reset() ;
      
	  // Store pointer to histo in result TObjArray
	  mResolution->AddAt(lResolution,indexKtY*mNumberOfResolutionHistos+index_MRMS_osl) ;
	}
    }
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
T49HbtMomentumResolutionCalculator::~T49HbtMomentumResolutionCalculator()
{
  mListOfListOfDisplacedMomenta->Clear() ;
  mResolution->Clear() ;
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void T49HbtMomentumResolutionCalculator::Use(TObjArray*  MCTrackList, TObjArray*  RETrackList)
{
  // Empty TObjArrays holding MC and Reconstructed tracks
  mSelectedMCTracks->Clear() ;
  mSelectedRETracks->Clear() ;

  // Select MC tracks and find best associated reconstructed tracks  
  SelectedAssocicatedTracks(MCTrackList, RETrackList) ;
  
  // After previous step, mSelectedMCTracks and mSelectedRETracks should exist and hold same number of tracks
  if( mSelectedMCTracks->GetEntries() == mSelectedRETracks->GetEntries() && mSelectedMCTracks->GetEntries() > 0 )
    {
      // Create pairs and fill reconstructed momentum differences in appropriate TClonesArrays 
      FillDisplacedMomenta() ;
    }
  else
    {
      // This should not happen
      cerr << "Serious error: Wrong number of MC or RC pairs " << endl ;
    }
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void T49HbtMomentumResolutionCalculator::SelectedAssocicatedTracks(TObjArray*  MCTrackList, TObjArray*  RETrackList)
{
  //
  // Select reconstructed tracks which are best associated with a given MC track -> fill lists with selected tracks
  //
  Int_t MCtrackcounter = 0 ;
  for(Int_t indexMCTracks = 0; indexMCTracks < MCTrackList->GetEntries(); indexMCTracks++) 
    {
      // Get MC track
      T49ParticleMCRoot *MCTrack = (T49ParticleMCRoot *) MCTrackList->At(indexMCTracks) ;
      if(MCTrack->GetPid() == 9 && MCTrack->GetStartVertex() == 0 && MCTrack->GetNPoint()>10)
	{
	  // Search best matched primary track
	  T49ParticleRoot* BestRecTrack = NULL ;
	  Float_t BestRatio = 0.0 ;
	  Float_t BestQinv = 0.0 ;

	  // Loop over matched RC tracks
	  for (Int_t iMatch = 0 ; iMatch < MCTrack->GetNPriMatched() ; iMatch++) 
	    { 
	      // Get RC track
	      Int_t iRec = MCTrack->GetPriMatched(iMatch) ;
	      T49ParticleRoot *RecTrack = (T49ParticleRoot *) RETrackList->At(iRec) ;
	      
	      // Ratio #Matched/#PointsOnTrack
	      Float_t Ratio = (Float_t) MCTrack->GetNPriMatchPoint(iMatch)/ (Float_t) RecTrack->GetNPoint() ;
	      
	      // Qinv
	      Float_t qx  =  MCTrack->GetPx() - RecTrack->GetPx() ;
	      Float_t qy  =  MCTrack->GetPy() - RecTrack->GetPy() ;
	      Float_t qz  =  MCTrack->GetPz() - RecTrack->GetPz() ;
	      Float_t qE  =  MCTrack->GetE() - RecTrack->GetE(0.139) ;
	      Float_t Qinv = sqrt(-(qE*qE-(qx*qx+qy*qy+qz*qz))) ;
	      
	      if (iMatch == 0)
		{
		  // start values
		  BestRecTrack = RecTrack ;
		  BestRatio = Ratio ;
		  BestQinv = Qinv ;
		}
	      else 
		{
		  // check ratio and qinv to choose the only one
		  if(Ratio>BestRatio && Qinv>BestQinv)
		    {
		      BestRecTrack = RecTrack ;
		    }
		}
	    } // loop matched tracks
	     
	  // Ratio > 0.8 is required as additional criterion for a track match
	  if(BestRatio>0.8 && BestRecTrack) 
	    {
	      mSelectedMCTracks->AddAtAndExpand(MCTrack,MCtrackcounter) ;
	      mSelectedRETracks->AddAtAndExpand(BestRecTrack,MCtrackcounter) ;
	      MCtrackcounter++ ;
	    }
	} // if good mc track
    } // loop mc tracks
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void T49HbtMomentumResolutionCalculator::FillDisplacedMomenta()
{
  //
  // Construct all possible pairs of MC tracks
  //
  T49ParticleMCRoot *Part1 ;
  T49ParticleMCRoot *Part2 ;
    
  Int_t IP1 = 0 ;
  Int_t IP2 = 0 ;
  
  TIter NextParticle(mSelectedMCTracks) ;
  while ((Part1 = ((T49ParticleMCRoot *) NextParticle())))
    { 
      IP1++ ; 
      IP2 = 0 ;
      TIter NextP2(mSelectedMCTracks) ;
      while ((Part2 = ((T49ParticleMCRoot *) NextP2()))) 
	{ 
	  IP2++ ;
	  if(IP1 > IP2) 
	    {
	      FillDisplacedMomenta(Part1,Part2);
	    }
	}
    }
}
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void T49HbtMomentumResolutionCalculator::FillDisplacedMomenta(T49ParticleMCRoot* MC1, T49ParticleMCRoot* MC2)
{
  //
  // Calculate MC and RE momentum difference, add difference MC-RE to  
  // vector with index according to the MCqout,MCside,MClong position in master TH3D 
  //

  // Calculate MC momentum difference
  Double_t MCBeta = (Double_t) ((MC1->GetPz() + MC2->GetPz())/(MC1->GetE() + MC2->GetE()));
  TLorentzVector M1 = T49Tool::zLorentz(MC1,MC1->GetMass(),MCBeta) ;
  TLorentzVector M2 = T49Tool::zLorentz(MC2,MC2->GetMass(),MCBeta) ;
  
  Double_t ktPAIR = TMath::Sqrt((M1.X()+M2.X())*(M1.X()+M2.X())
				+(M1.Y()+M2.Y())*(M1.Y()+M2.Y())) ;
  Double_t MCqside = ((M1.X()-M2.X())*(M1.Y()+M2.Y())
		      -(M1.Y()-M2.Y())*(M1.X()+M2.X())) / ktPAIR ; // Qside              
  Double_t MCqout  = (M1.X()*M1.X()-M2.X()*M2.X() 
		      +M1.Y()*M1.Y()-M2.Y()*M2.Y())     / ktPAIR ; // Qout
  Double_t MCqlong =  M1.Z() - M2.Z() ;                            // Qlong

  Double_t YPAIR   = 0.5 * TMath::Log( (MC1->GetE()+MC1->GetPz()+MC2->GetE()+MC2->GetPz()) / 
				       (MC1->GetE()-MC1->GetPz()+MC2->GetE()-MC2->GetPz()) ); 
  
  Int_t mReverseOrder = 0 ;
  if(MCqside < 0.0)
    {
      MCqside = MCqside * (-1.0) ;
      MCqout  = MCqout  * (-1.0) ;
      MCqlong = MCqlong * (-1.0) ;
      mReverseOrder = 1 ;
    }
  
  // Fill only if MCqside, MCqout & MCqlong are in the range of 3d MomentumDifferenceOfMCTracks
  if( MCqside > mMomentumDifferenceOfMCTracks->GetXaxis()->GetXmin() && MCqside < mMomentumDifferenceOfMCTracks->GetXaxis()->GetXmax() &&
      MCqout > mMomentumDifferenceOfMCTracks->GetYaxis()->GetXmin() && MCqout < mMomentumDifferenceOfMCTracks->GetYaxis()->GetXmax() &&
      MCqlong > mMomentumDifferenceOfMCTracks->GetZaxis()->GetXmin() && MCqlong < mMomentumDifferenceOfMCTracks->GetZaxis()->GetXmax() )
    {
      
      // Get RC tracks
      T49ParticleRoot* RE1 = (T49ParticleRoot*) mSelectedRETracks->At(mSelectedMCTracks->IndexOf(MC1)) ; 
      T49ParticleRoot* RE2 = (T49ParticleRoot*) mSelectedRETracks->At(mSelectedMCTracks->IndexOf(MC2)) ; 
      // Calculate RC momentum difference
      Double_t REBeta = (Double_t) ((RE1->GetPz() + RE2->GetPz())/(RE1->GetE(mPionMass) + RE2->GetE(mPionMass)));
      TLorentzVector R1 = T49Tool::zLorentz(RE1,mPionMass,REBeta) ;
      TLorentzVector R2 = T49Tool::zLorentz(RE2,mPionMass,REBeta) ;
      
      Double_t REktPAIR = TMath::Sqrt((R1.X()+R2.X())*(R1.X()+R2.X())
				      +(R1.Y()+R2.Y())*(R1.Y()+R2.Y())) ;
      Double_t REqside = ((R1.X()-R2.X())*(R1.Y()+R2.Y())
			  -(R1.Y()-R2.Y())*(R1.X()+R2.X())) / REktPAIR ; // Qside              
      Double_t REqout  = (R1.X()*R1.X()-R2.X()*R2.X() 
			  +R1.Y()*R1.Y()-R2.Y()*R2.Y())     / REktPAIR ; // Qout
      Double_t REqlong =  R1.Z() - R2.Z() ;                            // Qlong
      
      if(mReverseOrder == 1) 
	{
	  REqside = REqside  * (-1.0) ;
	  REqout  = REqout   * (-1.0) ;
	  REqlong = REqlong  * (-1.0) ;
	} 

      //
      // Fill in Clonesarray at the appropriate position
      //
      
      // Get appropriate position
      Int_t GeneralIndex = mMomentumDifferenceOfMCTracks->FindBin(MCqside,MCqout,MCqlong) ;
      
      // Create TVector3 with the momentum difference
      TVector3* MomDiffReconstructed = new TVector3((REqside-MCqside),(REqout-MCqout),(REqlong-MCqlong)) ;
      
      // If we do not have a TObjArray for this position yet, create it
      Int_t lKtYIndex = GetKtYIndex(ktPAIR,YPAIR) ;
      if(lKtYIndex>=0)
	{
	  if(!((TObjArray*)mListOfListOfDisplacedMomenta->At(lKtYIndex))->At(GeneralIndex))
	    {
	      TObjArray* aNewObjArray = new TObjArray(10) ;
	      ((TObjArray*)mListOfListOfDisplacedMomenta->At(lKtYIndex))->AddAt(aNewObjArray, GeneralIndex) ;
	    }
	  
	  // Add it
	  ((TObjArray*)((TObjArray*) mListOfListOfDisplacedMomenta->At(lKtYIndex))->At(GeneralIndex))->Add(MomDiffReconstructed) ;
	}
    } // if in range of 3d histo
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
Int_t T49HbtMomentumResolutionCalculator::GetKtYIndex(Double_t lKtPair, Double_t lYPair)
{
  // Kt
  Double_t ktIndex = -1.0 ;
  for(Int_t index = 0 ; index < mNKtBins ; index++) 
    {
      if(lKtPair>mKtBins[index] && lKtPair<mKtBins[index+1])
	{
	  ktIndex = index ;
	}
    }
  // Y
  Double_t yIndex = -1.0 ;
  for(Int_t index = 0 ; index < mNYBins ; index++) 
    {
      if(lYPair>mYBins[index] && lYPair<mYBins[index+1])
	{
	  yIndex = index ;
	}
    }
  // return kty bin, if out of bonds negative value
  return ((Int_t) (ktIndex+mNKtBins*yIndex)) ;
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
TObjArray* T49HbtMomentumResolutionCalculator::GetResolution3dHisots()
{
  // Calculate and return resulting resolution histograms

  // Loop over kt-Y bins
  for(Int_t lKtYIndex = 0 ; lKtYIndex < mMaxKtY ; lKtYIndex++ )
    {
      // Loop over all cells in the master 3d histo
      for(Int_t GeneralIndex = 0 ; GeneralIndex < mNumberOfCells ; GeneralIndex++)
	{
	  // If a TObjArray exists
	  if(((TObjArray*)mListOfListOfDisplacedMomenta->At(lKtYIndex))->At(GeneralIndex))
	    {
	      // Get TObjArray of displaxed momenta
	      TObjArray* currentObjArray = (TObjArray*) ((TObjArray*) mListOfListOfDisplacedMomenta->At(lKtYIndex))->At(GeneralIndex) ;

	      // Where are we in the master 3d histo ?
	      /*
		Int_t nbinsX = mMomentumDifferenceOfMCTracks->GetNbinsX()+2 ;
		Int_t nbinsY = mMomentumDifferenceOfMCTracks->GetNbinsY()+2 ;
		Int_t nbinsZ = mMomentumDifferenceOfMCTracks->GetNbinsZ()+2 ;
		Int_t binx = GeneralIndex%nbinsX;
		Int_t biny = ((GeneralIndex-binx)/nbinsX)%nbinsY ;
		Int_t binz = (((GeneralIndex-binx)/nbinsX)-biny)/nbinsZ ;
		cout << "working on bin:  qside: " << mMomentumDifferenceOfMCTracks->GetXaxis()->GetBinCenter(binx) << "t" ;
		cout << "qout: " << mMomentumDifferenceOfMCTracks->GetYaxis()->GetBinCenter(biny) << "t" ;
		cout << "qlong: " << mMomentumDifferenceOfMCTracks->GetZaxis()->GetBinCenter(binz) << "t" << endl;
	      */

	      // If vector for this bin has more than 3 entries 
	      if(currentObjArray->GetEntries() > 3)
		{
		  // Loop over vector to calculate mean shift
		  Double_t MeanQSide     = 0.0 ;
		  Double_t MeanQOut      = 0.0 ;
		  Double_t MeanQLong     = 0.0 ;
		  Double_t MeanCounter   = 0.0 ;
		  for(Int_t VectorIndex = 0 ; VectorIndex < currentObjArray->GetEntries() ; VectorIndex++)
		    {
		      MeanCounter++ ;
		      MeanQSide += ((TVector3*) (currentObjArray->At(VectorIndex)))->Px() ;
		      MeanQOut  += ((TVector3*) (currentObjArray->At(VectorIndex)))->Py() ;
		      MeanQLong += ((TVector3*) (currentObjArray->At(VectorIndex)))->Pz() ;
		      
		      // Some output, for debugging only
		      /*
			Double_t a = (Double_t) ((TVector3*) (currentObjArray->At(VectorIndex)))->Px() ;
			TVector3* c = ((TVector3*) (currentObjArray->At(VectorIndex))) ;
			Int_t b = currentObjArray->GetEntries() ;
			cout << "qsideshift "  << ((TVector3*) (currentObjArray->At(VectorIndex)))->Px() << "t" ; 
			cout << "qoutshift  " << ((TVector3*) (currentObjArray->At(VectorIndex)))->Py() << "t";
			cout << "qlongshift " << ((TVector3*) (currentObjArray->At(VectorIndex)))->Pz()<< endl  ;
		      */
		    } 
		  MeanQSide = MeanQSide / MeanCounter ;
		  MeanQOut  = MeanQOut  / MeanCounter ;
		  MeanQLong = MeanQLong / MeanCounter ;
		  
		  // Fill mean shift into 3d histo
		  ((TH3D*) mResolution->At(lKtYIndex*6))->SetBinContent(GeneralIndex,MeanQSide) ;
		  ((TH3D*) mResolution->At(lKtYIndex*6+1))->SetBinContent(GeneralIndex,MeanQOut) ;
		  ((TH3D*) mResolution->At(lKtYIndex*6+2))->SetBinContent(GeneralIndex,MeanQLong) ;
		  		  
		  // Loop over vector to calculate RMS
		  Double_t RMSQSide     = 0.0 ;
		  Double_t RMSQOut      = 0.0 ;
		  Double_t RMSQLong     = 0.0 ;
		  Double_t RMSCounter   = 0.0 ;
		  for(Int_t VectorIndex = 0 ; VectorIndex < currentObjArray->GetEntries() ; VectorIndex++)
		    {
		      RMSQSide += pow((((TVector3*) (currentObjArray->At(VectorIndex)))->Px() -  MeanQSide),2) ;
		      RMSQOut  += pow((((TVector3*) (currentObjArray->At(VectorIndex)))->Py() -  MeanQOut),2)  ; 
		      RMSQLong += pow((((TVector3*) (currentObjArray->At(VectorIndex)))->Pz() -  MeanQLong),2) ;
		      RMSCounter++ ;
		      // Some output, for debugging only
		      /*
			Double_t a = (Double_t) ((TVector3*) (currentObjArray->At(VectorIndex)))->Px() ;
			TVector3* c = ((TVector3*) (currentObjArray->At(VectorIndex))) ;
			Int_t b = currentObjArray->GetEntries() ;
			cout << "qsideRMS "  << pow((((TVector3*) (currentObjArray->At(VectorIndex)))->Px() -  MeanQSide),2)   << "t" ; 
			cout << "qoutRMS  " << pow((((TVector3*) (currentObjArray->At(VectorIndex)))->Py() -  MeanQOut),2)  << "t";
			cout << "qlongRMS " << pow((((TVector3*) (currentObjArray->At(VectorIndex)))->Pz() -  MeanQLong),2)  << endl  ;
		      */
		} 
		  RMSQSide = sqrt(RMSQSide / RMSCounter) ;
		  RMSQOut  = sqrt(RMSQOut  / RMSCounter) ;
		  RMSQLong = sqrt(RMSQLong / RMSCounter) ;
		  
		  // Fill mean into 3d histo
		  ((TH3D*) mResolution->At(lKtYIndex*6+3))->SetBinContent(GeneralIndex,RMSQSide) ;
		  ((TH3D*) mResolution->At(lKtYIndex*6+4))->SetBinContent(GeneralIndex,RMSQOut) ;
		  ((TH3D*) mResolution->At(lKtYIndex*6+5))->SetBinContent(GeneralIndex,RMSQLong) ;
		  
		} // if Entries > 0
	    } // if TObjArray
	} // loop over 3d master histo
    } // loop over kty index
  // Return ObjArray with the Result 3d histos
  return mResolution ;
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void T49HbtMomentumResolutionCalculator::PrintResolution3dHisots(Int_t lKtYIndex)
{
  // This member function plots for a given kty bin the mean and rms of each cell 
  // along with the cells position in the master TH3D.
  // If commented out, also the vector is plotted from which rms and mean are derived

  // Make sure the resolution histos are filled
  this->GetResolution3dHisots() ;
  
  // Loop over all cells in the master 3d histo
  for(Int_t GeneralIndex = 0 ; GeneralIndex < mNumberOfCells ; GeneralIndex++)
    {
      // If a TObjArray exists
      if(((TObjArray*) mListOfListOfDisplacedMomenta->At(lKtYIndex))->At(GeneralIndex))
	{
	  TObjArray* currentObjArray = (TObjArray*) ((TObjArray*)mListOfListOfDisplacedMomenta->At(lKtYIndex))->At(GeneralIndex) ;
	  if(currentObjArray->GetEntries()>3)
	    {
	      // Find the qout,qside,qlong value in the master 3d histo
	      Int_t nbinsX = mMomentumDifferenceOfMCTracks->GetNbinsX()+2 ;
	      Int_t nbinsY = mMomentumDifferenceOfMCTracks->GetNbinsY()+2 ;
	      Int_t nbinsZ = mMomentumDifferenceOfMCTracks->GetNbinsZ()+2 ;
	      Int_t binx = GeneralIndex%nbinsX;
	      Int_t biny = ((GeneralIndex-binx)/nbinsX)%nbinsY ;
	      Int_t binz = (((GeneralIndex-binx)/nbinsX)-biny)/nbinsZ ;
	      
	      cout << "General index : " << "t" <<  GeneralIndex 
		   << "t must equal t" << mMomentumDifferenceOfMCTracks->GetBin(binx,biny,binz) << "t" ;
	      cout << "t bin x,y,z: t " << binx  << "t" << biny  << "t" << binz  << "t"  << endl;
	      cout << "qside: " << mMomentumDifferenceOfMCTracks->GetXaxis()->GetBinCenter(binx) << "t" ;
	      cout << "qout: " << mMomentumDifferenceOfMCTracks->GetYaxis()->GetBinCenter(biny) << "t" ;
	      cout << "qlong: " << mMomentumDifferenceOfMCTracks->GetZaxis()->GetBinCenter(binz) << "t" << endl;
	     	     
	      cout << "Mean qside : " << ((TH3D*) mResolution->At(lKtYIndex*6))->GetBinContent(GeneralIndex) << "t"; 
	      cout << "Mean qOut : "  << ((TH3D*) mResolution->At(lKtYIndex*6+1))->GetBinContent(GeneralIndex) << "t"; 
	      cout << "Mean qLong : " << ((TH3D*) mResolution->At(lKtYIndex*6+2))->GetBinContent(GeneralIndex) << "t"; 
	      cout << "RMS qside : "  << ((TH3D*) mResolution->At(lKtYIndex*6+3))->GetBinContent(GeneralIndex) << "t";
	      cout << "RMS qOut : "   << ((TH3D*) mResolution->At(lKtYIndex*6+4))->GetBinContent(GeneralIndex) << "t";
	      cout << "RMS qLong : "  << ((TH3D*) mResolution->At(lKtYIndex*6+5))->GetBinContent(GeneralIndex) << "t" << endl ;
	     
	      // Show vector, which was used to calculate mean and rms
	      cout << "Number of entries : " << currentObjArray->GetEntries() << endl ;
	      for(Int_t VectorIndex = 0 ; VectorIndex < currentObjArray->GetEntries() ; VectorIndex++)
		{
		  //Double_t a = (Double_t) ((TVector3*) (currentObjArray->At(VectorIndex)))->Px() ;
		  //TVector3* c = ((TVector3*) (currentObjArray->At(VectorIndex))) ;
		  //Int_t b = currentObjArray->GetEntries() ;
		  cout << "t qside "  << ((TVector3*) (currentObjArray->At(VectorIndex)))->Px() << "t" ; 
		  cout << "qout  " << ((TVector3*) (currentObjArray->At(VectorIndex)))->Py() << "t";
		  cout <<"qlong " << ((TVector3*) (currentObjArray->At(VectorIndex)))->Pz()<< endl  ;
		}
	    }
	}
    }
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
TObjArray* T49HbtMomentumResolutionCalculator::GetResolutionDistributions(Int_t lKtYIndex)
{
  // This member function returns 6 TH1D's filled with the mean and the rms values for qSide, qOut & qLong of all cells 
  // of a given kty bin
  
  // Make sure the resolution histos are filled
  this->GetResolution3dHisots() ;

  // Create object array holding 6 distribution histograms
  TObjArray* mResolutionDistribution = new TObjArray(6) ;
  TH1D* MeanRsideDistribution = new TH1D("MeanRsideDistribution","MeanRsideDistribution",1000,-0.05,0.05) ;
  mResolutionDistribution->AddAt(MeanRsideDistribution,0) ;
  TH1D* MeanRoutDistribution = new TH1D("MeanRoutDistribution","MeanRoutDistribution",1000,-0.05,0.05) ;
  mResolutionDistribution->AddAt(MeanRoutDistribution,1) ;
  TH1D* MeanRlongDistribution = new TH1D("MeanRlongDistribution","MeanRlongDistribution",1000,-0.05,0.05) ;
  mResolutionDistribution->AddAt(MeanRlongDistribution,2) ;
  TH1D* RMSRsideDistribution = new TH1D("RMSRsideDistribution","RMSRsideDistribution",1000,0.0,0.05) ;
  mResolutionDistribution->AddAt(RMSRsideDistribution,3) ;
  TH1D* RMSRoutDistribution = new TH1D("RMSRoutDistribution","RMSRoutDistribution",1000,0.0,0.05) ;
  mResolutionDistribution->AddAt(RMSRoutDistribution,4) ;
  TH1D* RMSRlongDistribution = new TH1D("RMSRlongDistribution","RMSRlongDistribution",1000,0.0,0.05) ;
  mResolutionDistribution->AddAt(RMSRlongDistribution,5) ;
  
  // Loop over all cells in the master 3d histo and fill the 6 distribution histograms
  for(Int_t GeneralIndex = 0 ; GeneralIndex < mNumberOfCells ; GeneralIndex++)
    {
      MeanRsideDistribution->Fill(((TH3D*) mResolution->At(lKtYIndex*6))->GetBinContent(GeneralIndex)) ;
      MeanRoutDistribution->Fill(((TH3D*) mResolution->At(lKtYIndex*6+1))->GetBinContent(GeneralIndex)) ;
      MeanRlongDistribution->Fill(((TH3D*) mResolution->At(lKtYIndex*6+2))->GetBinContent(GeneralIndex)) ;
      
      RMSRsideDistribution->Fill(((TH3D*) mResolution->At(lKtYIndex*6+3))->GetBinContent(GeneralIndex)) ;
      RMSRoutDistribution->Fill(((TH3D*) mResolution->At(lKtYIndex*6+4))->GetBinContent(GeneralIndex)) ;
      RMSRlongDistribution->Fill(((TH3D*) mResolution->At(lKtYIndex*6+5))->GetBinContent(GeneralIndex)) ;
    }

  // Return Obj array
  return mResolutionDistribution ;
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void T49HbtMomentumResolutionCalculator::SaveResolution(Char_t* Name)
{
  //
  // Store the calculated resolution TH3Ds in a root file, 
  // to have easy access using class T49HbtMomentumResolutionProvider
  //

  // Create TGraphs holding information about kt-Y bins
  Double_t lDummy[mNKtBins+mNYBins] ;
  TGraph* lKtBins = new TGraph(mNKtBins+1,mKtBins,lDummy) ;
  lKtBins->SetName("KtBins");
  TGraph* lYBins =  new TGraph(mNYBins+1,mYBins,lDummy) ;
  lYBins->SetName("YBins");

  // Open resuls to file
  TFile mResultFile(Name,"recreate") ;
  lKtBins->Write() ;
  lYBins->Write() ;
  mMomentumDifferenceOfMCTracks->Write() ;
  mResolution->Write("resultArray",1) ;
  mResultFile.Close() ;
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 T49HbtMomentumResolutionProvider::T49HbtMomentumResolutionProvider(Char_t* Name)
{
  // Open result file
  TFile* mResultFile = new TFile(Name,"read") ;
  
  // Get number and size of kt-Y bins 
  // Kt
  TGraph* lgKtBins = (TGraph*) mResultFile->Get("KtBins") ;
  mNKtBins = lgKtBins->GetN()-1 ;
  mKtBins = new Double_t[mNKtBins+1] ; 
  Double_t* lKtBins = lgKtBins->GetX() ;
  for(Int_t index = 0 ; index < mNKtBins+1 ; index++)  mKtBins[index] = lKtBins[index] ;
  // Y
  TGraph* lgYBins = (TGraph*) mResultFile->Get("YBins") ;
  mNYBins = lgYBins->GetN()-1 ;
  mYBins = new Double_t[mNYBins+1] ; 
  Double_t* lYBins = lgYBins->GetX() ;
  for(Int_t index = 0 ; index < mNYBins+1 ; index++)  mYBins[index] = lYBins[index] ;
 
  // Get TObjArray with resolution histos 
  mResolution = (TObjArray*) mResultFile->Get("resultArray") ;

  // Get master 3d Histo
  mMaster3d = (TH3D*) mResultFile->Get("Master3d") ;
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 Int_t T49HbtMomentumResolutionProvider::GetKtYIndex(Double_t lKtPair, Double_t lYPair)
{
  // Kt
  Double_t ktIndex = -1.0 ;
  for(Int_t index = 0 ; index < mNKtBins ; index++) 
    {
      if(lKtPair>=mKtBins[index] && lKtPair<mKtBins[index+1])
	{
	  ktIndex = index ;
	}
    }
  // Y
  Double_t yIndex = -1.0 ;
  for(Int_t index = 0 ; index < mNYBins ; index++) 
    {
      if(lYPair>=mYBins[index] && lYPair<mYBins[index+1])
	{
	  yIndex = index ;
	}
    }
  // return kty bin, if out of bonds negative value
  if((ktIndex+mNKtBins*yIndex)>24)
    {
      ktIndex = -1.0 ;
      cerr <<"This should not happen !" << endl ; 
    }

  return ((Int_t) (ktIndex+mNKtBins*yIndex)) ;

}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 void T49HbtMomentumResolutionProvider::GetMeanAndRMS(Double_t ktPair, Double_t YPair, 
						     Double_t qside, Double_t qout, Double_t qlong,
						     Double_t& MeanQside, Double_t& MeanQout, Double_t& MeanQlong,
						     Double_t& RMSQside , Double_t& RMSQout , Double_t& RMSQlong ) 
  
{
  // Get index of this qout-qside-qlong bin
  Int_t lGeneralIndex = mMaster3d->FindBin(qside,qout,qlong) ;
  // Get kt-Y index of the pair
  Int_t lKtYIndex = this->GetKtYIndex(ktPair,YPair) ;
  if(lKtYIndex>=0)
    {
      // Fill the requested Mean and RMS values
      MeanQside = ((TH3D*) mResolution->At(lKtYIndex*6))->GetBinContent(lGeneralIndex) ;
      MeanQout  = ((TH3D*) mResolution->At(lKtYIndex*6+1))->GetBinContent(lGeneralIndex) ;
      MeanQlong = ((TH3D*) mResolution->At(lKtYIndex*6+2))->GetBinContent(lGeneralIndex) ;
      RMSQside  = ((TH3D*) mResolution->At(lKtYIndex*6+3))->GetBinContent(lGeneralIndex) ;
      RMSQout   = ((TH3D*) mResolution->At(lKtYIndex*6+4))->GetBinContent(lGeneralIndex) ;
      RMSQlong  = ((TH3D*) mResolution->At(lKtYIndex*6+5))->GetBinContent(lGeneralIndex) ;
      if(RMSQside<0 || RMSQout<0 || RMSQlong<0)
	{
	  cerr << "This should never happen ! "<< endl ;
	}
    }
  else
    {
      MeanQside = 0.0 ;
      MeanQout  = 0.0 ;
      MeanQlong = 0.0 ;
      RMSQside  = 0.0 ;
      RMSQout   = 0.0 ;
      RMSQlong  = 0.0 ;
    }
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////



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.