[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemtoUser / AliFemtoCorrFctn3DSphericalEMCIC.cxx

/***************************************************************************
 *
 * $Id: AliFemtoCorrFctn3DSphercicalEMCIC.h  $
 *
 * Author: Nicolas Bock, Ohio State University, bock@mps.ohio-state.edu
 ***************************************************************************
 *
 * Description: Calculates of the 3D Correlation Function in Spherical
 *              coordinates, and also produces histograms to calculate 
 *              Energy Momentum Conservation Induced Correlations  (EMCICs)
 *
 * This Class produces the following histograms as function of Q, theta, phi
 * (for both real and mixed pairs):
 *        1)   E1 + E2
 *        2)   E1 * E2
 *        3)   Pt1*Pt2
 *        4)   Pz1*Pz2
 *  
 * This class is similar to AliFemtoCorrFctn3DSpherical, but it uses Q
 * instead of K to do the binning. 
 * 
 * NOTE: The EMCIC histograms are not averaged in this class, to obtain 
 * the average, the user needs to divide the real pair histograms by 
 * the numerator, and the mixed pair histograms by the denominator.
 *
 ***************************************************************************
 *
 **************************************************************************/

#include "AliFemtoCorrFctn3DSphericalEMCIC.h"
#include <TMath.h>
#include <TVector2.h>
#include <cstdio>

#ifdef __ROOT__ 
ClassImp(AliFemtoCorrFctn3DSphericalEMCIC)
#endif


//____________________________
AliFemtoCorrFctn3DSphericalEMCIC::AliFemtoCorrFctn3DSphericalEMCIC(char* title, const int& nqbins, const float& QLo, const float& QHi, const int& nphibins, const int& ncthetabins):
  AliFemtoCorrFctn(),
  fNumerator(0),
  fDenominator(0),
  fEnergyTotalReal(0),  
  fEnergyMultReal(0),        
  fPzMultReal(0),      
  fPtMultReal(0),            
  fEnergyTotalMix (0),      
  fEnergyMultMix (0),      
  fPzMultMix(0),            
  fPtMultMix(0),
  fPairCut(0x0)
{
  
  // To better sample on phi shift bin low edge by binsize/2 = Pi/numBins.
  Double_t shiftPhi=TMath::Pi()/nphibins;
  // set up numerator
  char tTitNum[100] = "Num";
  strcat(tTitNum,title);
  fNumerator = new TH3D(tTitNum,title,nqbins,QLo,QHi,nphibins,
			-TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,ncthetabins,-1.0,1.0);
  // set up denominator
  char tTitDen[100] = "Den";
  strcat(tTitDen,title);
  fDenominator = new TH3D(tTitDen,title,nqbins,QLo,QHi,nphibins,
			  -TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,ncthetabins,-1.0,1.0);

  //Added histograms to calculate EMCICs , Nicolas Bock 19.01.2010
  //Setup EnergyTotalReal
  char tTitNum1[100] = "ESumReal";
  strcat(tTitNum1,title);
  fEnergyTotalReal = new TH3D(tTitNum1,title,nqbins,QLo,QHi,nphibins,
			      -TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,ncthetabins,-1.0,1.0);
 
  //Setup EnergyMultReal
  char tTitNum2[100] = "EMultReal";
  strcat(tTitNum2,title);
  fEnergyMultReal = new TH3D(tTitNum2,title,nqbins,QLo,QHi,nphibins,
			     -TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,ncthetabins,-1.0,1.0);
  
  //Setup Pz MultReal
  char tTitNum3[100] = "PzMultReal";
  strcat(tTitNum3,title);
  fPzMultReal = new TH3D(tTitNum3,title,nqbins,QLo,QHi,nphibins,
			 -TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,ncthetabins,-1.0,1.0);

  //Setup Pt MultReal
  char tTitNum4[100] = "PtMultReal";
  strcat(tTitNum4,title);
  fPtMultReal = new TH3D(tTitNum4,title,nqbins,QLo,QHi,nphibins,
			 -TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,ncthetabins,-1.0,1.0);
  

  //Setup EnergyTotalMix
  char tTitNum5[100] = "ESumMix";
  strcat(tTitNum5,title);
  fEnergyTotalMix = new TH3D(tTitNum5,title,nqbins,QLo,QHi,nphibins,
			     -TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,ncthetabins,-1.0,1.0);
  
  //Setup EnergyMultMix
  char tTitNum6[100] = "EMultMix";
  strcat(tTitNum6,title);
  fEnergyMultMix = new TH3D(tTitNum6,title,nqbins,QLo,QHi,nphibins,
			    -TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,ncthetabins,-1.0,1.0);
  
  //Setup Pz MultMix
  char tTitNum7[100] = "PzMultMix";
  strcat(tTitNum7,title);
  fPzMultMix = new TH3D(tTitNum7,title,nqbins,QLo,QHi,nphibins,
			-TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,ncthetabins,-1.0,1.0);

  //Setup Pt MultMix
  char tTitNum8[100] = "PtMultMix";
  strcat(tTitNum8,title);
  fPtMultMix = new TH3D(tTitNum8,title,nqbins,QLo,QHi,nphibins,
			-TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,
			ncthetabins,-1.0,1.0);

  // to enable error bar calculation...
  fNumerator->Sumw2();
  fDenominator->Sumw2();
  fEnergyTotalReal->Sumw2();
  fEnergyMultReal->Sumw2();
  fPzMultReal->Sumw2();
  fPtMultReal->Sumw2();
  fEnergyTotalMix->Sumw2();
  fEnergyMultMix->Sumw2();
  fPzMultMix->Sumw2();
  fPtMultMix->Sumw2();


}

AliFemtoCorrFctn3DSphericalEMCIC::AliFemtoCorrFctn3DSphericalEMCIC(const AliFemtoCorrFctn3DSphericalEMCIC& aCorrFctn) :
  AliFemtoCorrFctn(),
  fNumerator(0),
  fDenominator(0),
  fEnergyTotalReal(0),
  fEnergyMultReal(0),        
  fPzMultReal(0),      
  fPtMultReal(0),            
  fEnergyTotalMix (0),      
  fEnergyMultMix (0),      
  fPzMultMix(0),            
  fPtMultMix(0),
  fPairCut(0x0)
{
  // Copy constructor
  fNumerator = new TH3D(*aCorrFctn.fNumerator);
  fDenominator = new TH3D(*aCorrFctn.fDenominator);
  fEnergyTotalReal = new TH3D(*aCorrFctn.fEnergyTotalReal);
  fEnergyMultReal = new TH3D(*aCorrFctn.fEnergyMultReal);
  fPzMultReal = new TH3D(*aCorrFctn.fPzMultReal);
  fPtMultReal = new TH3D(*aCorrFctn.fPtMultReal);
  fEnergyTotalMix = new TH3D(*aCorrFctn.fEnergyTotalMix);
  fEnergyMultMix = new TH3D(*aCorrFctn.fEnergyMultMix);
  fPzMultMix = new TH3D(*aCorrFctn.fPzMultMix);
  fPtMultMix = new TH3D(*aCorrFctn.fPtMultMix);
  fPairCut = aCorrFctn.fPairCut;
}
//____________________________
AliFemtoCorrFctn3DSphericalEMCIC::~AliFemtoCorrFctn3DSphericalEMCIC(){
  // Destructor
  delete fNumerator;
  delete fDenominator;
  delete fEnergyTotalReal;
  delete fEnergyMultReal;        
  delete fPzMultReal;     
  delete fPtMultReal;            
  delete fEnergyTotalMix;      
  delete fEnergyMultMix; 
  delete fPzMultMix;   
  delete fPtMultMix;
}
//_________________________
AliFemtoCorrFctn3DSphericalEMCIC& AliFemtoCorrFctn3DSphericalEMCIC::operator=(const AliFemtoCorrFctn3DSphericalEMCIC& aCorrFctn)
{
  // assignment operator
  if (this == &aCorrFctn)
    return *this;

  if (fNumerator) delete fNumerator;
  fNumerator = new TH3D(*aCorrFctn.fNumerator);
  if (fDenominator) delete fDenominator;
  fDenominator = new TH3D(*aCorrFctn.fDenominator);
  if (fEnergyTotalReal) delete fEnergyTotalReal;
  fEnergyTotalReal = new TH3D(*aCorrFctn.fEnergyTotalReal);
  if (fEnergyMultReal) delete fEnergyMultReal;
  fEnergyMultReal = new TH3D(*aCorrFctn.fEnergyMultReal);
  if (fPzMultReal) delete fPzMultReal;
  fPzMultReal = new TH3D(*aCorrFctn.fPzMultReal);
  if (fPtMultReal) delete fPtMultReal;
  fPtMultReal = new TH3D(*aCorrFctn.fPtMultReal);
  if (fEnergyTotalMix) delete fEnergyTotalMix;
  fEnergyTotalMix = new TH3D(*aCorrFctn.fEnergyTotalMix);
  if (fEnergyMultMix) delete fEnergyMultMix;
  fEnergyMultMix = new TH3D(*aCorrFctn.fEnergyMultMix);
  if (fPzMultMix) delete fPzMultMix;
  fPzMultMix = new TH3D(*aCorrFctn.fPzMultMix);
  if (fPtMultMix) delete fPtMultMix;
  fPtMultMix = new TH3D(*aCorrFctn.fPtMultMix);
  fPairCut = aCorrFctn.fPairCut;
  
  return *this;
}

//_________________________
void AliFemtoCorrFctn3DSphericalEMCIC::WriteOutHistos(){
  // Write out all histograms to file
  fNumerator->Write();
  fDenominator->Write();
  fEnergyTotalReal->Write();
  fEnergyMultReal->Write();        
  fPzMultReal->Write();      
  fPtMultReal->Write();            
  fEnergyTotalMix->Write();      
  fEnergyMultMix->Write();      
  fPzMultMix->Write();            
  fPtMultMix->Write();
}
//______________________________
TList* AliFemtoCorrFctn3DSphericalEMCIC::GetOutputList()
{
  // Prepare the list of objects to be written to the output
  TList *tOutputList = new TList();

  tOutputList->Add(fNumerator); 
  tOutputList->Add(fDenominator);  
  tOutputList->Add(fEnergyTotalReal);
  tOutputList->Add(fEnergyMultReal);        
  tOutputList->Add(fPzMultReal);      
  tOutputList->Add(fPtMultReal);            
  tOutputList->Add(fEnergyTotalMix );      
  tOutputList->Add(fEnergyMultMix );      
  tOutputList->Add(fPzMultMix);            
  tOutputList->Add(fPtMultMix);
  return tOutputList;
}

//_________________________
void AliFemtoCorrFctn3DSphericalEMCIC::Finish(){
  // here is where we should normalize, fit, etc...
}

//____________________________
AliFemtoString AliFemtoCorrFctn3DSphericalEMCIC::Report(){
  // Construct the report
  string stemp = "PRF Frame SphericalEMCIC 3D Correlation Function Report:\n";
  char ctemp[100];
  sprintf(ctemp,"Number of entries in numerator:\t%E\n",fNumerator->GetEntries());
  stemp += ctemp;
  sprintf(ctemp,"Number of entries in denominator:\t%E\n",fDenominator->GetEntries());
  stemp += ctemp;

  if (fPairCut){
    sprintf(ctemp,"Here is the PairCut specific to this CorrFctn\n");
    stemp += ctemp;
    stemp += fPairCut->Report();
  }
  else{
    sprintf(ctemp,"No PairCut specific to this CorrFctn\n");
    stemp += ctemp;
  }

  //  
  AliFemtoString returnThis = stemp;
  return returnThis;
}
//____________________________
void AliFemtoCorrFctn3DSphericalEMCIC::AddRealPair( AliFemtoPair* pair){
  // perform operations on real pairs
  if (fPairCut){
    if (!(fPairCut->Pass(pair))) return;
  }

  //                          //Changed K to Q to be in LCMS, N. Bock
  double tKO = pair->QOutCMS();  //KOut();
  double tKS = pair->QSideCMS();  //KSide(); 
  double tKL = pair->QLongCMS();  //  KLong();

  double tKR = sqrt(tKO*tKO + tKS*tKS + tKL*tKL);
  double tKC;
  if ( fabs(tKR) < 1e-10 ) tKC = 0.0;
  else tKC=tKL/tKR;
  double tKP=atan2(tKS,tKO);

  fNumerator->Fill(tKR,tKP,tKC);
  
  // EMCICs  
  AliFemtoLorentzVector tMom1 = pair->Track1()->FourMomentum();
  AliFemtoLorentzVector tMom2 = pair->Track2()->FourMomentum();
  double tE1 = tMom1.e();
  double tE2 = tMom2.e();
  double tPz1 = tMom1.pz();
  double tPz2 = tMom2.pz();
  TVector2 tPt1;
  TVector2 tPt2;
  cout <<"AddRealPair" << "E1,E2, K,Ks,KL" <<  
    tE1 <<" " << tE2<<" " <<  tKR <<" " <<  tKS <<" " << tKL<<endl;
  tPt1.Set(tMom1.px(),tMom1.py());
  tPt2.Set(tMom2.px(),tMom2.py());
  double tPt1DotPt2=tPt1*tPt2;
  
  fEnergyTotalReal->Fill(tKR,tKP,tKC,tE1+tE2);
  fEnergyMultReal->Fill(tKR,tKP,tKC,tE1*tE2);
  fPzMultReal->Fill(tKR,tKP,tKC,tPz1*tPz2);
  fPtMultReal->Fill(tKR,tKP,tKC,tPt1DotPt2);
   

}
//____________________________
void AliFemtoCorrFctn3DSphericalEMCIC::AddMixedPair( AliFemtoPair* pair){
  // perform operations on mixed pairs
  if (fPairCut){
    if (!(fPairCut->Pass(pair))) return;
  }
 //                          //Changed K to Q to be in LCMS, N. Bock
  double tKO = pair->QOutCMS();  //KOut();
  double tKS = pair->QSideCMS();  //KSide(); 
  double tKL = pair->QLongCMS();  //  KLong();

  double tKR = sqrt(tKO*tKO + tKS*tKS + tKL*tKL);
  double tKC;
  if ( fabs(tKR) < 1e-10 ) tKC = 0.0;
  else tKC=tKL/tKR;
  double tKP=atan2(tKS,tKO);

  fDenominator->Fill(tKR,tKP,tKC);

  // EMCICs   
  AliFemtoLorentzVector tMom1 = pair->Track1()->FourMomentum();
  AliFemtoLorentzVector tMom2 = pair->Track2()->FourMomentum();
  double tE1 = tMom1.e();
  double tE2 = tMom2.e();
  double tPz1 = tMom1.pz();
  double tPz2 = tMom2.pz();
  TVector2 tPt1;
  TVector2 tPt2;
  tPt1.Set(tMom1.px(),tMom1.py());
  tPt2.Set(tMom2.px(),tMom2.py());
  double tPt1DotPt2=tPt1*tPt2;
  
  fEnergyTotalMix->Fill(tKR,tKP,tKC,tE1+tE2);
  fEnergyMultMix->Fill(tKR,tKP,tKC,tE1*tE2);
  fPzMultMix->Fill(tKR,tKP,tKC,tPz1*tPz2);
  fPtMultMix->Fill(tKR,tKP,tKC,tPt1DotPt2);
  
}
Commit	Line	Data
fa35b516	1	/***************************************************************************
	2	*
	3	* $Id: AliFemtoCorrFctn3DSphercicalEMCIC.h $
	4	*
	5	* Author: Nicolas Bock, Ohio State University, bock@mps.ohio-state.edu
	6	***************************************************************************
	7	*
	8	* Description: Calculates of the 3D Correlation Function in Spherical
	9	* coordinates, and also produces histograms to calculate
	10	* Energy Momentum Conservation Induced Correlations (EMCICs)
	11	*
	12	* This Class produces the following histograms as function of Q, theta, phi
	13	* (for both real and mixed pairs):
	14	* 1) E1 + E2
	15	* 2) E1 * E2
	16	* 3) Pt1*Pt2
	17	* 4) Pz1*Pz2
	18	*
	19	* This class is similar to AliFemtoCorrFctn3DSpherical, but it uses Q
	20	* instead of K to do the binning.
	21	*
	22	* NOTE: The EMCIC histograms are not averaged in this class, to obtain
	23	* the average, the user needs to divide the real pair histograms by
	24	* the numerator, and the mixed pair histograms by the denominator.
	25	*
	26	***************************************************************************
	27	*
	28	**************************************************************************/
	29
	30	#include "AliFemtoCorrFctn3DSphericalEMCIC.h"
	31	#include <TMath.h>
	32	#include <TVector2.h>
	33	#include <cstdio>
	34
	35	#ifdef __ROOT__
	36	ClassImp(AliFemtoCorrFctn3DSphericalEMCIC)
	37	#endif
	38
	39
	40	//____________________________
	41	AliFemtoCorrFctn3DSphericalEMCIC::AliFemtoCorrFctn3DSphericalEMCIC(char* title, const int& nqbins, const float& QLo, const float& QHi, const int& nphibins, const int& ncthetabins):
	42	AliFemtoCorrFctn(),
	43	fNumerator(0),
	44	fDenominator(0),
	45	fEnergyTotalReal(0),
	46	fEnergyMultReal(0),
	47	fPzMultReal(0),
	48	fPtMultReal(0),
	49	fEnergyTotalMix (0),
	50	fEnergyMultMix (0),
	51	fPzMultMix(0),
	52	fPtMultMix(0),
	53	fPairCut(0x0)
	54	{
	55
	56	// To better sample on phi shift bin low edge by binsize/2 = Pi/numBins.
	57	Double_t shiftPhi=TMath::Pi()/nphibins;
	58	// set up numerator
	59	char tTitNum[100] = "Num";
	60	strcat(tTitNum,title);
	61	fNumerator = new TH3D(tTitNum,title,nqbins,QLo,QHi,nphibins,
	62	-TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,ncthetabins,-1.0,1.0);
	63	// set up denominator
	64	char tTitDen[100] = "Den";
65	strcat(tTitDen,title);
66	fDenominator = new TH3D(tTitDen,title,nqbins,QLo,QHi,nphibins,
67	-TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,ncthetabins,-1.0,1.0);
68
69	//Added histograms to calculate EMCICs , Nicolas Bock 19.01.2010
70	//Setup EnergyTotalReal
71	char tTitNum1[100] = "ESumReal";
72	strcat(tTitNum1,title);
73	fEnergyTotalReal = new TH3D(tTitNum1,title,nqbins,QLo,QHi,nphibins,
74	-TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,ncthetabins,-1.0,1.0);
75
76	//Setup EnergyMultReal
77	char tTitNum2[100] = "EMultReal";
78	strcat(tTitNum2,title);
79	fEnergyMultReal = new TH3D(tTitNum2,title,nqbins,QLo,QHi,nphibins,
80	-TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,ncthetabins,-1.0,1.0);
81
82	//Setup Pz MultReal
83	char tTitNum3[100] = "PzMultReal";
84	strcat(tTitNum3,title);
85	fPzMultReal = new TH3D(tTitNum3,title,nqbins,QLo,QHi,nphibins,
86	-TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,ncthetabins,-1.0,1.0);
87
88	//Setup Pt MultReal
89	char tTitNum4[100] = "PtMultReal";
90	strcat(tTitNum4,title);
91	fPtMultReal = new TH3D(tTitNum4,title,nqbins,QLo,QHi,nphibins,
92	-TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,ncthetabins,-1.0,1.0);
93
94
95
96	//Setup EnergyTotalMix
97	char tTitNum5[100] = "ESumMix";
98	strcat(tTitNum5,title);
99	fEnergyTotalMix = new TH3D(tTitNum5,title,nqbins,QLo,QHi,nphibins,
100	-TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,ncthetabins,-1.0,1.0);
101
102	//Setup EnergyMultMix
103	char tTitNum6[100] = "EMultMix";
104	strcat(tTitNum6,title);
105	fEnergyMultMix = new TH3D(tTitNum6,title,nqbins,QLo,QHi,nphibins,
106	-TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,ncthetabins,-1.0,1.0);
107
108	//Setup Pz MultMix
109	char tTitNum7[100] = "PzMultMix";
110	strcat(tTitNum7,title);
111	fPzMultMix = new TH3D(tTitNum7,title,nqbins,QLo,QHi,nphibins,
112	-TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,ncthetabins,-1.0,1.0);
113
114	//Setup Pt MultMix
115	char tTitNum8[100] = "PtMultMix";
116	strcat(tTitNum8,title);
117	fPtMultMix = new TH3D(tTitNum8,title,nqbins,QLo,QHi,nphibins,
118	-TMath::Pi()-shiftPhi,TMath::Pi()-shiftPhi,
119	ncthetabins,-1.0,1.0);
120
121	// to enable error bar calculation...
122	fNumerator->Sumw2();
123	fDenominator->Sumw2();
124	fEnergyTotalReal->Sumw2();
125	fEnergyMultReal->Sumw2();
126	fPzMultReal->Sumw2();
127	fPtMultReal->Sumw2();
128	fEnergyTotalMix->Sumw2();
129	fEnergyMultMix->Sumw2();
130	fPzMultMix->Sumw2();
131	fPtMultMix->Sumw2();
132
133
134	}
135
136	AliFemtoCorrFctn3DSphericalEMCIC::AliFemtoCorrFctn3DSphericalEMCIC(const AliFemtoCorrFctn3DSphericalEMCIC& aCorrFctn) :
137	AliFemtoCorrFctn(),
138	fNumerator(0),
139	fDenominator(0),
140	fEnergyTotalReal(0),
141	fEnergyMultReal(0),
142	fPzMultReal(0),
143	fPtMultReal(0),
144	fEnergyTotalMix (0),
145	fEnergyMultMix (0),
146	fPzMultMix(0),
147	fPtMultMix(0),
148	fPairCut(0x0)
149	{
150	// Copy constructor
151	fNumerator = new TH3D(*aCorrFctn.fNumerator);
152	fDenominator = new TH3D(*aCorrFctn.fDenominator);
153	fEnergyTotalReal = new TH3D(*aCorrFctn.fEnergyTotalReal);
154	fEnergyMultReal = new TH3D(*aCorrFctn.fEnergyMultReal);
155	fPzMultReal = new TH3D(*aCorrFctn.fPzMultReal);
156	fPtMultReal = new TH3D(*aCorrFctn.fPtMultReal);
157	fEnergyTotalMix = new TH3D(*aCorrFctn.fEnergyTotalMix);
158	fEnergyMultMix = new TH3D(*aCorrFctn.fEnergyMultMix);
159	fPzMultMix = new TH3D(*aCorrFctn.fPzMultMix);
160	fPtMultMix = new TH3D(*aCorrFctn.fPtMultMix);
161	fPairCut = aCorrFctn.fPairCut;
162	}
163	//____________________________
164	AliFemtoCorrFctn3DSphericalEMCIC::~AliFemtoCorrFctn3DSphericalEMCIC(){
165	// Destructor
166	delete fNumerator;
167	delete fDenominator;
168	delete fEnergyTotalReal;
169	delete fEnergyMultReal;
170	delete fPzMultReal;
171	delete fPtMultReal;
172	delete fEnergyTotalMix;
173	delete fEnergyMultMix;
174	delete fPzMultMix;
175	delete fPtMultMix;
176	}
177	//_________________________
178	AliFemtoCorrFctn3DSphericalEMCIC& AliFemtoCorrFctn3DSphericalEMCIC::operator=(const AliFemtoCorrFctn3DSphericalEMCIC& aCorrFctn)
179	{
180	// assignment operator
181	if (this == &aCorrFctn)
182	return *this;
183
184	if (fNumerator) delete fNumerator;
185	fNumerator = new TH3D(*aCorrFctn.fNumerator);
186	if (fDenominator) delete fDenominator;
187	fDenominator = new TH3D(*aCorrFctn.fDenominator);
188	if (fEnergyTotalReal) delete fEnergyTotalReal;
189	fEnergyTotalReal = new TH3D(*aCorrFctn.fEnergyTotalReal);
190	if (fEnergyMultReal) delete fEnergyMultReal;
191	fEnergyMultReal = new TH3D(*aCorrFctn.fEnergyMultReal);
192	if (fPzMultReal) delete fPzMultReal;
193	fPzMultReal = new TH3D(*aCorrFctn.fPzMultReal);
194	if (fPtMultReal) delete fPtMultReal;
195	fPtMultReal = new TH3D(*aCorrFctn.fPtMultReal);
196	if (fEnergyTotalMix) delete fEnergyTotalMix;
197	fEnergyTotalMix = new TH3D(*aCorrFctn.fEnergyTotalMix);
198	if (fEnergyMultMix) delete fEnergyMultMix;
199	fEnergyMultMix = new TH3D(*aCorrFctn.fEnergyMultMix);
200	if (fPzMultMix) delete fPzMultMix;
201	fPzMultMix = new TH3D(*aCorrFctn.fPzMultMix);
202	if (fPtMultMix) delete fPtMultMix;
203	fPtMultMix = new TH3D(*aCorrFctn.fPtMultMix);
204	fPairCut = aCorrFctn.fPairCut;
205
206	return *this;
207	}
208
209	//_________________________
210	void AliFemtoCorrFctn3DSphericalEMCIC::WriteOutHistos(){
211	// Write out all histograms to file
212	fNumerator->Write();
213	fDenominator->Write();
214	fEnergyTotalReal->Write();
215	fEnergyMultReal->Write();
216	fPzMultReal->Write();
217	fPtMultReal->Write();
218	fEnergyTotalMix->Write();
219	fEnergyMultMix->Write();
220	fPzMultMix->Write();
221	fPtMultMix->Write();
222	}
223	//______________________________
224	TList* AliFemtoCorrFctn3DSphericalEMCIC::GetOutputList()
225	{
226	// Prepare the list of objects to be written to the output
227	TList *tOutputList = new TList();
228
229	tOutputList->Add(fNumerator);
230	tOutputList->Add(fDenominator);
231	tOutputList->Add(fEnergyTotalReal);
232	tOutputList->Add(fEnergyMultReal);
233	tOutputList->Add(fPzMultReal);
234	tOutputList->Add(fPtMultReal);
235	tOutputList->Add(fEnergyTotalMix );
236	tOutputList->Add(fEnergyMultMix );
237	tOutputList->Add(fPzMultMix);
238	tOutputList->Add(fPtMultMix);
239	return tOutputList;
240	}
241
242	//_________________________
243	void AliFemtoCorrFctn3DSphericalEMCIC::Finish(){
244	// here is where we should normalize, fit, etc...
245	}
246
247	//____________________________
248	AliFemtoString AliFemtoCorrFctn3DSphericalEMCIC::Report(){
249	// Construct the report
250	string stemp = "PRF Frame SphericalEMCIC 3D Correlation Function Report:\n";
251	char ctemp[100];
252	sprintf(ctemp,"Number of entries in numerator:\t%E\n",fNumerator->GetEntries());
253	stemp += ctemp;
254	sprintf(ctemp,"Number of entries in denominator:\t%E\n",fDenominator->GetEntries());
255	stemp += ctemp;
256
257	if (fPairCut){
258	sprintf(ctemp,"Here is the PairCut specific to this CorrFctn\n");
259	stemp += ctemp;
260	stemp += fPairCut->Report();
261	}
262	else{
263	sprintf(ctemp,"No PairCut specific to this CorrFctn\n");
264	stemp += ctemp;
265	}
266
267	//
268	AliFemtoString returnThis = stemp;
269	return returnThis;
270	}
271	//____________________________
272	void AliFemtoCorrFctn3DSphericalEMCIC::AddRealPair( AliFemtoPair* pair){
273	// perform operations on real pairs
274	if (fPairCut){
275	if (!(fPairCut->Pass(pair))) return;
276	}
277
278	// //Changed K to Q to be in LCMS, N. Bock
279	double tKO = pair->QOutCMS(); //KOut();
280	double tKS = pair->QSideCMS(); //KSide();
281	double tKL = pair->QLongCMS(); // KLong();
282
283	double tKR = sqrt(tKOtKO + tKStKS + tKL*tKL);
284	double tKC;
285	if ( fabs(tKR) < 1e-10 ) tKC = 0.0;
286	else tKC=tKL/tKR;
287	double tKP=atan2(tKS,tKO);
288
289	fNumerator->Fill(tKR,tKP,tKC);
290
291	// EMCICs
292	AliFemtoLorentzVector tMom1 = pair->Track1()->FourMomentum();
293	AliFemtoLorentzVector tMom2 = pair->Track2()->FourMomentum();
294	double tE1 = tMom1.e();
295	double tE2 = tMom2.e();
296	double tPz1 = tMom1.pz();
297	double tPz2 = tMom2.pz();
298	TVector2 tPt1;
299	TVector2 tPt2;
300	cout <<"AddRealPair" << "E1,E2, K,Ks,KL" <<
301	tE1 <<" " << tE2<<" " << tKR <<" " << tKS <<" " << tKL<<endl;
302	tPt1.Set(tMom1.px(),tMom1.py());
303	tPt2.Set(tMom2.px(),tMom2.py());
304	double tPt1DotPt2=tPt1*tPt2;
305
306	fEnergyTotalReal->Fill(tKR,tKP,tKC,tE1+tE2);
307	fEnergyMultReal->Fill(tKR,tKP,tKC,tE1*tE2);
308	fPzMultReal->Fill(tKR,tKP,tKC,tPz1*tPz2);
309	fPtMultReal->Fill(tKR,tKP,tKC,tPt1DotPt2);
310
311
312	}
313	//____________________________
314	void AliFemtoCorrFctn3DSphericalEMCIC::AddMixedPair( AliFemtoPair* pair){
315	// perform operations on mixed pairs
316	if (fPairCut){
317	if (!(fPairCut->Pass(pair))) return;
318	}
319	// //Changed K to Q to be in LCMS, N. Bock
320	double tKO = pair->QOutCMS(); //KOut();
321	double tKS = pair->QSideCMS(); //KSide();
322	double tKL = pair->QLongCMS(); // KLong();
323
324	double tKR = sqrt(tKOtKO + tKStKS + tKL*tKL);
325	double tKC;
326	if ( fabs(tKR) < 1e-10 ) tKC = 0.0;
327	else tKC=tKL/tKR;
328	double tKP=atan2(tKS,tKO);
329
330	fDenominator->Fill(tKR,tKP,tKC);
331
332	// EMCICs
333	AliFemtoLorentzVector tMom1 = pair->Track1()->FourMomentum();
334	AliFemtoLorentzVector tMom2 = pair->Track2()->FourMomentum();
335	double tE1 = tMom1.e();
336	double tE2 = tMom2.e();
337	double tPz1 = tMom1.pz();
338	double tPz2 = tMom2.pz();
339	TVector2 tPt1;
340	TVector2 tPt2;
341	tPt1.Set(tMom1.px(),tMom1.py());
342	tPt2.Set(tMom2.px(),tMom2.py());
343	double tPt1DotPt2=tPt1*tPt2;
344
345	fEnergyTotalMix->Fill(tKR,tKP,tKC,tE1+tE2);
346	fEnergyMultMix->Fill(tKR,tKP,tKC,tE1*tE2);
347	fPzMultMix->Fill(tKR,tKP,tKC,tPz1*tPz2);
348	fPtMultMix->Fill(tKR,tKP,tKC,tPt1DotPt2);
349
350	}