////////////////////////////////////////////////////////////////////////////////
// //
// File T49Proj2.C //
// Version 2.0 R.Ganz Jan.20 1999 //
// //
// The class T49Proj2BASE provides functionality to generate two particle //
// projections such as Inv mass spectra, HBT correlation .. //
// What type of projection is used, dependents on the derived class //
// which one uses: //
// (eg T49Proj2BP for Bertsch Pratt type HBT Qside,Qout,Qlong 3-dimensional) //
// (eg T49Proj2INVmass for invariant mass spectra 1-dimensional) //
// The T49Proj2BASE class allows for analyses eg in kt and y bins by //
// booking the arrording and managing the appropriate histogram arrays //
// //
////////////////////////////////////////////////////////////////////////////////
/*
$Log: T49Proj2.C,v $
Revision 1.12 2002/02/04 16:16:05 cblume
Fix bug in SetCMS()
Revision 1.11 2001/12/12 17:43:11 cblume
Caluclate k in the pair cm
Revision 1.10 2001/11/27 16:14:52 cblume
Change dimensions of the character arrays
Revision 1.9 2001/11/21 17:02:24 cblume
Add T49Proj2K
Revision 1.8 2001/11/12 12:12:14 cblume
Update for ROOT v3.01
* Revision 1.7 2000/07/07 14:13:29 cblume
* Made ROOT v2.25 compatible
*
* Revision 1.6 1999/12/01 09:39:27 cblume
* Changed to root 2.23-09. Updated the makefiles. Changes in T49CutTrack, T49Dedx, T49Proj1, T49Proj2, TRootDS
*
* Revision 1.5 1999/11/24 16:03:50 cblume
* Fixes in gteval_nt_root.h and T49Mixer. Addrd T49Dedx again in T49ANA. Some cosmetics.
*
* Revision 1.4 1999/11/23 13:54:13 cblume
* Remove files
*
*/
#include "T49Proj2.h"
ClassImp(T49Proj2BP)
ClassImp(T49Proj2YKP)
ClassImp(T49Proj2QtQl)
ClassImp(T49Proj2QINV2)
ClassImp(T49Proj2QINV)
ClassImp(T49Proj2INVmass)
ClassImp(T49Proj2Armenteros)
ClassImp(T49Proj2K)
ClassImp(T49Proj2BASE)
//______________________________________________________________________________
T49Proj2BASE::T49Proj2BASE()
{
//
// Constructor creating default histogram array
//
fSysID = 0;
fBeta = 0.0;
fPart1Mass = 0.0;
fPart2Mass = 0.0;
fVerbose = kFALSE;
fNPairs = 0;
fQinvFlag = 0;
fMirrorFlag = 0;
fDimArray = 0;
fQx = 0.0;
fQy = 0.0;
fQz = 0.0;
fNXBin = 1;
fNYBin = 1;
fNZBin = 1;
fXMin = 0.0;
fYMin = 0.0;
fZMin = 0.0;
fXMax = 0.0;
fYMax = 0.0;
fZMax = 0.0;
SetCMS(158.0);
fBinStat = NULL;
fKtAxis = NULL;
fAverKtHist = NULL;
fYpAxis = NULL;
fAverYpHist = NULL;
fYAxis = NULL;
fAverYHist = NULL;
fPtAxis = NULL;
fAverPtHist = NULL;
fPhiAxis = NULL;
fAverPhiHist = NULL;
fXFAxis = NULL;
fAverXFHist = NULL;
fHist = NULL;
fXHist = NULL;
fYHist = NULL;
fZHist = NULL;
fAHist = NULL;
}
//______________________________________________________________________________
T49Proj2BASE::~T49Proj2BASE()
{
if (fKtAxis) { delete fKtAxis; };
if (fYpAxis) { delete fYpAxis; };
if (fYAxis) { delete fYAxis; };
if (fPtAxis) { delete fPtAxis; };
if (fPhiAxis) { delete fPhiAxis; };
if (fXFAxis) { delete fXFAxis; };
if (fHist) {
fHist->Delete();
delete fHist;
}
if (fXHist) {
fXHist->Delete();
delete fXHist;
}
if (fYHist) {
fYHist->Delete();
delete fYHist;
}
if (fZHist) {
fZHist->Delete();
delete fZHist;
}
if (fAHist) {
fAHist->Delete();
delete fAHist;
}
}
//______________________________________________________________________________
void T49Proj2BASE::SetCMS(Double_t pLab)
{
//
// Set the center of mass system
//
fBetaCMS = pLab / (kNUCLMASS+sqrt(kNUCLMASS2+pLab*pLab));
Double_t gamma = 1.0 / TMath::Sqrt(1.0 + fBetaCMS*fBetaCMS);
fYcms = TMath::Log(gamma+gamma * fBetaCMS);
fPprojCMS = gamma * (pLab - fBetaCMS*TMath::Sqrt(pLab*pLab +kNUCLMASS2));
}
//______________________________________________________________________________
void T49Proj2BASE::SetParticleMass(Int_t a, Int_t m)
{
//
// Set the particle masses of the lists (a = 1,2)
// to m = 1 proton
// 2 kaon
// 3 pion
// 4 electron
//
Float_t mass = 0.0;
switch (TMath::Abs(m)) {
case 1:
mass = kPMASS;
break;
case 2:
mass = kKMASS;
break;
case 3:
mass = kPIMASS;
break;
case 4:
mass = kEMASS;
break;
default:
printf(" <E> Don't know this particle %dn",m);
break;
};
SetParticleMass(a,mass);
}
//______________________________________________________________________________
void T49Proj2BASE::SetParticleMass(Int_t a, Float_t mass)
{
//
// Set the particle masses of the lists (a = 1,2)
//
switch (a) {
case 1:
fPart1Mass = mass;
break;
case 2:
fPart2Mass = mass;
break;
default:
printf(" <E> Don't know where to put this particle %dn",a);
break;
};
}
//______________________________________________________________________________
void T49Proj2BASE::Book(char* nam, char* tit,
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)
{
//
// Books one set of histogram or an array of sets of histogram
// depending on whether SetKtArray(...) or/and SetYpArray(...)
// was activated BEFORE booking
// It calls BookOne(..) for each set
//
Int_t nKt, nYp, nY, nPt, nPhi, nXF ,ID, OF;
Int_t nKtB, nYpB, nYB, nPtB, nPhiB, nXFB;
char TAGname[120] = "";
char TAGAxisName[120] = "";
char TAGtitl[120] = "";
// Further loop for more array types (but keep the order)
nKtB = (fKtAxis ? fKtAxis->GetNbins() : 1);
nYpB = (fYpAxis ? fYpAxis->GetNbins() : 1);
nYB = (fYAxis ? fYAxis->GetNbins() : 1);
nPtB = (fPtAxis ? fPtAxis->GetNbins() : 1);
nPhiB = (fPhiAxis? fPhiAxis->GetNbins() : 1);
nXFB = (fXFAxis ? fXFAxis->GetNbins() : 1);
nKt = 0;
nYp = 0;
nY = 0;
nPt = 0;
nPhi = 0;
nXF = 0;
OF = 0;
fDimArray = nKtB*nYpB*nYB*nPtB*nPhiB*nXFB;
for (ID = 0; ID < fDimArray; ID++) {
OF = 1;
if (fVerbose) {
printf(" <I> Booking %d %d %d %d %d %dn" ,nKt, nYp, nY, nPt, nPhi, nXF);
}
sprintf(TAGname,"%s%d",nam,ID);
sprintf(TAGtitl," %s ",tit);
if (fKtAxis) {
sprintf(TAGAxisName,"%sKtAxis",nam);
fKtAxis->SetName(TAGAxisName);
char ktChar[40] = "";
sprintf(ktChar," %.2f < k_{t} < %.2f",fKtAxis->GetBinLowEdge(nKt+1)
,fKtAxis->GetBinLowEdge(nKt+2));
strcat(TAGtitl,ktChar);
if (OF == 1) {
OF = 0;
if (nKt < nKtB-1) {
nKt++;
}
else {
nKt = 0;
OF = 1;
}
}
}
if (fYpAxis) {
sprintf(TAGAxisName,"%sYpAxis",nam);
fYpAxis->SetName(TAGAxisName);
char ypChar[40] = "";
sprintf(ypChar," %.1f < y < %.1f",fYpAxis->GetBinLowEdge(nYp+1)
,fYpAxis->GetBinLowEdge(nYp+2));
strcat(TAGtitl,ypChar);
if (OF == 1) {
OF = 0;
if (nYp < nYpB-1) {
nYp++;
}
else {
nYp = 0;
OF = 1;
}
}
}
if (fYAxis) {
sprintf(TAGAxisName,"%sYAxis",nam);
fYAxis->SetName(TAGAxisName);
char yChar[40] = "";
sprintf(yChar," %.1f < y < %.1f",fYAxis->GetBinLowEdge(nY+1)
,fYAxis->GetBinLowEdge(nY+2));
strcat(TAGtitl,yChar);
if (OF == 1) {
OF = 0;
if (nY < nYB-1) {
nY++;
}
else {
nY = 0;
OF = 1;
}
}
}
if (fPtAxis) {
sprintf(TAGAxisName,"%sPtAxis",nam);
fPtAxis->SetName(TAGAxisName);
char PtChar[40] = "";
sprintf(PtChar," %.2f < p_{t} < %.2f",fPtAxis->GetBinLowEdge(nPt+1)
,fPtAxis->GetBinLowEdge(nPt+2));
strcat(TAGtitl,PtChar);
if (OF == 1) {
OF = 0;
if (nPt < nPtB-1) {
nPt++;
}
else {
nPt = 0;
OF = 1;
}
}
}
if (fPhiAxis) {
sprintf(TAGAxisName,"%sPhiAxis",nam);
fPhiAxis->SetName(TAGAxisName);
char PhiChar[40] = "";
sprintf(PhiChar," %.3f < #phi < %.3f",fPhiAxis->GetBinLowEdge(nPhi+1)
,fPhiAxis->GetBinLowEdge(nPhi+2));
strcat(TAGtitl,PhiChar);
if (OF == 1) {
OF = 0;
if (nPhi < nPhiB-1) {
nPhi++;
}
else {
nPhi = 0;
OF = 1;
}
}
}
if (fXFAxis) {
sprintf(TAGAxisName,"%sXFAxis",nam);
fXFAxis->SetName(TAGAxisName);
char XFChar[40] = "";
sprintf(XFChar," %.3f < x_{F} < %.3f",fXFAxis->GetBinLowEdge(nXF+1)
,fXFAxis->GetBinLowEdge(nXF+2));
strcat(TAGtitl,XFChar);
if (OF == 1) {
OF = 0;
if (nXF < nXFB-1) {
nXF++;
}
else {
nXF = 0;
OF = 1;
}
}
}
BookOne(TAGname,TAGtitl,nx,xmin,xmax,ny,ymin,ymax,nz,zmin,zmax);
}
char AverName[40];
sprintf(AverName,"%sBinStat",nam);
fBinStat = new TH1D(AverName,"Number of entries distribution"
,fDimArray,-0.5,1.0*(fDimArray-0.5));
if (fKtAxis) {
sprintf(AverName,"%sAverKt",nam);
fAverKtHist = new TProfile(AverName,"Average K_{t}",fDimArray
,-0.5,1.0*(fDimArray-0.5),0.0,5.0);
}
if (fYpAxis) {
sprintf(AverName,"%sAverYp",nam);
fAverYpHist = new TProfile(AverName,"Average y_{p}",fDimArray
,-0.5,1.0*(fDimArray-0.5),0.0,6.0);
}
if (fYAxis) {
sprintf(AverName,"%sAverY",nam);
fAverYHist = new TProfile(AverName,"Average y",fDimArray
,-0.5,1.0*(fDimArray-0.5),0.0,6.0);
}
if (fPtAxis) {
sprintf(AverName,"%sAverPt",nam);
fAverPtHist = new TProfile(AverName,"Average p_{t}",fDimArray
,-0.5,1.0*(fDimArray-0.5),0.0,5.0);
}
if (fPhiAxis) {
sprintf(AverName,"%sAverPhi",nam);
fAverPhiHist= new TProfile(AverName,"Average #phi",fDimArray
,-0.5,1.0*(fDimArray-0.5),0.0,360.0);
}
if (fXFAxis) {
sprintf(AverName,"%sAverXF",nam);
fAverXFHist = new TProfile(AverName,"Average x_{F}",fDimArray
,-0.5,1.0*(fDimArray-0.5),0.0,1.5);
}
}
//______________________________________________________________________________
void T49Proj2BASE::BookOne(char* nam, char* tit,
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)
{
//
// Books one set of histograms
// Each set contains dependent of the dimension of the projection:
// H Histogram for the projected quantiti
// X Histogram accumulating the X component of the projected quantity (eg Qside for BP HBT)
// Y Histogram accumulating the Y component (only for 2-D and 3-D projections)
// Z Histogram accumulating the Z component (only 3-D projections)
//
// If fQinvFlag was set to 1 it also contains a histogram accumulating Qinv
//
// Dividing eg. X by H results in the average X component in each histgram bin
// of the projection
//
// The input parameter <nam> and <tit> will be part of the title and the ROOT name
//
fNXBin = nx;
fXMin = xmin;
fXMax = xmax;
// Create the histogram arrays
if (!fHist) fHist = new TObjArray();
if (!fXHist) fXHist = new TObjArray();
if (!fYHist) fYHist = new TObjArray();
if (!fZHist) fZHist = new TObjArray();
if (!fAHist) fAHist = new TObjArray();
if (fDim > 1) {
fNYBin = ny;
fYMin = ymin;
fYMax = ymax;
}
else {
fNYBin = 1;
fYMin = 0.0;
fYMax = 0.0;
}
if (fDim > 2) {
fNZBin = nz;
fZMin = zmin;
fZMax = zmax;
}
else {
fNZBin = 1;
fZMin = 0.0;
fZMax = 0.0;
}
if (fDim > 3) {
printf(" <E> Error defining dimension of projection. n");
}
char title[120] = "";
char name[120] = "";
switch (fDim) {
case 1:
if (fVerbose) printf(" <I> 1D object %s %s n",nam,fHistTitle);
sprintf(name,"H%s",nam);
sprintf(title,"%s %s",fHistTitle,tit);
fHist->Add(new TH1D(name,title,fNXBin,fXMin,fXMax));
((TH1D*) (fHist->At(fHist->GetLast())))->SetXTitle(fXLabel);
((TH1D*) (fHist->At(fHist->GetLast())))->Sumw2();
sprintf(name,"X%s", nam);
sprintf(title,"<%s> %s %s",fXLabel,fHistTitle,tit);
fXHist->Add(new TH1D(name,title,fNXBin,fXMin,fXMax));
if (fQinvFlag) {
sprintf(name,"A%s",nam);
sprintf(title,"<Q_{inv}> %s %s",fHistTitle,tit);
fAHist->Add(new TH1D(name,title,fNXBin,fXMin,fXMax));
}
break;
case 2:
if (fVerbose) printf(" <I> 2D object %s %s n",nam,fHistTitle);
sprintf(name,"H%s",nam);
sprintf(title,"%s %s",fHistTitle,tit);
fHist->Add(new TH2D(name,title,fNXBin,fXMin,fXMax
,fNYBin,fYMin,fYMax));
((TH2D*) (fHist->At(fHist->GetLast())))->SetXTitle(fXLabel);
((TH2D*) (fHist->At(fHist->GetLast())))->SetYTitle(fYLabel);
((TH2D*) (fHist->At(fHist->GetLast())))->Sumw2();
sprintf(name,"X%s",nam);
sprintf(title,"<%s> %s %s",fXLabel,fHistTitle,tit);
fXHist->Add(new TH2D(name,title,fNXBin,fXMin,fXMax
,fNYBin,fYMin,fYMax));
sprintf(name,"Y%s",nam);
sprintf(title,"<%s> %s %s",fYLabel,fHistTitle,tit);
fYHist->Add(new TH2D(name,title,fNXBin,fXMin,fXMax
,fNYBin,fYMin,fYMax));
if (fQinvFlag) {
sprintf(name,"A%s",nam);
sprintf(title,"<Q_{inv}> %s %s",fHistTitle,tit);
fAHist->Add(new TH2D(name,title,fNXBin,fXMin,fXMax
,fNYBin,fYMin,fYMax));
}
break;
case 3:
if (fVerbose) printf(" <I> 3D object %s %s n",nam,fHistTitle);
sprintf(name,"H%s",nam);
sprintf(title,"%s %s",fHistTitle,tit);
fHist->Add(new TH3D(name,title,fNXBin,fXMin,fXMax
,fNYBin,fYMin,fYMax
,fNZBin,fZMin,fZMax));
((TH3D*) (fHist->At(fHist->GetLast())))->SetXTitle(fXLabel);
((TH3D*) (fHist->At(fHist->GetLast())))->SetYTitle(fYLabel);
((TH3D*) (fHist->At(fHist->GetLast())))->SetZTitle(fZLabel);
((TH3D*) (fHist->At(fHist->GetLast())))->Sumw2();;
sprintf(name,"X%s",nam);
sprintf(title,"<%s> %s %s",fXLabel,fHistTitle,tit);
fXHist->Add(new TH3D(name,title,fNXBin,fXMin,fXMax
,fNYBin,fYMin,fYMax
,fNZBin,fZMin,fZMax));
sprintf(name,"Y%s",nam);
sprintf(title,"<%s> %s %s",fYLabel,fHistTitle,tit);
fYHist->Add(new TH3D(name,title,fNXBin,fXMin,fXMax
,fNYBin,fYMin,fYMax
,fNZBin,fZMin,fZMax));
sprintf(name,"Z%s",nam);
sprintf(title,"<%s> %s %s",fZLabel,fHistTitle,tit);
fZHist->Add(new TH3D(name,title,fNXBin,fXMin,fXMax
,fNYBin,fYMin,fYMax
,fNZBin,fZMin,fZMax));
if (fQinvFlag) {
sprintf(name,"A%s",nam);
sprintf(title,"<Q_{inv}> %s %s",fHistTitle,tit);
fAHist->Add(new TH3D(name,title,fNXBin,fXMin,fXMax
,fNYBin,fYMin,fYMax
,fNZBin,fZMin,fZMax));
}
break;
};
}
//______________________________________________________________________________
Int_t T49Proj2BASE::CheckMirror(T49ParticleRoot* t1,T49ParticleRoot* t2)
{
//
// Checks whether the pair is below or above cms-rapidity
//
Double_t y = 0.25 * (TMath::Log((t1->GetE(fPart1Mass)+t1->GetPz())
/(t1->GetE(fPart1Mass)-t1->GetPz()))
+ TMath::Log((t2->GetE(fPart2Mass)+t2->GetPz())
/(t2->GetE(fPart2Mass)-t2->GetPz())));
return ((y < fYcms) ? 1 : 0);
}
//______________________________________________________________________________
void T49Proj2BASE::Fill(TObjArray* ListA, TObjArray* ListB, T49Cut2Track* T2Cut)
{
//
// Filling the projection by combining each particle I
// from list A with each particle J of list B
// If the lists are identical only I<J pairs are combined
// Careful ! If the lists are not identical
// but contain the same particle double counting will appear !!
//
T49ParticleRoot *Part1;
T49ParticleRoot *Part2;
Int_t SameListFlag = 0;
Int_t IP1 = 0;
Int_t IP2 = 0;
if (ListA == ListB) {
SameListFlag = 1;
}
IP1 = 0;
TIter NextParticle(ListA);
while ((Part1 = ((T49ParticleRoot *) NextParticle()))) {
IP1++;
IP2 = 0;
TIter NextP2(ListB);
while ((Part2 = ((T49ParticleRoot *) NextP2()))) {
if ((SameListFlag == 0) && (Part1 == Part2)) {
printf(" <E> Identical track in two different lists!n");
printf(" Double counting in projection!n");
}
IP2++;
if ((SameListFlag == 0) || (IP1 > IP2)) {
if ((T2Cut == NULL) || (T2Cut->CheckPair(Part1,Part2))) {
if ((Part1->GetPx() == Part2->GetPx()) &&
(Part1->GetPy() == Part2->GetPy()) &&
(Part1->GetPz() == Part2->GetPz()) &&
(SameListFlag == 0)) {
printf(" <E> Identical track in two different lists!n");
printf(" Double counting in projection!n");
}
Fill(Part1,Part2);
}
}
}
}
}
//______________________________________________________________________________
void T49Proj2BASE::Fill(T49ParticleRoot* t1, T49ParticleRoot* t2)
{
//
// Fill are single pair into projection
// (at the proper index of the projection)
//
Double_t Qinv = 0.0;
Int_t MirrorThis = 0;
Int_t ID;
if (fMirrorFlag == 1) {
MirrorThis = CheckMirror(t1,t2);
}
// Get the Index for Hist arrays
ID = GetArrayIndex(t1,t2,MirrorThis);
if ((ID > -1) && (ID < fDimArray)) {
// fSysID == 2 is LCMS: beta derived from pair
if (fSysID == 2) {
fBeta = (Double_t) ((t1->GetPz() + t2->GetPz())
/(t1->GetE(fPart1Mass) + t2->GetE(fPart2Mass)));
if (fVerbose == 2) printf(" <I> beta LCMS: %f n",fBeta);
}
QCalc(t1,t2,MirrorThis);
// Here for non id part a more general qinv expr. is needed
if (fQinvFlag) {
Qinv = T49Tool::Qinv(t1,t2,fPart1Mass,fPart2Mass);
}
switch (fDim) {
case 1:
if ((fXMin < fQx ) &&
(fQx < fXMax)) {
fNPairs++;
((TH1D*) ( fHist->At(ID)))->Fill(fQx,1.0);
((TH1D*) (fXHist->At(ID)))->Fill(fQx,fQx);
if (fQinvFlag) {
((TH1D*) (fAHist->At(ID)))->Fill(fQx,Qinv);
}
}
break;
case 2:
if ((fXMin < fQx ) &&
(fQx < fXMax) &&
(fYMin < fQy ) &&
(fQy < fYMax)) {
fNPairs++;
((TH2D*) ( fHist->At(ID)))->Fill(fQx,fQy,1.0);
((TH2D*) (fXHist->At(ID)))->Fill(fQx,fQy,fQx);
((TH2D*) (fYHist->At(ID)))->Fill(fQx,fQy,fQy);
if (fQinvFlag) {
((TH2D*) (fAHist->At(ID)))->Fill(fQx,fQy,Qinv);
}
}
break;
case 3:
if ((fXMin < fQx ) &&
(fQx < fXMax) &&
(fYMin < fQy ) &&
(fQy < fYMax) &&
(fZMin < fQz ) &&
(fQz < fZMax)) {
fNPairs++;
((TH3D*) ( fHist->At(ID)))->Fill(fQx,fQy,fQz,1.0);
((TH3D*) (fXHist->At(ID)))->Fill(fQx,fQy,fQz,fQx);
((TH3D*) (fYHist->At(ID)))->Fill(fQx,fQy,fQz,fQy);
((TH3D*) (fZHist->At(ID)))->Fill(fQx,fQy,fQz,fQz);
if (fQinvFlag) {
((TH3D*) (fAHist->At(ID)))->Fill(fQx,fQy,fQz,Qinv);
}
}
break;
default:
printf(" <E> Fill - Some error with dimension occured. n");
break;
};
}
}
//______________________________________________________________________________
Int_t T49Proj2BASE::GetArrayIndex(T49ParticleRoot *t1,T49ParticleRoot *t2
,Int_t MirrorThis)
{
//
// Calculates the histogram index depending on kt or/and y of the pair
// (also y(P), Pt(P), Phi(P) with P=p1+p2 are supported)
// depending on which have been activated (eg. via SetKtArray(...)
// returns -1 if axis was defined but pairs gives value out of range
//
Int_t Bin = 0;
Int_t InLimit = 1;
Float_t Kt = 0;
Float_t Yp = 0;
Float_t Y = 0;
Float_t Pt = 0;
Float_t Phi = 0;
Float_t XF = 0;
// Inverse order compared to booking -> slowest changing of index
if (fXFAxis && InLimit) {
XF = (T49Tool::zLorentz(t1,fPart1Mass,fBetaCMS).Z()
+ T49Tool::zLorentz(t2,fPart2Mass,fBetaCMS).Z()) / fPprojCMS;
if ((fXFAxis->GetXmin() < XF ) &&
(XF < fXFAxis->GetXmax())) {
Bin = Bin * fXFAxis->GetNbins() + (fXFAxis->FindBin(XF)-1);
if (fVerbose == 2) {
printf(" <i> GetArrayIndex : %d XF=%f n",Bin,XF);
}
}
else {
InLimit = 0;
return -1;
}
}
if (fPhiAxis && InLimit) {
if ((t1->GetPx() + t2->GetPx()) == 0.0) {
Phi = ((t1->GetPy() + t2->GetPy()) > 0.0 ? 90.0 :270.0) ;
}
else {
Phi = TMath::ATan((t1->GetPy()+t2->GetPy())
/ (t1->GetPx()+t2->GetPx())) * 57.29578;
if ((t1->GetPx()+t2->GetPx()) < 0.0) {
Phi = Phi + 180.0;
}
if (Phi < 0.0) {
Phi = Phi + 360.0;
}
}
if ((fPhiAxis->GetXmin() < Phi ) &&
(Phi < fPhiAxis->GetXmax())) {
Bin = Bin * fPhiAxis->GetNbins() + (fPhiAxis->FindBin(Phi)-1);
if (fVerbose == 2) {
printf(" <i> GetArrayIndex : %d Phi=%f n",Bin,Phi);
}
}
else {
InLimit = 0;
return -1;
}
}
if (fPtAxis && InLimit) {
Pt = TMath::Sqrt(T49Tool::Px(t1,t2)*T49Tool::Px(t1,t2)
+ T49Tool::Py(t1,t2)*T49Tool::Py(t1,t2));
if ((fPtAxis->GetXmin() < Pt ) &&
(Pt < fPtAxis->GetXmax())) {
Bin = Bin * fPtAxis->GetNbins() + (fPtAxis->FindBin(Pt)-1);
if (fVerbose == 2) {
printf(" <i> GetArrayIndex : %d Pt=%f n",Bin,Pt);
}
}
else {
InLimit = 0;
}
}
if (fYAxis && InLimit) {
Y = 0.5 * TMath::Log((t1->GetE(fPart1Mass)+t2->GetE(fPart2Mass)
+t1->GetPz()+t2->GetPz())
/(t1->GetE(fPart1Mass)+t2->GetE(fPart2Mass)
-t1->GetPz()-t2->GetPz()));
if (MirrorThis) {
Y = 2.0*fYcms - Y;
}
if ((fYAxis->GetXmin() < Y ) &&
(Y < fYAxis->GetXmax())) {
Bin = Bin * fYAxis->GetNbins() + (fYAxis->FindBin(Y)-1);
if (fVerbose == 2) {
printf(" <i> GetArrayIndex : %d Y=%f n",Bin,Y);
}
}
else {
InLimit = 0;
return -1;
}
}
if (fYpAxis && InLimit) {
Yp = 0.25 * (TMath::Log((t1->GetE(fPart1Mass)+t1->GetPz())
/(t1->GetE(fPart1Mass)-t1->GetPz()))
+TMath::Log((t2->GetE(fPart2Mass)+t2->GetPz())
/(t2->GetE(fPart2Mass)-t2->GetPz())));
if (MirrorThis) {
Yp = 2.0*fYcms - Yp;
}
if ((fYpAxis->GetXmin() < Yp) &&
(Yp < fYpAxis->GetXmax())) {
Bin = Bin * fYpAxis->GetNbins() + (fYpAxis->FindBin(Yp)-1);
if (fVerbose == 2) {
printf(" <i> GetArrayIndex : %d yp=%f n",Bin,Yp);
}
}
else {
InLimit = 0;
return -1;
}
}
if (fKtAxis && InLimit) {
Kt = 0.5*TMath::Sqrt(T49Tool::Px(t1,t2)*T49Tool::Px(t1,t2)
+T49Tool::Py(t1,t2)*T49Tool::Py(t1,t2));
if ((fKtAxis->GetXmin() < Kt) &&
(Kt < fKtAxis->GetXmax())) {
Bin = Bin * fKtAxis->GetNbins() + (fKtAxis->FindBin(Kt)-1);
if (fVerbose == 2) {
printf(" <i> GetArrayIndex : %d kt=%f n",Bin,Kt);
}
}
else {
InLimit = 0;
return -1;
}
}
if (Bin > (fDimArray - 1)) {
printf(" <E> Found Array Index OUT OF LIMIT %d (%d)n",Bin,fDimArray);
Bin = -1;
}
// Here add further axis but keep reverse order as in Book()
if (InLimit == 0) {
Bin = -1;
return -1;
}
else {
fBinStat->Fill((Double_t) Bin,1.0);
if (fKtAxis) {
fAverKtHist->Fill((Double_t) Bin,(Double_t) Kt);
}
if (fYpAxis) {
fAverYpHist->Fill((Double_t) Bin,(Double_t) Yp);
}
if (fYAxis) {
fAverYHist->Fill((Double_t) Bin,(Double_t) Y);
}
if (fPtAxis) {
fAverPtHist->Fill((Double_t) Bin,(Double_t) Pt);
}
if (fPhiAxis) {
fAverPhiHist->Fill((Double_t) Bin,(Double_t) Phi);
}
if (fXFAxis) {
fAverXFHist->Fill((Double_t) Bin,(Double_t) XF);
}
}
return Bin;
}
//______________________________________________________________________________
void T49Proj2BASE::WriteHist()
{
//
// Write histgrams to file
// (which has to be open before using this function)
//
fHist->Write();
fXHist->Write();
if (fYHist) {
fYHist->Write();
}
if (fZHist) {
fZHist->Write();
}
if (fAHist) {
fAHist->Write();
}
fBinStat->Write();
if (fKtAxis) {
fAverKtHist->Write();
fKtAxis->Write();
}
if (fYpAxis) {
fAverYpHist->Write();
fYpAxis->Write();
}
if (fYAxis) {
fAverYHist->Write();
fYAxis->Write();
}
if (fPtAxis) {
fAverPtHist->Write();
fPtAxis->Write();
}
if (fPhiAxis) {
fAverPhiHist->Write();
fPhiAxis->Write();
}
if (fXFAxis) {
fAverXFHist->Write();
fXFAxis->Write();
}
}
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.