From dd82cadcabd7d7498b1652c0136ecaf406c3c50c Mon Sep 17 00:00:00 2001 From: skowron Date: Tue, 28 Oct 2003 17:59:19 +0000 Subject: [PATCH] CRAB added --- HBTAN/AliHBTCorrFitFctn.cxx | 2 +- HBTAN/AliHBTCorrelFctn.h | 6 +- HBTAN/AliHBTCrab.cxx | 347 ++++++++ HBTAN/AliHBTCrab.h | 66 ++ HBTAN/AliHBTFunction.cxx | 57 +- HBTAN/AliHBTFunction.h | 4 +- HBTAN/AliHBTLLWeightTheorFctn.h | 82 -- HBTAN/AliHBTLLWeights.cxx | 18 +- HBTAN/AliHBTLLWeights.h | 5 +- HBTAN/AliHBTPair.cxx | 20 +- HBTAN/AliHBTPair.h | 8 +- HBTAN/AliHBTParticle.h | 3 + HBTAN/AliHBTReaderESD.cxx | 5 +- ...TLLWeightFctn.cxx => AliHBTWeightFctn.cxx} | 193 +++-- ...liHBTLLWeightFctn.h => AliHBTWeightFctn.h} | 70 +- ...heorFctn.cxx => AliHBTWeightTheorFctn.cxx} | 85 +- HBTAN/AliHBTWeightTheorFctn.h | 100 +++ HBTAN/AliHBTWeights.cxx | 11 + HBTAN/AliHBTWeights.h | 23 + ...TLLWeightsPID.cxx => AliHBTWeightsPID.cxx} | 14 +- ...liHBTLLWeightsPID.h => AliHBTWeightsPID.h} | 24 +- HBTAN/HBTAnalysisLinkDef.h | 35 +- HBTAN/libHBTAN.pkg | 6 +- HBTAN/volya_complex.h | 799 ++++++++++++++++++ 24 files changed, 1668 insertions(+), 315 deletions(-) create mode 100644 HBTAN/AliHBTCrab.cxx create mode 100644 HBTAN/AliHBTCrab.h delete mode 100644 HBTAN/AliHBTLLWeightTheorFctn.h rename HBTAN/{AliHBTLLWeightFctn.cxx => AliHBTWeightFctn.cxx} (64%) rename HBTAN/{AliHBTLLWeightFctn.h => AliHBTWeightFctn.h} (67%) rename HBTAN/{AliHBTLLWeightTheorFctn.cxx => AliHBTWeightTheorFctn.cxx} (57%) create mode 100644 HBTAN/AliHBTWeightTheorFctn.h create mode 100644 HBTAN/AliHBTWeights.cxx create mode 100644 HBTAN/AliHBTWeights.h rename HBTAN/{AliHBTLLWeightsPID.cxx => AliHBTWeightsPID.cxx} (95%) rename HBTAN/{AliHBTLLWeightsPID.h => AliHBTWeightsPID.h} (66%) create mode 100644 HBTAN/volya_complex.h diff --git a/HBTAN/AliHBTCorrFitFctn.cxx b/HBTAN/AliHBTCorrFitFctn.cxx index 8c44d77eb52..36614021c0d 100644 --- a/HBTAN/AliHBTCorrFitFctn.cxx +++ b/HBTAN/AliHBTCorrFitFctn.cxx @@ -34,7 +34,7 @@ void AliHBTCorrFitFctn::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTP Double_t q = trackpair->GetQInv(); Bool_t fill = kFALSE; - Double_t weight = partpair->GetLLWeight(); + Double_t weight = partpair->GetWeight(); fNumerator->Fill(q,weight); if ( (q < 0.15) && (fNPairsFitArea < 2.e+5)) diff --git a/HBTAN/AliHBTCorrelFctn.h b/HBTAN/AliHBTCorrelFctn.h index bf9ea38f884..88dce65b475 100644 --- a/HBTAN/AliHBTCorrelFctn.h +++ b/HBTAN/AliHBTCorrelFctn.h @@ -37,7 +37,7 @@ class AliHBTQOutCMSLCCorrelFctn: public AliHBTOnePairFctn1D virtual ~AliHBTQOutCMSLCCorrelFctn(){}; TH1* GetResult(); protected: - Double_t GetValue(AliHBTPair * pair){return TMath::Abs(pair->GetQOutCMSLC());} + Double_t GetValue(AliHBTPair * pair){return pair->GetQOutCMSLC();} public: ClassDef(AliHBTQOutCMSLCCorrelFctn,1) @@ -53,7 +53,7 @@ class AliHBTQLongCMSLCCorrelFctn: public AliHBTOnePairFctn1D virtual ~AliHBTQLongCMSLCCorrelFctn(){}; TH1* GetResult(); protected: - Double_t GetValue(AliHBTPair * pair){return TMath::Abs(pair->GetQLongCMSLC());} + Double_t GetValue(AliHBTPair * pair){return pair->GetQLongCMSLC();} public: ClassDef(AliHBTQLongCMSLCCorrelFctn,1) @@ -69,7 +69,7 @@ class AliHBTQSideCMSLCCorrelFctn: public AliHBTOnePairFctn1D virtual ~AliHBTQSideCMSLCCorrelFctn(){} TH1* GetResult(); protected: - Double_t GetValue(AliHBTPair * pair){return TMath::Abs(pair->GetQSideCMSLC());} + Double_t GetValue(AliHBTPair * pair){return pair->GetQSideCMSLC();} public: ClassDef(AliHBTQSideCMSLCCorrelFctn,1) diff --git a/HBTAN/AliHBTCrab.cxx b/HBTAN/AliHBTCrab.cxx new file mode 100644 index 00000000000..db1070ca60f --- /dev/null +++ b/HBTAN/AliHBTCrab.cxx @@ -0,0 +1,347 @@ +#include "AliHBTCrab.h" +//______________________________________________________________________ +///////////////////////////////////////////////////////////////////////// +// // +// AliRoot wrapper to CRAB // +// taken from http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html // +// written by Scott Pratt // +// // +///////////////////////////////////////////////////////////////////////// + +#include "AliHBTPair.h" + +#include +#include + +#include "volya_complex.h" + +//#include +//using namespace std; +AliHBTCrab* AliHBTCrab::fgCrab = 0x0; + +const Double_t AliHBTCrab::fgkWcons = 1./0.1973; +const Double_t AliHBTCrab::fgkROOT2=1.41421356237309504880; +const double_complex AliHBTCrab::ci(0.0,1.0); + +/************************************************************/ + +AliHBTCrab* AliHBTCrab::Instance() +{ +// returns instance of class + if (fgCrab == 0x0) + { + fgCrab = new AliHBTCrab(); + } + return fgCrab; +} +//=================================================================== +void AliHBTCrab::Set() +{ + //sets this as weighitng class + Info("Set","Setting CRAB as Weighing Class"); + + if ( fgWeights == 0x0 ) + { + fgWeights = AliHBTCrab::Instance(); + return; + } + if ( fgWeights == AliHBTCrab::Instance() ) return; + delete fgWeights; + fgWeights = AliHBTCrab::Instance(); +} +//=================================================================== + +AliHBTCrab::AliHBTCrab(): +fBreitWigner(kFALSE), +fReducedMom(kTRUE), +fMaxMomentum(100.0) +{ + //ctor +} +//=================================================================== +void AliHBTCrab::Init(Int_t pid1,Int_t pid2) +{ + MASS1 = TDatabasePDG::Instance()->GetParticle(pid1)->Mass(); + MASS2 = TDatabasePDG::Instance()->GetParticle(pid2)->Mass(); + INTERACTION_WSYM = 1.0; + INTERACTION_WANTI = 0.0; + INTERACTION_WNOSYM = 0.0; + INTERACTION_DELK = 1.0; + INTERACTION_NKMAX = 100; + + fPid1 = pid1; + fPid2 = pid2; + +} + +Bool_t AliHBTCrab::SetConfig(const AliHBTPair* pair) +{ + Int_t pdg1 = pair->Particle1()->GetPdgCode(); + Int_t pdg2 = pair->Particle2()->GetPdgCode(); + + if ( ( pdg1 == fPid1) && ( pdg2 == fPid2) ) return kFALSE; + else Init (pdg1,pdg2); + + return kTRUE; +} +//=================================================================== + +Double_t AliHBTCrab::GetWeight(const AliHBTPair* partpair) +{ + Double_t qred, r, qdotr, mom; + Int_t test; + + SetConfig(partpair); + + get_com_quantities(partpair, &qred, &r, &qdotr, &mom, &test); + if(test==0) + { + Info("GetWeight","Test is 0"); + } + Double_t corr = corrcalc(qred,qdotr,r); + + return corr; +} +//=================================================================== + +void AliHBTCrab::get_com_quantities(const AliHBTPair* pair, + double *qred,double *r,double *qdotr,double *mom, int *test) + { +//************************************ +// ALICE // + + double p1[4]; + double p2[4]; + double r1[4]; + double r2[4]; + + static const Double_t cmtofm = 1.e13; + static const Double_t cmtoOneOverGeV = cmtofm*fgkWcons; + + AliHBTParticle *part1 = pair->Particle1(); + AliHBTParticle *part2 = pair->Particle2(); + + p1[0] = part1->Energy()*1000.0; + p1[1] = part1->Px()*1000.0; + p1[2] = part1->Py()*1000.0; + p1[3] = part1->Pz()*1000.0; + + p2[0] = part2->Energy()*1000.0; + p2[1] = part2->Px()*1000.0; + p2[2] = part2->Py()*1000.0; + p2[3] = part2->Pz()*1000.0; + + r1[0] = part1->T(); + r1[1] = part1->Vx()*cmtofm; + r1[2] = part1->Vy()*cmtofm; + r1[3] = part1->Vz()*cmtofm; + + r2[0] = part2->T(); + r2[1] = part2->Vx()*cmtofm; + r2[2] = part2->Vy()*cmtofm; + r2[3] = part2->Vz()*cmtofm; + +// END OF ALICE STUFF + +// This code is written by Scott Pratt +// taken from http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html + int alpha; + double kdotr; + if (fReducedMom) + { + const double momtest=4.0*fMaxMomentum*fMaxMomentum; + } + else + { + const double momtest=fMaxMomentum*fMaxMomentum; + } + + double ptot2,pdotr,pp,rr; + + if ( part1->GetPdgCode() == part2->GetPdgCode() ) + { + *test=1; + *mom=-(p2[0]-p1[0])*(p2[0]-p1[0]); + for(alpha=1;alpha<4;alpha++){ + *mom=*mom+(p2[alpha]-p1[alpha])*(p2[alpha]-p1[alpha]); + } + //#if ! defined MIXED_PAIRS_FOR_DENOM + // if(*mom>momtest){ + // *test=0; + // return; + // } + //#endif + pp=(p1[0]+p2[0]); + rr=(r2[0]-r1[0]); + pdotr=pp*rr; + kdotr=(p2[0]-p1[0])*rr; + ptot2=pp*pp; + *r=-rr*rr; + for(alpha=1;alpha<4;alpha++){ + pp=(p1[alpha]+p2[alpha]); + rr=(r2[alpha]-r1[alpha]); + pdotr=pdotr-pp*rr; + kdotr=kdotr-(p2[alpha]-p1[alpha])*rr; + ptot2=ptot2-pp*pp; + *r=*r+rr*rr; + } + *mom=sqrt(*mom); + *qred=0.5**mom; + + if (fReducedMom) + { + *mom=*qred; + } + + *qdotr=0.5*kdotr; + *r=sqrt(*r+pdotr*pdotr/ptot2); + } + else //identical + { + // const double kdotp=MASS2*MASS2-MASS1*MASS1; + const double kdotp = part2->GetMass()*part2->GetMass()- part1->GetMass()*part1->GetMass(); + *test=1; + *mom=-(p2[0]-p1[0])*(p2[0]-p1[0]); + ptot2=(p1[0]+p2[0])*(p1[0]+p2[0]); + for(alpha=1;alpha<4;alpha++){ + *mom=*mom+(p2[alpha]-p1[alpha])*(p2[alpha]-p1[alpha]); + ptot2=ptot2-(p1[alpha]+p2[alpha])*(p1[alpha]+p2[alpha]); + } + *mom=*mom+kdotp*kdotp/ptot2; + //#if ! defined MIXED_PAIRS_FOR_DENOM + // if(*mom>momtest){ + // *test=0; + // return; + // } + //#endif + pp=(p1[0]+p2[0]); + rr=(r2[0]-r1[0]); + pdotr=pp*rr; + kdotr=(p2[0]-p1[0])*rr; + *r=-rr*rr; + for(alpha=1;alpha<4;alpha++){ + pp=(p1[alpha]+p2[alpha]); + rr=(r2[alpha]-r1[alpha]); + pdotr=pdotr-pp*rr; + kdotr=kdotr-(p2[alpha]-p1[alpha])*rr; + *r=*r+rr*rr; + } + kdotr=(-kdotr+kdotp*pdotr/ptot2); + *mom=sqrt(*mom); + *qred=0.5**mom; + + if (fReducedMom) + { + *mom=*qred; + } + + *qdotr=0.5*kdotr; + *r=sqrt(*r+pdotr*pdotr/ptot2); + }//not identical + + return; +} + + +//=================================================================== + +double AliHBTCrab::corrcalc(double trueqred,double trueqdotr,double truer) +{ +//#define REDUCED_MOM +// This code is written by Scott Pratt +// taken from http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html + double eta,arg,krmax,corr0; + double xx,xxprime,xxjj,p1,zk; + int jj,kk,ipart,ipartcount,ispin; + double pfactor,wsym_leftover,wanti_leftover,wnosym_leftover; + double qred,qdotr,r; + const double rmass=MASS1*MASS2/(MASS1+MASS2); + double_complex cphi1,cphi2,cphis,cphia; + + arg=trueqdotr/197.323-2.0*TMath::Pi()*floor(trueqdotr/(197.323*2.0*TMath::Pi())); + cphi1=exp(ci*arg); + cphis=fgkROOT2*real(cphi1); + cphia=ci*fgkROOT2*imag(cphi1); + corr0=real(INTERACTION_WSYM*cphis*conj(cphis) + +INTERACTION_WANTI*cphia*conj(cphia) + +INTERACTION_WNOSYM*cphi1*conj(cphi1)); + goto OUTSIDE_INTERACTION_RANGE; + +#ifdef REDUCED_MOM + kk=(int)floor(trueqred/INTERACTION_DELK); + qred=(0.5+kk)*INTERACTION_DELK; +#else + kk=(int)floor(2.0*trueqred/INTERACTION_DELK); + qred=(0.5+kk)*INTERACTION_DELK/2.0; +#endif + qdotr=trueqdotr*qred/trueqred; + if(kk>=INTERACTION_NKMAX){ + corr0=1.0; + goto OUTSIDE_INTERACTION_RANGE; + } + r=truer; + + eta=0.0; + arg=qdotr/197.323-2.0*TMath::Pi()*floor(qdotr/(197.323*2.0*TMath::Pi())); + cphi1=exp(ci*arg); + cphi2=conj(cphi1); + + cphis=(cphi1+cphi2)/fgkROOT2; + cphia=(cphi1-cphi2)/fgkROOT2; + corr0=0.0; + /* If there are corrections for strong interactions, add the + change for each partial wave. If npartial = 0 then there + are no strong int. corrections. */ + wsym_leftover=INTERACTION_WSYM; + wanti_leftover=INTERACTION_WANTI; + wnosym_leftover=INTERACTION_WNOSYM; + + corr0=corr0+real(wsym_leftover*cphis*conj(cphis) + +wanti_leftover*cphia*conj(cphia) + +wnosym_leftover*cphi1*conj(cphi1)); +OUTSIDE_INTERACTION_RANGE: + +#ifdef BREIT_WIGNER + corr0=corr0+bwcalc(trueqred,truer); +#endif + + return corr0; +} + +double_complex AliHBTCrab::cgamma(double_complex c){ + /* This calc.s gamma functions which are in the form gamma(n+i*y) + where n is an int and y is real. */ +// This code is written by Scott Pratt +// taken from http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html + double_complex cg,cphase; + int mm,j; + double x,y,phase,delp,cgmag; + x=real(c); + y=imag(c); + phase=-TMath::E()*y; + for(j=1;j<=100000;j++){ + delp=(y/(double)j)-atan(y/(double)j); + phase=phase+delp; + if(fabs(delp)<1E-10) goto CGAMMA_ESCAPE; + } + printf("oops not accurate enough, increase jmax\n"); +CGAMMA_ESCAPE: + phase=phase-2.0*TMath::Pi()*floor(phase/(2.0*TMath::Pi())); + cphase=exp(ci*phase); + cgmag=sqrt(TMath::Pi()*y/sinh(TMath::Pi()*y)); + mm=(int)floor(x+0.5); + cg=cgmag*cphase; + if(mm<1){ + for(j=1;j<=-mm+1;j++){ + cg=cg/(1.0+(double)(-j)+ci*y); + } + } + if(mm>1) { + for(j=1;j<=mm-1;j++){ + cg=cg*((double)(j)+ci*y); + } + } + return cg; +} +//=================================================================== + diff --git a/HBTAN/AliHBTCrab.h b/HBTAN/AliHBTCrab.h new file mode 100644 index 00000000000..35ccbca7be7 --- /dev/null +++ b/HBTAN/AliHBTCrab.h @@ -0,0 +1,66 @@ +/* $Id$ */ + +// This class introduces the weight's calculation +// according to the Lednicky's algorithm. +// The detailed description of the algorithm can be found +// in comments to fortran code: +// fsiw.f, fsiini.f + +#ifndef ALIHBTCrab_H +#define ALIHBTCrab_H + +#include "AliHBTWeights.h" + +class Complex; +typedef Complex double_complex; + +class AliHBTPair; + +class AliHBTCrab: public AliHBTWeights + { + public: + + virtual ~AliHBTCrab(){fgCrab =0x0;} + static AliHBTCrab* Instance(); + void Set(); + + Double_t GetWeight(const AliHBTPair* partpair); + void Init(Int_t pid1,Int_t pid2); //put the initial values in fortran commons fsiini, led_bldata + private: + AliHBTCrab(); + AliHBTCrab(const AliHBTCrab &/*source*/); + AliHBTCrab & operator=(const AliHBTCrab& /*source*/); + + void get_com_quantities(const AliHBTPair* pair, double *qred,double *r,double *qdotr,double *mom, int *test); + double corrcalc(double trueqred,double trueqdotr,double truer); + + Bool_t fBreitWigner; + Bool_t fReducedMom; + Float_t fMaxMomentum; + + Bool_t SetConfig(const AliHBTPair* pair); + + Int_t fPid1; + Int_t fPid2; + + Double_t MASS1; + Double_t MASS2; + + Float_t INTERACTION_WSYM;/* fractions of symmetric and antisym weights of the various spin channels */ + Float_t INTERACTION_WANTI; + Float_t INTERACTION_WNOSYM; + + Float_t INTERACTION_DELK; + Int_t INTERACTION_NKMAX;/* number of momentum points in mesh for strong/coul. interaction */ + + static const double_complex ci; + static const Double_t fgkROOT2;//! some const + static const Double_t fgkWcons; //constant for fm->GeV conversion 1/0.1973 + + double_complex cgamma(double_complex c); + + static AliHBTCrab* fgCrab; + ClassDef(AliHBTCrab,1) + }; + +#endif diff --git a/HBTAN/AliHBTFunction.cxx b/HBTAN/AliHBTFunction.cxx index 767834a3fdc..aceda9fbcc6 100644 --- a/HBTAN/AliHBTFunction.cxx +++ b/HBTAN/AliHBTFunction.cxx @@ -342,9 +342,8 @@ Double_t AliHBTFunction1D::Scale(TH1D* num,TH1D* den) } if (AliHBTParticle::GetDebug()>0) Info("Scale","No errors detected"); - Double_t ratio; - Double_t sum = 0; - Int_t n = 0; + Double_t densum = 0.0; + Double_t numsum = 0.0; Int_t offset = nbins - fNBinsToScale - 1; @@ -352,16 +351,16 @@ Double_t AliHBTFunction1D::Scale(TH1D* num,TH1D* den) { if ( num->GetBinContent(i) > 0.0 ) { - ratio = den->GetBinContent(i)/num->GetBinContent(i); - sum += ratio; - n++; + densum = fDenominator->GetBinContent(i); + numsum = fNumerator->GetBinContent(i); } } - if(AliHBTParticle::GetDebug() > 0) Info("Scale","sum=%f fNBinsToScale=%d n=%d",sum,fNBinsToScale,n); + if(AliHBTParticle::GetDebug() > 0) + Info("Scale","numsum=%f densum=%f fNBinsToScaleX=%d",numsum,densum,fNBinsToScale); - if (n == 0) return 0.0; - Double_t ret = sum/((Double_t)n); + if (numsum == 0) return 0.0; + Double_t ret = densum/numsum; if(AliHBTParticle::GetDebug() > 0) Info("Scale","returning %f",ret); return ret; @@ -538,25 +537,24 @@ Double_t AliHBTFunction2D::Scale() Int_t offsetX = nbinsX - fNBinsToScaleX - 1; //bin that we start loop over bins in axis X Int_t offsetY = nbinsY - fNBinsToScaleY - 1; //bin that we start loop over bins in axis X - Double_t ratio; - Double_t sum = 0; - Int_t n = 0; + Double_t densum = 0.0; + Double_t numsum = 0.0; for (UInt_t j = offsetY; j< nbinsY; j++) for (UInt_t i = offsetX; i< nbinsX; i++) { if ( fNumerator->GetBinContent(i,j) > 0.0 ) { - ratio = fDenominator->GetBinContent(i,j)/fNumerator->GetBinContent(i,j); - sum += ratio; - n++; + densum = fDenominator->GetBinContent(i,j); + numsum = fNumerator->GetBinContent(i,j); } } - if(AliHBTParticle::GetDebug() > 0) Info("Scale","sum=%f fNBinsToScaleX=%d fNBinsToScaleY=%d n=%d",sum,fNBinsToScaleX,fNBinsToScaleY,n); + if(AliHBTParticle::GetDebug() > 0) + Info("Scale","numsum=%f densum=%f fNBinsToScaleX=%d fNBinsToScaleY=%d",numsum,densum,fNBinsToScaleX,fNBinsToScaleY); - if (n == 0) return 0.0; - Double_t ret = sum/((Double_t)n); + if (numsum == 0) return 0.0; + Double_t ret = densum/numsum; if(AliHBTParticle::GetDebug() > 0) Info("Scale","returning %f",ret); return ret; @@ -692,10 +690,10 @@ void AliHBTFunction3D::BuildHistos(Int_t nxbins, Float_t xmax, Float_t xmin, TString denstr = fName + " Denominator";//title and name of the //denominator histogram - fNumerator = new TH3D(numstr.Data(),numstr.Data(), + fNumerator = new TH3F(numstr.Data(),numstr.Data(), nxbins,xmin,xmax,nybins,ymin,ymax,nzbins,zmin,zmax); - fDenominator = new TH3D(denstr.Data(),denstr.Data(), + fDenominator = new TH3F(denstr.Data(),denstr.Data(), nxbins,xmin,xmax,nybins,ymin,ymax,nzbins,zmin,zmax); fNumerator->Sumw2(); @@ -753,9 +751,8 @@ Double_t AliHBTFunction3D::Scale() Int_t offsetY = nbinsY - fNBinsToScaleY - 1; //bin that we start loop over bins in axis Y Int_t offsetZ = nbinsZ - fNBinsToScaleZ - 1; //bin that we start loop over bins in axis Z - Double_t ratio; - Double_t sum = 0; - Int_t n = 0; + Double_t densum = 0.0; + Double_t numsum = 0.0; for (UInt_t k = offsetZ; kGetBinContent(i,j,k) > 0.0 ) { - ratio = fDenominator->GetBinContent(i,j,k)/fNumerator->GetBinContent(i,j,k); - sum += ratio; - n++; + + densum = fDenominator->GetBinContent(i,j,k); + numsum = fNumerator->GetBinContent(i,j,k); } } if(AliHBTParticle::GetDebug() > 0) - Info("Scale","sum=%f fNBinsToScaleX=%d fNBinsToScaleY=%d fNBinsToScaleZ=%d n=%d", - sum,fNBinsToScaleX,fNBinsToScaleY,fNBinsToScaleZ,n); + Info("Scale","numsum=%f densum=%f fNBinsToScaleX=%d fNBinsToScaleY=%d fNBinsToScaleZ=%d", + numsum,densum,fNBinsToScaleX,fNBinsToScaleY,fNBinsToScaleZ); - if (n == 0) return 0.0; - Double_t ret = sum/((Double_t)n); + if (numsum == 0) return 0.0; + Double_t ret = densum/numsum; if(AliHBTParticle::GetDebug() > 0) Info("Scale","returning %f",ret); return ret; diff --git a/HBTAN/AliHBTFunction.h b/HBTAN/AliHBTFunction.h index 1132c405ead..c5f91b3e89d 100644 --- a/HBTAN/AliHBTFunction.h +++ b/HBTAN/AliHBTFunction.h @@ -323,8 +323,8 @@ class AliHBTFunction3D: public AliHBTFunction Int_t nzbins, Float_t zmax, Float_t zmin); virtual void BuildHistos(); - TH3D* fNumerator; // Numerator histogram - TH3D* fDenominator; // Denominator histogram + TH3F* fNumerator; // Numerator histogram + TH3F* fDenominator; // Denominator histogram //definition of area used for scaling -> Scale is calculated this //way that after division tale is on 1 diff --git a/HBTAN/AliHBTLLWeightTheorFctn.h b/HBTAN/AliHBTLLWeightTheorFctn.h deleted file mode 100644 index 7acd78ecd87..00000000000 --- a/HBTAN/AliHBTLLWeightTheorFctn.h +++ /dev/null @@ -1,82 +0,0 @@ -#ifndef ALIHBTLLWEIGHTTHEORFCTN_H -#define ALIHBTLLWEIGHTTHEORFCTN_H -/* $Id$ */ - -//This function allows to obtain Q_inv correlation function with weights -//calculated by Lednicky's alghorithm. -//Numerator is filled with weighted events. Weights are attributed to simulated particles. -//Weights are calculated with corresponding simulated particles momenta. -//Denominator is filled with mixing unweighted simulated particles. -//One needs only simulated pairs, so -//this function is of class AliHBTOnePairFctn1D. -//Author Ludmila Malinina JINR (malinina@sunhe.jinr.ru) - -#include "AliHBTFunction.h" - -class AliHBTLLWeights; - -class AliHBTLLWeightTheorQInvFctn: public AliHBTOnePairFctn1D -{ - - public: - AliHBTLLWeightTheorQInvFctn(Int_t nbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0); - virtual ~AliHBTLLWeightTheorQInvFctn(){} - - TH1* GetResult(); - void ProcessSameEventParticles(AliHBTPair* partpair); - - Double_t GetValue(AliHBTPair* partpair) - { return partpair->GetQInv();} //isn't used - - ClassDef(AliHBTLLWeightTheorQInvFctn,1) -}; - -class AliHBTLLWeightTheorQOutFctn: public AliHBTOnePairFctn1D -{ - - public: - AliHBTLLWeightTheorQOutFctn(Int_t nbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0); - virtual ~AliHBTLLWeightTheorQOutFctn(){} - - TH1* GetResult(); - void ProcessSameEventParticles(AliHBTPair* partpair); - - Double_t GetValue(AliHBTPair* partpair) - { return partpair->GetQOutCMSLC();} //isn't used - - ClassDef(AliHBTLLWeightTheorQOutFctn,1) -}; - -class AliHBTLLWeightTheorQSideFctn: public AliHBTOnePairFctn1D -{ - - public: - AliHBTLLWeightTheorQSideFctn(Int_t nbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0); - virtual ~AliHBTLLWeightTheorQSideFctn(){} - - TH1* GetResult(); - void ProcessSameEventParticles(AliHBTPair* partpair); - - Double_t GetValue(AliHBTPair* partpair) - { return partpair->GetQSideCMSLC();} //isn't used - - ClassDef(AliHBTLLWeightTheorQSideFctn,1) -}; - -class AliHBTLLWeightTheorQLongFctn: public AliHBTOnePairFctn1D -{ - - public: - AliHBTLLWeightTheorQLongFctn(Int_t nbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0); - virtual ~AliHBTLLWeightTheorQLongFctn(){} - - TH1* GetResult(); - void ProcessSameEventParticles(AliHBTPair* partpair); - - Double_t GetValue(AliHBTPair* partpair) - { return partpair->GetQLongCMSLC();} //isn't used - - ClassDef(AliHBTLLWeightTheorQLongFctn,1) -}; - -#endif diff --git a/HBTAN/AliHBTLLWeights.cxx b/HBTAN/AliHBTLLWeights.cxx index 49e1868077f..32c19c7f266 100644 --- a/HBTAN/AliHBTLLWeights.cxx +++ b/HBTAN/AliHBTLLWeights.cxx @@ -197,7 +197,7 @@ AliHBTLLWeights::AliHBTLLWeights(): /**************************************************************/ AliHBTLLWeights::AliHBTLLWeights(const AliHBTLLWeights &/*source*/): - TObject(), + AliHBTWeights(), fTest(kTRUE), fColoumbSwitch(kTRUE), fQuantStatSwitch(kTRUE), @@ -240,6 +240,22 @@ AliHBTLLWeights* AliHBTLLWeights::Instance() } /************************************************************/ +void AliHBTLLWeights::Set() +{ + //sets this as weighitng class + Info("Set","Setting Lednicky-Lyuboshitz as Weighing Class"); + + if ( fgWeights == 0x0 ) + { + fgWeights = AliHBTLLWeights::Instance(); + return; + } + if ( fgWeights == AliHBTLLWeights::Instance() ) return; + delete fgWeights; + fgWeights = AliHBTLLWeights::Instance(); +} +/************************************************************/ + Double_t AliHBTLLWeights::GetWeight(const AliHBTPair* partpair) { // calculates weight for a pair diff --git a/HBTAN/AliHBTLLWeights.h b/HBTAN/AliHBTLLWeights.h index 6a6d6f665e7..b625d54178a 100644 --- a/HBTAN/AliHBTLLWeights.h +++ b/HBTAN/AliHBTLLWeights.h @@ -9,14 +9,15 @@ #ifndef ALIHBTLLWEIGHTS_H #define ALIHBTLLWEIGHTS_H -#include +#include "AliHBTWeights.h" class AliHBTPair; -class AliHBTLLWeights: public TObject +class AliHBTLLWeights: public AliHBTWeights { public: virtual ~AliHBTLLWeights(){;} static AliHBTLLWeights* Instance(); + void Set(); Double_t GetWeight(const AliHBTPair* partpair); //get weight calculated by Lednicky's algorithm diff --git a/HBTAN/AliHBTPair.cxx b/HBTAN/AliHBTPair.cxx index a385fb6edee..fe5c020bd9f 100644 --- a/HBTAN/AliHBTPair.cxx +++ b/HBTAN/AliHBTPair.cxx @@ -1,6 +1,6 @@ #include "AliHBTPair.h" #include "AliHBTParticle.h" -#include "AliHBTLLWeights.h" +#include "AliHBTWeights.h" ClassImp(AliHBTPair) @@ -33,8 +33,8 @@ AliHBTPair::AliHBTPair(Bool_t rev): fMassSqrNotCalc(kTRUE), fQInvL(0.0), fQInvLNotCalc(kTRUE), - fLLWeight(0.0), - ffLLWeightNotCalc(kTRUE), + fWeight(0.0), + fWeightNotCalc(kTRUE), fPxSum(0.0), fPySum(0.0), fPzSum(0.0), @@ -86,8 +86,8 @@ AliHBTPair::AliHBTPair(AliHBTParticle* part1, AliHBTParticle* part2, Bool_t rev) fMassSqrNotCalc(kTRUE), fQInvL(0.0), fQInvLNotCalc(kTRUE), - fLLWeight(0.0), - ffLLWeightNotCalc(kTRUE), + fWeight(0.0), + fWeightNotCalc(kTRUE), fPxSum(0.0), fPySum(0.0), fPzSum(0.0), @@ -305,12 +305,12 @@ Double_t AliHBTPair::GetMt() } /************************************************************************/ -Double_t AliHBTPair::GetLLWeight() +Double_t AliHBTPair::GetWeight() { - if (ffLLWeightNotCalc) + if (fWeightNotCalc) { - fLLWeight = AliHBTLLWeights::Instance()->GetWeight(this); - ffLLWeightNotCalc = kFALSE; + fWeight = AliHBTWeights::Weight(this); + fWeightNotCalc = kFALSE; } - return fLLWeight; + return fWeight; } diff --git a/HBTAN/AliHBTPair.h b/HBTAN/AliHBTPair.h index f01c5922f67..7c79d6752ea 100644 --- a/HBTAN/AliHBTPair.h +++ b/HBTAN/AliHBTPair.h @@ -41,7 +41,7 @@ class AliHBTPair: public TObject virtual Double_t GetDeltaPz(); virtual Double_t GetGammaToCMSLC(); - Double_t GetLLWeight(); + Double_t GetWeight(); protected: AliHBTParticle* fPart1; //pointer to first particle AliHBTParticle* fPart2; //pointer to second particle @@ -92,8 +92,8 @@ class AliHBTPair: public TObject Bool_t fQInvLNotCalc; void CalculateQInvL(); - Double_t fLLWeight; - Bool_t ffLLWeightNotCalc; + Double_t fWeight; + Bool_t fWeightNotCalc; Double_t fPxSum; Double_t fPySum; @@ -153,7 +153,7 @@ void AliHBTPair::Changed() fKStarNotCalc = kTRUE; fQInvLNotCalc = kTRUE; fGammaCMSLCNotCalc = kTRUE; - ffLLWeightNotCalc = kTRUE; + fWeightNotCalc = kTRUE; } /****************************************************************/ inline diff --git a/HBTAN/AliHBTParticle.h b/HBTAN/AliHBTParticle.h index 7f3f9d90d60..3fa2047b023 100644 --- a/HBTAN/AliHBTParticle.h +++ b/HBTAN/AliHBTParticle.h @@ -99,6 +99,9 @@ public: void SetProductionVertex(const TLorentzVector& v) {SetProductionVertex(v.X(),v.Y(),v.Z(),v.T());} void SetCalcMass(Double_t mass) {fCalcMass = mass;} + + void SetUID(Int_t id){fIdxInEvent = id;} + const Char_t* GetName() const; void Print() const; diff --git a/HBTAN/AliHBTReaderESD.cxx b/HBTAN/AliHBTReaderESD.cxx index 9ea9657f0fc..232b67e8232 100644 --- a/HBTAN/AliHBTReaderESD.cxx +++ b/HBTAN/AliHBTReaderESD.cxx @@ -198,7 +198,8 @@ Int_t AliHBTReaderESD::ReadNext() for (Int_t s = 0; sGetNumberOfParticles(), fTracksEvent->GetNumberOfParticles(), + fTracksEvent->GetNumberOfParticles(), fParticlesEvent->GetNumberOfParticles(), fNEventsRead,fCurrentEvent,fCurrentDir); fCurrentEvent++; diff --git a/HBTAN/AliHBTLLWeightFctn.cxx b/HBTAN/AliHBTWeightFctn.cxx similarity index 64% rename from HBTAN/AliHBTLLWeightFctn.cxx rename to HBTAN/AliHBTWeightFctn.cxx index 0991d5de201..2d2b91fd42e 100644 --- a/HBTAN/AliHBTLLWeightFctn.cxx +++ b/HBTAN/AliHBTWeightFctn.cxx @@ -1,4 +1,4 @@ -#include "AliHBTLLWeightFctn.h" +#include "AliHBTWeightFctn.h" /* $Id$ */ //_________________________________________________________________________ // @@ -12,15 +12,15 @@ //One needs both pairs //(simulated and recontructed), thus function is of class AliHBTTwoPairFctn1D. //Author: Ludmila Malinina, JINR (malinina@sunhe.jinr.ru) -#include "AliHBTLLWeightsPID.h" +#include "AliHBTWeightsPID.h" -//--for test--AliHBTLLWeightQInvFctn* yyy= new AliHBTLLWeightQInvFctn(); +//--for test--AliHBTWeightQInvFctn* yyy= new AliHBTWeightQInvFctn(); -ClassImp( AliHBTLLWeightQInvFctn ) +ClassImp( AliHBTWeightQInvFctn ) /****************************************************************/ -AliHBTLLWeightQInvFctn::AliHBTLLWeightQInvFctn(Int_t nbins, Double_t maxXval, Double_t minXval): +AliHBTWeightQInvFctn::AliHBTWeightQInvFctn(Int_t nbins, Double_t maxXval, Double_t minXval): AliHBTTwoPairFctn1D(nbins,maxXval,minXval) { //ctor @@ -28,22 +28,27 @@ AliHBTLLWeightQInvFctn::AliHBTLLWeightQInvFctn(Int_t nbins, Double_t maxXval, Do Rename("wqinvcf","Q_{inv} Weight Correlation Function"); } /****************************************************************/ -void AliHBTLLWeightQInvFctn::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) +void AliHBTWeightQInvFctn::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) { //process particles from same events (fills numerator) trackpair = CheckPair(trackpair); partpair = CheckPair(partpair); if ( trackpair && partpair) { - Double_t weightPID=1.; - Double_t weightHBT=partpair->GetLLWeight(); - Double_t weight=weightHBT*weightPID; - if(TMath::Abs(weight)<=10.)fNumerator->Fill(trackpair->GetQInv(),weight); +// Double_t weightPID=1.; + Double_t weight = 1.0; + if ( ( trackpair->Particle1()->GetPdgCode() == partpair->Particle1()->GetPdgCode()) && + ( trackpair->Particle2()->GetPdgCode() == partpair->Particle2()->GetPdgCode()) ) + { + weight=partpair->GetWeight(); + } +// Double_t weight=weightHBT*weightPID; + fNumerator->Fill(trackpair->GetQInv(),weight); } } /****************************************************************/ -void AliHBTLLWeightQInvFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) +void AliHBTWeightQInvFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) { // Fills the denominator using mixed pairs trackpair = CheckPair(trackpair); @@ -54,7 +59,7 @@ void AliHBTLLWeightQInvFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, A } } /**************************************************************/ -TH1* AliHBTLLWeightQInvFctn::GetResult() +TH1* AliHBTWeightQInvFctn::GetResult() { //returns ratio of numerator and denominator return GetRatio(Scale()); @@ -65,9 +70,9 @@ TH1* AliHBTLLWeightQInvFctn::GetResult() /**************************************************************************************/ /**************************************************************************************/ -ClassImp(AliHBTLLWeightQOutFctn) +ClassImp(AliHBTWeightQOutFctn) -AliHBTLLWeightQOutFctn::AliHBTLLWeightQOutFctn(Int_t nbins, Double_t maxXval, Double_t minXval): +AliHBTWeightQOutFctn::AliHBTWeightQOutFctn(Int_t nbins, Double_t maxXval, Double_t minXval): AliHBTTwoPairFctn1D(nbins,maxXval,minXval) { //ctor @@ -75,22 +80,27 @@ AliHBTLLWeightQOutFctn::AliHBTLLWeightQOutFctn(Int_t nbins, Double_t maxXval, Do Rename("wqoutcf","Q_{out} Weight Correlation Function"); } /****************************************************************/ -void AliHBTLLWeightQOutFctn::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) +void AliHBTWeightQOutFctn::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) { //process particles from same events (fills numerator) trackpair = CheckPair(trackpair); partpair = CheckPair(partpair); if ( trackpair && partpair) { - Double_t weightPID=1.; - Double_t weightHBT=partpair->GetLLWeight(); - Double_t weight=weightHBT*weightPID; - if(TMath::Abs(weight)<=10.) fNumerator->Fill(trackpair->GetQOutCMSLC(),weight); +// Double_t weightPID=1.; + Double_t weight = 1.0; + if ( ( trackpair->Particle1()->GetPdgCode() == partpair->Particle1()->GetPdgCode()) && + ( trackpair->Particle2()->GetPdgCode() == partpair->Particle2()->GetPdgCode()) ) + { + weight=partpair->GetWeight(); + } +// Double_t weight=weightHBT*weightPID; + fNumerator->Fill(trackpair->GetQOutCMSLC(),weight); } } /****************************************************************/ -void AliHBTLLWeightQOutFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) +void AliHBTWeightQOutFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) { //process particles from diff events (fills denominator) trackpair = CheckPair(trackpair); @@ -101,7 +111,7 @@ void AliHBTLLWeightQOutFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, Al } } /**************************************************************/ -TH1* AliHBTLLWeightQOutFctn::GetResult() +TH1* AliHBTWeightQOutFctn::GetResult() { //returns ratio of numerator and denominator @@ -112,8 +122,8 @@ TH1* AliHBTLLWeightQOutFctn::GetResult() /*************************************************************************************/ /*************************************************************************************/ -ClassImp(AliHBTLLWeightQLongFctn) -AliHBTLLWeightQLongFctn::AliHBTLLWeightQLongFctn(Int_t nbins, Double_t maxXval, Double_t minXval): +ClassImp(AliHBTWeightQLongFctn) +AliHBTWeightQLongFctn::AliHBTWeightQLongFctn(Int_t nbins, Double_t maxXval, Double_t minXval): AliHBTTwoPairFctn1D(nbins,maxXval,minXval) { //ctor @@ -121,22 +131,27 @@ AliHBTLLWeightQLongFctn::AliHBTLLWeightQLongFctn(Int_t nbins, Double_t maxXval, Rename("wqlongcf","Q_{long} Weight Correlation Function"); } /****************************************************************/ -void AliHBTLLWeightQLongFctn::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) +void AliHBTWeightQLongFctn::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) { //process particles from same events (fills numerator) trackpair = CheckPair(trackpair); partpair = CheckPair(partpair); if ( trackpair && partpair) { - Double_t weightPID=1.; - Double_t weightHBT=partpair->GetLLWeight(); - Double_t weight=weightHBT*weightPID; - if(TMath::Abs(weight)<=10.) fNumerator->Fill(trackpair->GetQLongCMSLC(),weight); +// Double_t weightPID=1.; + Double_t weight = 1.0; + if ( ( trackpair->Particle1()->GetPdgCode() == partpair->Particle1()->GetPdgCode()) && + ( trackpair->Particle2()->GetPdgCode() == partpair->Particle2()->GetPdgCode()) ) + { + weight=partpair->GetWeight(); + } +// Double_t weight=weightHBT*weightPID; + fNumerator->Fill(trackpair->GetQLongCMSLC(),weight); } } /****************************************************************/ -void AliHBTLLWeightQLongFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) +void AliHBTWeightQLongFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) { //process particles from diff events (fills denominator) trackpair = CheckPair(trackpair); @@ -147,7 +162,7 @@ void AliHBTLLWeightQLongFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, A } } /**************************************************************/ -TH1* AliHBTLLWeightQLongFctn::GetResult() +TH1* AliHBTWeightQLongFctn::GetResult() { //returns ratio of numerator and denominator @@ -158,10 +173,10 @@ TH1* AliHBTLLWeightQLongFctn::GetResult() /*************************************************************************************/ /*************************************************************************************/ -ClassImp(AliHBTLLWeightQSideFctn) +ClassImp(AliHBTWeightQSideFctn) /*************************************************************************************/ -AliHBTLLWeightQSideFctn::AliHBTLLWeightQSideFctn(Int_t nbins, Double_t maxXval, Double_t minXval): +AliHBTWeightQSideFctn::AliHBTWeightQSideFctn(Int_t nbins, Double_t maxXval, Double_t minXval): AliHBTTwoPairFctn1D(nbins,maxXval,minXval) { //ctor @@ -169,22 +184,27 @@ AliHBTLLWeightQSideFctn::AliHBTLLWeightQSideFctn(Int_t nbins, Double_t maxXval, Rename("wqsidecf","Q_{side} Weight Correlation Function"); } /****************************************************************/ -void AliHBTLLWeightQSideFctn::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) +void AliHBTWeightQSideFctn::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) { //process particles from same events (fills numerator) trackpair = CheckPair(trackpair); partpair = CheckPair(partpair); if ( trackpair && partpair) { - Double_t weightPID=1.; - Double_t weightHBT=partpair->GetLLWeight(); - Double_t weight=weightHBT*weightPID; - if(TMath::Abs(weight)<=10.) fNumerator->Fill(trackpair->GetQSideCMSLC(),weight); +// Double_t weightPID=1.; + Double_t weight = 1.0; + if ( ( trackpair->Particle1()->GetPdgCode() == partpair->Particle1()->GetPdgCode()) && + ( trackpair->Particle2()->GetPdgCode() == partpair->Particle2()->GetPdgCode()) ) + { + weight=partpair->GetWeight(); + } +// Double_t weight=weightHBT*weightPID; + fNumerator->Fill(trackpair->GetQSideCMSLC(),weight); } } /****************************************************************/ -void AliHBTLLWeightQSideFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) +void AliHBTWeightQSideFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) { //process particles from diff events (fills denominator) trackpair = CheckPair(trackpair); @@ -195,7 +215,7 @@ void AliHBTLLWeightQSideFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, } } /**************************************************************/ -TH1* AliHBTLLWeightQSideFctn::GetResult() +TH1* AliHBTWeightQSideFctn::GetResult() { //returns ratio of numerator and denominator @@ -206,9 +226,9 @@ TH1* AliHBTLLWeightQSideFctn::GetResult() /*************************************************************************************/ /*************************************************************************************/ -ClassImp(AliHBTLLWeightTwoKStarFctn) +ClassImp(AliHBTWeightTwoKStarFctn) /*************************************************************************************/ -AliHBTLLWeightTwoKStarFctn::AliHBTLLWeightTwoKStarFctn(Int_t nbins, Double_t maxXval, Double_t minXval): +AliHBTWeightTwoKStarFctn::AliHBTWeightTwoKStarFctn(Int_t nbins, Double_t maxXval, Double_t minXval): AliHBTTwoPairFctn1D(nbins,maxXval,minXval) { //ctor @@ -216,22 +236,27 @@ AliHBTLLWeightTwoKStarFctn::AliHBTLLWeightTwoKStarFctn(Int_t nbins, Double_t max Rename("wtwokstarcf","2*K^{*} Weight Correlation Function"); } /****************************************************************/ -void AliHBTLLWeightTwoKStarFctn::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) +void AliHBTWeightTwoKStarFctn::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) { //process particles from same events (fills numerator) trackpair = CheckPair(trackpair); partpair = CheckPair(partpair); if ( trackpair && partpair) { - Double_t weightPID=1.; - Double_t weightHBT=partpair->GetLLWeight(); - Double_t weight=weightHBT*weightPID; - if(TMath::Abs(weight)<=10.) fNumerator->Fill(2.0*(trackpair->GetKStar()),weight); +// Double_t weightPID=1.; + Double_t weight = 1.0; + if ( ( trackpair->Particle1()->GetPdgCode() == partpair->Particle1()->GetPdgCode()) && + ( trackpair->Particle2()->GetPdgCode() == partpair->Particle2()->GetPdgCode()) ) + { + weight=partpair->GetWeight(); + } +// Double_t weight=weightHBT*weightPID; + fNumerator->Fill(2.0*(trackpair->GetKStar()),weight); } } /****************************************************************/ -void AliHBTLLWeightTwoKStarFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) +void AliHBTWeightTwoKStarFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) { //process particles from diff events (fills denominator) trackpair = CheckPair(trackpair); @@ -242,7 +267,7 @@ void AliHBTLLWeightTwoKStarFctn::ProcessDiffEventParticles(AliHBTPair* trackpai } } /**************************************************************/ -TH1* AliHBTLLWeightTwoKStarFctn::GetResult() +TH1* AliHBTWeightTwoKStarFctn::GetResult() { //returns ratio of numerator and denominator @@ -253,10 +278,10 @@ TH1* AliHBTLLWeightTwoKStarFctn::GetResult() /*************************************************************************************/ /*************************************************************************************/ -ClassImp(AliHBTLLWeightQOutQSideFctn) +ClassImp(AliHBTWeightQOutQSideFctn) /*************************************************************************************/ -AliHBTLLWeightQOutQSideFctn::AliHBTLLWeightQOutQSideFctn(Int_t nxbins, Double_t maxXval, Double_t minXval, +AliHBTWeightQOutQSideFctn::AliHBTWeightQOutQSideFctn(Int_t nxbins, Double_t maxXval, Double_t minXval, Int_t nybins, Double_t maxYval, Double_t minYval): AliHBTTwoPairFctn2D(nxbins,maxXval,minXval,nybins,maxYval,minYval) { @@ -265,23 +290,27 @@ AliHBTLLWeightQOutQSideFctn::AliHBTLLWeightQOutQSideFctn(Int_t nxbins, Double_t Rename("wqoutqsidecf","Q_{out} Q_{side} Weight Correlation Function 2D"); } /*************************************************************************************/ -void AliHBTLLWeightQOutQSideFctn::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) +void AliHBTWeightQOutQSideFctn::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) { //process particles from same events (fills numerator) trackpair = CheckPair(trackpair); partpair = CheckPair(partpair); if ( trackpair && partpair) { - Double_t weightPID=1.; - Double_t weightHBT=partpair->GetLLWeight(); - Double_t weight=weightHBT*weightPID; - if(TMath::Abs(weight)<=10.) - fNumerator->Fill(trackpair->GetQOutCMSLC(),trackpair->GetQSideCMSLC(),weight); +// Double_t weightPID=1.; + Double_t weight = 1.0; + if ( ( trackpair->Particle1()->GetPdgCode() == partpair->Particle1()->GetPdgCode()) && + ( trackpair->Particle2()->GetPdgCode() == partpair->Particle2()->GetPdgCode()) ) + { + weight=partpair->GetWeight(); + } +// Double_t weight=weightHBT*weightPID; + fNumerator->Fill(trackpair->GetQOutCMSLC(),trackpair->GetQSideCMSLC(),weight); } } /****************************************************************/ -void AliHBTLLWeightQOutQSideFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) +void AliHBTWeightQOutQSideFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) { //process particles from diff events (fills denominator) trackpair = CheckPair(trackpair); @@ -292,7 +321,7 @@ void AliHBTLLWeightQOutQSideFctn::ProcessDiffEventParticles(AliHBTPair* trackpai } } /**************************************************************/ -TH1* AliHBTLLWeightQOutQSideFctn::GetResult() +TH1* AliHBTWeightQOutQSideFctn::GetResult() { //returns result return GetRatio(Scale()); @@ -302,10 +331,10 @@ TH1* AliHBTLLWeightQOutQSideFctn::GetResult() /*************************************************************************************/ /*************************************************************************************/ -ClassImp(AliHBTLLWeightQOutQLongFctn) +ClassImp(AliHBTWeightQOutQLongFctn) /*************************************************************************************/ -AliHBTLLWeightQOutQLongFctn::AliHBTLLWeightQOutQLongFctn(Int_t nxbins, Double_t maxXval, Double_t minXval, +AliHBTWeightQOutQLongFctn::AliHBTWeightQOutQLongFctn(Int_t nxbins, Double_t maxXval, Double_t minXval, Int_t nybins, Double_t maxYval, Double_t minYval): AliHBTTwoPairFctn2D(nxbins,maxXval,minXval,nybins,maxYval,minYval) { @@ -314,23 +343,27 @@ AliHBTLLWeightQOutQLongFctn::AliHBTLLWeightQOutQLongFctn(Int_t nxbins, Double_t Rename("wqoutqlongcf","Q_{out} Q_{long} Weight Correlation Function 2D"); } /*************************************************************************************/ -void AliHBTLLWeightQOutQLongFctn::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) +void AliHBTWeightQOutQLongFctn::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) { //process particles from same events (fills numerator) trackpair = CheckPair(trackpair); partpair = CheckPair(partpair); if ( trackpair && partpair) { - Double_t weightPID=1.; - Double_t weightHBT=partpair->GetLLWeight(); - Double_t weight=weightHBT*weightPID; - if(TMath::Abs(weight)<=10.) - fNumerator->Fill(trackpair->GetQOutCMSLC(),trackpair->GetQLongCMSLC(),weight); +// Double_t weightPID=1.; + Double_t weight = 1.0; + if ( ( trackpair->Particle1()->GetPdgCode() == partpair->Particle1()->GetPdgCode()) && + ( trackpair->Particle2()->GetPdgCode() == partpair->Particle2()->GetPdgCode()) ) + { + weight=partpair->GetWeight(); + } +// Double_t weight=weightHBT*weightPID; + fNumerator->Fill(trackpair->GetQOutCMSLC(),trackpair->GetQLongCMSLC(),weight); } } /****************************************************************/ -void AliHBTLLWeightQOutQLongFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) +void AliHBTWeightQOutQLongFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) { //process particles from diff events (fills denominator) trackpair = CheckPair(trackpair); @@ -342,7 +375,7 @@ void AliHBTLLWeightQOutQLongFctn::ProcessDiffEventParticles(AliHBTPair* trackpai } /**************************************************************/ -TH1* AliHBTLLWeightQOutQLongFctn::GetResult() +TH1* AliHBTWeightQOutQLongFctn::GetResult() { //returns result return GetRatio(Scale()); @@ -352,10 +385,10 @@ TH1* AliHBTLLWeightQOutQLongFctn::GetResult() /*************************************************************************************/ /*************************************************************************************/ -ClassImp(AliHBTLLWeightQSideQLongFctn) +ClassImp(AliHBTWeightQSideQLongFctn) /*************************************************************************************/ -AliHBTLLWeightQSideQLongFctn::AliHBTLLWeightQSideQLongFctn(Int_t nxbins, Double_t maxXval, Double_t minXval, +AliHBTWeightQSideQLongFctn::AliHBTWeightQSideQLongFctn(Int_t nxbins, Double_t maxXval, Double_t minXval, Int_t nybins, Double_t maxYval, Double_t minYval): AliHBTTwoPairFctn2D(nxbins,maxXval,minXval,nybins,maxYval,minYval) { @@ -364,23 +397,27 @@ AliHBTLLWeightQSideQLongFctn::AliHBTLLWeightQSideQLongFctn(Int_t nxbins, Double_ Rename("wqsideqlongcf","Q_{side} Q_{long} Weight Correlation Function 2D"); } /*************************************************************************************/ -void AliHBTLLWeightQSideQLongFctn::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) +void AliHBTWeightQSideQLongFctn::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) { //process particles from same events (fills numerator) trackpair = CheckPair(trackpair); partpair = CheckPair(partpair); if ( trackpair && partpair) { - Double_t weightPID=1.; - Double_t weightHBT=partpair->GetLLWeight(); - Double_t weight=weightHBT*weightPID; - if(TMath::Abs(weight)<=10.) - fNumerator->Fill(trackpair->GetQSideCMSLC(),trackpair->GetQLongCMSLC(),weight); +// Double_t weightPID=1.; + Double_t weight = 1.0; + if ( ( trackpair->Particle1()->GetPdgCode() == partpair->Particle1()->GetPdgCode()) && + ( trackpair->Particle2()->GetPdgCode() == partpair->Particle2()->GetPdgCode()) ) + { + weight=partpair->GetWeight(); + } +// Double_t weight=weightHBT*weightPID; + fNumerator->Fill(trackpair->GetQSideCMSLC(),trackpair->GetQLongCMSLC(),weight); } } /****************************************************************/ -void AliHBTLLWeightQSideQLongFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) +void AliHBTWeightQSideQLongFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) { //process particles from diff events (fills denominator) trackpair = CheckPair(trackpair); @@ -391,7 +428,7 @@ void AliHBTLLWeightQSideQLongFctn::ProcessDiffEventParticles(AliHBTPair* trackpa } } /**************************************************************/ -TH1* AliHBTLLWeightQSideQLongFctn::GetResult() +TH1* AliHBTWeightQSideQLongFctn::GetResult() { //returns result return GetRatio(Scale()); diff --git a/HBTAN/AliHBTLLWeightFctn.h b/HBTAN/AliHBTWeightFctn.h similarity index 67% rename from HBTAN/AliHBTLLWeightFctn.h rename to HBTAN/AliHBTWeightFctn.h index 97c00d0c45f..278eb5c115e 100644 --- a/HBTAN/AliHBTLLWeightFctn.h +++ b/HBTAN/AliHBTWeightFctn.h @@ -1,5 +1,5 @@ -#ifndef ALIHBTLLWEIGHTQINVFCTN_H -#define ALIHBTLLWEIGHTQINVFCTN_H +#ifndef ALIHBTWeightQINVFCTN_H +#define ALIHBTWeightQINVFCTN_H /* $Id$ */ @@ -7,12 +7,12 @@ #include "AliHBTFunction.h" -class AliHBTLLWeights; -class AliHBTLLWeightQInvFctn: public AliHBTTwoPairFctn1D +class AliHBTWeights; +class AliHBTWeightQInvFctn: public AliHBTTwoPairFctn1D { public: - AliHBTLLWeightQInvFctn(Int_t nbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0); - virtual ~AliHBTLLWeightQInvFctn(){}; + AliHBTWeightQInvFctn(Int_t nbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0); + virtual ~AliHBTWeightQInvFctn(){}; TH1* GetResult(); void ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair); @@ -21,17 +21,17 @@ class AliHBTLLWeightQInvFctn: public AliHBTTwoPairFctn1D protected: Double_t GetValue(AliHBTPair* trackpair, AliHBTPair* partpair) { return trackpair->GetQInv()-partpair->GetQInv();} //isn't use - ClassDef(AliHBTLLWeightQInvFctn,1) + ClassDef(AliHBTWeightQInvFctn,1) }; /*************************************************************************************/ -class AliHBTLLWeightQOutFctn: public AliHBTTwoPairFctn1D +class AliHBTWeightQOutFctn: public AliHBTTwoPairFctn1D { // friend class AliHBTOnePairFctn1D; public: - AliHBTLLWeightQOutFctn(Int_t nbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0); - virtual ~AliHBTLLWeightQOutFctn(){}; + AliHBTWeightQOutFctn(Int_t nbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0); + virtual ~AliHBTWeightQOutFctn(){}; TH1* GetResult(); void ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair); void ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair); @@ -39,17 +39,17 @@ class AliHBTLLWeightQOutFctn: public AliHBTTwoPairFctn1D protected: Double_t GetValue(AliHBTPair* trackpair, AliHBTPair* partpair) { return trackpair->GetQOutCMSLC()-partpair->GetQOutCMSLC();} //isn't use - ClassDef(AliHBTLLWeightQOutFctn,1) + ClassDef(AliHBTWeightQOutFctn,1) }; /*************************************************************************************/ -class AliHBTLLWeightQLongFctn: public AliHBTTwoPairFctn1D +class AliHBTWeightQLongFctn: public AliHBTTwoPairFctn1D { // friend class AliHBTOnePairFctn1D; public: - AliHBTLLWeightQLongFctn(Int_t nbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0); - virtual ~AliHBTLLWeightQLongFctn(){}; + AliHBTWeightQLongFctn(Int_t nbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0); + virtual ~AliHBTWeightQLongFctn(){}; TH1* GetResult(); void ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair); void ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair); @@ -58,17 +58,17 @@ class AliHBTLLWeightQLongFctn: public AliHBTTwoPairFctn1D Double_t GetValue(AliHBTPair* trackpair, AliHBTPair* partpair) { return trackpair->GetQLongCMSLC()-partpair->GetQLongCMSLC();} //isn't used - ClassDef(AliHBTLLWeightQLongFctn,1) + ClassDef(AliHBTWeightQLongFctn,1) }; /*************************************************************************************/ -class AliHBTLLWeightQSideFctn: public AliHBTTwoPairFctn1D +class AliHBTWeightQSideFctn: public AliHBTTwoPairFctn1D { // friend class AliHBTOnePairFctn1D; public: - AliHBTLLWeightQSideFctn(Int_t nbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0); - virtual ~AliHBTLLWeightQSideFctn(){}; + AliHBTWeightQSideFctn(Int_t nbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0); + virtual ~AliHBTWeightQSideFctn(){}; TH1* GetResult(); void ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair); void ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair); @@ -77,16 +77,16 @@ class AliHBTLLWeightQSideFctn: public AliHBTTwoPairFctn1D Double_t GetValue(AliHBTPair* trackpair, AliHBTPair* partpair) { return trackpair->GetQLongCMSLC()-partpair->GetQLongCMSLC();} //isn't used - ClassDef(AliHBTLLWeightQSideFctn,1) + ClassDef(AliHBTWeightQSideFctn,1) }; /*************************************************************************************/ -class AliHBTLLWeightTwoKStarFctn: public AliHBTTwoPairFctn1D +class AliHBTWeightTwoKStarFctn: public AliHBTTwoPairFctn1D { // friend class AliHBTOnePairFctn1D; public: - AliHBTLLWeightTwoKStarFctn(Int_t nbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0); - virtual ~AliHBTLLWeightTwoKStarFctn(){}; + AliHBTWeightTwoKStarFctn(Int_t nbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0); + virtual ~AliHBTWeightTwoKStarFctn(){}; TH1* GetResult(); void ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair); void ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair); @@ -94,65 +94,65 @@ class AliHBTLLWeightTwoKStarFctn: public AliHBTTwoPairFctn1D protected: Double_t GetValue(AliHBTPair* trackpair, AliHBTPair* partpair) { return trackpair->GetKStar()-partpair->GetKStar();} //isn't used - ClassDef(AliHBTLLWeightTwoKStarFctn,1) + ClassDef(AliHBTWeightTwoKStarFctn,1) }; /*************************************************************************************/ -class AliHBTLLWeightQOutQSideFctn: public AliHBTTwoPairFctn2D +class AliHBTWeightQOutQSideFctn: public AliHBTTwoPairFctn2D { // friend class AliHBTOnePairFctn1D; public: - AliHBTLLWeightQOutQSideFctn(Int_t nxbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0, + AliHBTWeightQOutQSideFctn(Int_t nxbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0, Int_t nybins = 100, Double_t maxYval = 0.15, Double_t minYval = 0.0); - virtual ~AliHBTLLWeightQOutQSideFctn(){}; + virtual ~AliHBTWeightQOutQSideFctn(){}; TH1* GetResult(); void ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair); void ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair); protected: void GetValues(AliHBTPair* /*trackpair*/, AliHBTPair* /*partpair*/, Double_t& /*x*/, Double_t& /*y*/){} - ClassDef(AliHBTLLWeightQOutQSideFctn,1) + ClassDef(AliHBTWeightQOutQSideFctn,1) }; /*************************************************************************************/ -class AliHBTLLWeightQOutQLongFctn: public AliHBTTwoPairFctn2D +class AliHBTWeightQOutQLongFctn: public AliHBTTwoPairFctn2D { // friend class AliHBTOnePairFctn1D; public: - AliHBTLLWeightQOutQLongFctn(Int_t nxbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0, + AliHBTWeightQOutQLongFctn(Int_t nxbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0, Int_t nybins = 100, Double_t maxYval = 0.15, Double_t minYval = 0.0); - virtual ~AliHBTLLWeightQOutQLongFctn(){}; + virtual ~AliHBTWeightQOutQLongFctn(){}; TH1* GetResult(); void ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair); void ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair); protected: void GetValues(AliHBTPair* /*trackpair*/, AliHBTPair* /*partpair*/, Double_t& /*x*/, Double_t& /*y*/){} - ClassDef(AliHBTLLWeightQOutQLongFctn,1) + ClassDef(AliHBTWeightQOutQLongFctn,1) }; /*************************************************************************************/ -class AliHBTLLWeightQSideQLongFctn: public AliHBTTwoPairFctn2D +class AliHBTWeightQSideQLongFctn: public AliHBTTwoPairFctn2D { // friend class AliHBTOnePairFctn1D; public: - AliHBTLLWeightQSideQLongFctn(Int_t nxbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0, + AliHBTWeightQSideQLongFctn(Int_t nxbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0, Int_t nybins = 100, Double_t maxYval = 0.15, Double_t minYval = 0.0); - virtual ~AliHBTLLWeightQSideQLongFctn(){}; + virtual ~AliHBTWeightQSideQLongFctn(){}; TH1* GetResult(); void ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair); void ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair); protected: void GetValues(AliHBTPair* /*trackpair*/, AliHBTPair* /*partpair*/, Double_t& /*x*/, Double_t& /*y*/){} - ClassDef(AliHBTLLWeightQSideQLongFctn,1) + ClassDef(AliHBTWeightQSideQLongFctn,1) }; diff --git a/HBTAN/AliHBTLLWeightTheorFctn.cxx b/HBTAN/AliHBTWeightTheorFctn.cxx similarity index 57% rename from HBTAN/AliHBTLLWeightTheorFctn.cxx rename to HBTAN/AliHBTWeightTheorFctn.cxx index 515ab870f97..5bd7a34a9da 100644 --- a/HBTAN/AliHBTLLWeightTheorFctn.cxx +++ b/HBTAN/AliHBTWeightTheorFctn.cxx @@ -7,18 +7,17 @@ // Author: Ludmila Malinina, JINR (malinina@sunhe.jinr.ru) //----------------------------------------------------------- -#include "AliHBTLLWeightTheorFctn.h" +#include "AliHBTWeightTheorFctn.h" -//--for test--AliHBTLLWeightQInvFctn* yyy= new AliHBTLLWeightQInvFctn(); +//--for test--AliHBTWeightQInvFctn* yyy= new AliHBTWeightQInvFctn(); /*************************************************************/ /*************************************************************/ /*************************************************************/ -ClassImp(AliHBTLLWeightTheorQInvFctn) +ClassImp(AliHBTWeightTheorQInvFctn) /*************************************************************/ -AliHBTLLWeightTheorQInvFctn:: -AliHBTLLWeightTheorQInvFctn(Int_t nbins, Double_t maxXval, Double_t minXval): +AliHBTWeightTheorQInvFctn::AliHBTWeightTheorQInvFctn(Int_t nbins, Double_t maxXval, Double_t minXval): AliHBTOnePairFctn1D(nbins,maxXval,minXval) { //ctor @@ -26,17 +25,17 @@ AliHBTLLWeightTheorQInvFctn(Int_t nbins, Double_t maxXval, Double_t minXval): Rename("wqinvtheorcf","Q_{inv} Weight Theoretical Correlation Function"); } /****************************************************************/ -void AliHBTLLWeightTheorQInvFctn::ProcessSameEventParticles(AliHBTPair* partpair) +void AliHBTWeightTheorQInvFctn::ProcessSameEventParticles(AliHBTPair* partpair) { - //Processes Particles and tracks Same different even + //Processes Particles and tracks Same different event partpair = CheckPair(partpair); if (partpair == 0x0) return; - Double_t weight = partpair->GetLLWeight(); + Double_t weight = partpair->GetWeight(); if(TMath::Abs(weight)<=10.) fNumerator->Fill(partpair->GetQInv(),weight); } /**************************************************************/ -TH1* AliHBTLLWeightTheorQInvFctn::GetResult() +TH1* AliHBTWeightTheorQInvFctn::GetResult() { //returns ratio of numerator and denominator return GetRatio(Scale()); @@ -46,11 +45,10 @@ TH1* AliHBTLLWeightTheorQInvFctn::GetResult() /*************************************************************/ /*************************************************************/ -ClassImp(AliHBTLLWeightTheorQOutFctn) +ClassImp(AliHBTWeightTheorQOutFctn) /*************************************************************/ -AliHBTLLWeightTheorQOutFctn:: -AliHBTLLWeightTheorQOutFctn(Int_t nbins, Double_t maxXval, Double_t minXval): +AliHBTWeightTheorQOutFctn::AliHBTWeightTheorQOutFctn(Int_t nbins, Double_t maxXval, Double_t minXval): AliHBTOnePairFctn1D(nbins,maxXval,minXval) { //ctor @@ -58,17 +56,17 @@ AliHBTLLWeightTheorQOutFctn(Int_t nbins, Double_t maxXval, Double_t minXval): Rename("wqouttheorcf","Q_{out} Weight Theoretical Correlation Function"); } /****************************************************************/ -void AliHBTLLWeightTheorQOutFctn::ProcessSameEventParticles(AliHBTPair* partpair) +void AliHBTWeightTheorQOutFctn::ProcessSameEventParticles(AliHBTPair* partpair) { //Processes Particles and tracks Same different even partpair = CheckPair(partpair); if (partpair == 0x0) return; - Double_t weight = partpair->GetLLWeight(); + Double_t weight = partpair->GetWeight(); if(TMath::Abs(weight)<=10.) fNumerator->Fill(partpair->GetQOutCMSLC(),weight); } /**************************************************************/ -TH1* AliHBTLLWeightTheorQOutFctn::GetResult() +TH1* AliHBTWeightTheorQOutFctn::GetResult() { //returns ratio of numerator and denominator return GetRatio(Scale()); @@ -78,11 +76,10 @@ TH1* AliHBTLLWeightTheorQOutFctn::GetResult() /*************************************************************/ /*************************************************************/ -ClassImp(AliHBTLLWeightTheorQSideFctn) +ClassImp(AliHBTWeightTheorQSideFctn) /*************************************************************/ -AliHBTLLWeightTheorQSideFctn:: -AliHBTLLWeightTheorQSideFctn(Int_t nbins, Double_t maxXval, Double_t minXval): +AliHBTWeightTheorQSideFctn::AliHBTWeightTheorQSideFctn(Int_t nbins, Double_t maxXval, Double_t minXval): AliHBTOnePairFctn1D(nbins,maxXval,minXval) { //ctor @@ -90,17 +87,17 @@ AliHBTLLWeightTheorQSideFctn(Int_t nbins, Double_t maxXval, Double_t minXval): Rename("wqsidetheorcf","Q_{side} Weight Theoretical Correlation Function"); } /****************************************************************/ -void AliHBTLLWeightTheorQSideFctn::ProcessSameEventParticles(AliHBTPair* partpair) +void AliHBTWeightTheorQSideFctn::ProcessSameEventParticles(AliHBTPair* partpair) { //Processes Particles and tracks Same different even partpair = CheckPair(partpair); if (partpair == 0x0) return; - Double_t weight = partpair->GetLLWeight(); + Double_t weight = partpair->GetWeight(); if(TMath::Abs(weight)<=10.) fNumerator->Fill(partpair->GetQSideCMSLC(),weight); } /**************************************************************/ -TH1* AliHBTLLWeightTheorQSideFctn::GetResult() +TH1* AliHBTWeightTheorQSideFctn::GetResult() { //returns ratio of numerator and denominator return GetRatio(Scale()); @@ -110,11 +107,10 @@ TH1* AliHBTLLWeightTheorQSideFctn::GetResult() /*************************************************************/ /*************************************************************/ -ClassImp(AliHBTLLWeightTheorQLongFctn) +ClassImp(AliHBTWeightTheorQLongFctn) /*************************************************************/ -AliHBTLLWeightTheorQLongFctn:: -AliHBTLLWeightTheorQLongFctn(Int_t nbins, Double_t maxXval, Double_t minXval): +AliHBTWeightTheorQLongFctn::AliHBTWeightTheorQLongFctn(Int_t nbins, Double_t maxXval, Double_t minXval): AliHBTOnePairFctn1D(nbins,maxXval,minXval) { //ctor @@ -122,17 +118,52 @@ AliHBTLLWeightTheorQLongFctn(Int_t nbins, Double_t maxXval, Double_t minXval): Rename("wqlongtheorcf","Q_{long} Weight Theoretical Correlation Function"); } /****************************************************************/ -void AliHBTLLWeightTheorQLongFctn::ProcessSameEventParticles(AliHBTPair* partpair) +void AliHBTWeightTheorQLongFctn::ProcessSameEventParticles(AliHBTPair* partpair) { //Processes Particles and tracks Same different even partpair = CheckPair(partpair); if (partpair == 0x0) return; - Double_t weight = partpair->GetLLWeight(); + Double_t weight = partpair->GetWeight(); if(TMath::Abs(weight)<=10.) fNumerator->Fill(partpair->GetQLongCMSLC(),weight); } /**************************************************************/ -TH1* AliHBTLLWeightTheorQLongFctn::GetResult() +TH1* AliHBTWeightTheorQLongFctn::GetResult() +{ + //returns ratio of numerator and denominator + return GetRatio(Scale()); +} + +/*************************************************************/ +/*************************************************************/ +/*************************************************************/ + +ClassImp(AliHBTWeightTheorOSLFctn) + +AliHBTWeightTheorOSLFctn::AliHBTWeightTheorOSLFctn(Int_t nXbins, Double_t maxXval, Double_t minXval, + Int_t nYbins, Double_t maxYval, Double_t minYval, + Int_t nZbins, Double_t maxZval, Double_t minZval): + AliHBTOnePairFctn3D(nXbins,maxXval,minXval,nYbins,maxYval,minYval,nZbins,maxZval,minZval) +{ + fWriteNumAndDen = kTRUE;//change default behaviour + Rename("wqosltheorcf","Q_{out}-Q_{side}-Q_{long} Weight Theoretical Correlation Fctn"); +} +/*************************************************************/ + +void AliHBTWeightTheorOSLFctn::ProcessSameEventParticles(AliHBTPair* partpair) +{ +//Fills numerator + partpair = CheckPair(partpair); + if (partpair == 0x0) return; + Double_t weight = partpair->GetWeight(); + Double_t out = TMath::Abs(partpair->GetQOutCMSLC()); + Double_t side = TMath::Abs(partpair->GetQSideCMSLC()); + Double_t lon = TMath::Abs(partpair->GetQLongCMSLC()); + fNumerator->Fill(out,side,lon,weight); +} +/*************************************************************/ + +TH1* AliHBTWeightTheorOSLFctn::GetResult() { //returns ratio of numerator and denominator return GetRatio(Scale()); diff --git a/HBTAN/AliHBTWeightTheorFctn.h b/HBTAN/AliHBTWeightTheorFctn.h new file mode 100644 index 00000000000..6ec31844eb9 --- /dev/null +++ b/HBTAN/AliHBTWeightTheorFctn.h @@ -0,0 +1,100 @@ +#ifndef ALIHBTWeightTHEORFCTN_H +#define ALIHBTWeightTHEORFCTN_H +/* $Id$ */ + +//This function allows to obtain Q_inv correlation function with weights +//calculated by Lednicky's alghorithm. +//Numerator is filled with weighted events. Weights are attributed to simulated particles. +//Weights are calculated with corresponding simulated particles momenta. +//Denominator is filled with mixing unweighted simulated particles. +//One needs only simulated pairs, so +//this function is of class AliHBTOnePairFctn1D. +//Author Ludmila Malinina JINR (malinina@sunhe.jinr.ru) + +#include "AliHBTFunction.h" + +class AliHBTWeights; + +class AliHBTWeightTheorQInvFctn: public AliHBTOnePairFctn1D +{ + + public: + AliHBTWeightTheorQInvFctn(Int_t nbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0); + virtual ~AliHBTWeightTheorQInvFctn(){} + + TH1* GetResult(); + void ProcessSameEventParticles(AliHBTPair* partpair); + + Double_t GetValue(AliHBTPair* partpair) + { return partpair->GetQInv();} + + ClassDef(AliHBTWeightTheorQInvFctn,1) +}; + +class AliHBTWeightTheorQOutFctn: public AliHBTOnePairFctn1D +{ + + public: + AliHBTWeightTheorQOutFctn(Int_t nbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0); + virtual ~AliHBTWeightTheorQOutFctn(){} + + TH1* GetResult(); + void ProcessSameEventParticles(AliHBTPair* partpair); + + Double_t GetValue(AliHBTPair* partpair) + { return partpair->GetQOutCMSLC();} + + ClassDef(AliHBTWeightTheorQOutFctn,1) +}; + +class AliHBTWeightTheorQSideFctn: public AliHBTOnePairFctn1D +{ + + public: + AliHBTWeightTheorQSideFctn(Int_t nbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0); + virtual ~AliHBTWeightTheorQSideFctn(){} + + TH1* GetResult(); + void ProcessSameEventParticles(AliHBTPair* partpair); + + Double_t GetValue(AliHBTPair* partpair) + { return partpair->GetQSideCMSLC();} + + ClassDef(AliHBTWeightTheorQSideFctn,1) +}; + +class AliHBTWeightTheorQLongFctn: public AliHBTOnePairFctn1D +{ + + public: + AliHBTWeightTheorQLongFctn(Int_t nbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0); + virtual ~AliHBTWeightTheorQLongFctn(){} + + TH1* GetResult(); + void ProcessSameEventParticles(AliHBTPair* partpair); + + Double_t GetValue(AliHBTPair* partpair) + { return partpair->GetQLongCMSLC();} + + ClassDef(AliHBTWeightTheorQLongFctn,1) +}; + +class AliHBTWeightTheorOSLFctn: public AliHBTOnePairFctn3D +{ + + public: + AliHBTWeightTheorOSLFctn(Int_t nXbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0, + Int_t nYbins = 100, Double_t maxYval = 0.15, Double_t minYval = 0.0, + Int_t nZbins = 100, Double_t maxZval = 0.15, Double_t minZval = 0.0); + virtual ~AliHBTWeightTheorOSLFctn(){} + + TH1* GetResult(); + void ProcessSameEventParticles(AliHBTPair* partpair); + + void GetValues(AliHBTPair* pair, Double_t& x, Double_t& y, Double_t& z) + { x=TMath::Abs(pair->GetQOutCMSLC()); y=TMath::Abs(pair->GetQSideCMSLC()); z=TMath::Abs(pair->GetQLongCMSLC());} + + ClassDef(AliHBTWeightTheorOSLFctn,1) +}; + +#endif diff --git a/HBTAN/AliHBTWeights.cxx b/HBTAN/AliHBTWeights.cxx new file mode 100644 index 00000000000..9952c118f06 --- /dev/null +++ b/HBTAN/AliHBTWeights.cxx @@ -0,0 +1,11 @@ +#include "AliHBTWeights.h" + +ClassImp(AliHBTWeights) + +AliHBTWeights* AliHBTWeights::fgWeights = 0x0; + +AliHBTWeights::~AliHBTWeights() + { + delete fgWeights; + fgWeights = 0x0; + } diff --git a/HBTAN/AliHBTWeights.h b/HBTAN/AliHBTWeights.h new file mode 100644 index 00000000000..a9c7dd5807e --- /dev/null +++ b/HBTAN/AliHBTWeights.h @@ -0,0 +1,23 @@ +#ifndef ALIHBTWEIGHTS_H +#define ALIHBTWEIGHTS_H + +#include + +class AliHBTPair; + +class AliHBTWeights: public TObject + { + public: + virtual ~AliHBTWeights(); + static Double_t Weight(const AliHBTPair* partpair){return (fgWeights)?fgWeights->GetWeight(partpair):0.0;} + virtual Double_t GetWeight(const AliHBTPair* partpair) = 0; + virtual void Set() = 0; + static AliHBTWeights* Instance() {return fgWeights;} + + protected: + static AliHBTWeights* fgWeights; + + ClassDef(AliHBTWeights,1) + }; + +#endif diff --git a/HBTAN/AliHBTLLWeightsPID.cxx b/HBTAN/AliHBTWeightsPID.cxx similarity index 95% rename from HBTAN/AliHBTLLWeightsPID.cxx rename to HBTAN/AliHBTWeightsPID.cxx index 0a5d4e63b48..9b439b42c1b 100644 --- a/HBTAN/AliHBTLLWeightsPID.cxx +++ b/HBTAN/AliHBTWeightsPID.cxx @@ -7,7 +7,7 @@ //Author: Ludmila Malinina, JINR (malinina@sunhe.jinr.ru) //----------------------------------------------------------- -#include "AliHBTLLWeightsPID.h" +#include "AliHBTWeightsPID.h" #include "AliHBTPair.h" #include "AliHBTParticle.h" #include @@ -15,11 +15,11 @@ #include #include -ClassImp(AliHBTLLWeightsPID) +ClassImp(AliHBTWeightsPID) -AliHBTLLWeightsPID* AliHBTLLWeightsPID::fgWeightsPID=NULL; +AliHBTWeightsPID* AliHBTWeightsPID::fgWeightsPID=NULL; -AliHBTLLWeightsPID::AliHBTLLWeightsPID() +AliHBTWeightsPID::AliHBTWeightsPID() { //ctor //initial parameters of model @@ -101,20 +101,20 @@ AliHBTLLWeightsPID::AliHBTLLWeightsPID() } -AliHBTLLWeightsPID* AliHBTLLWeightsPID::Instance() +AliHBTWeightsPID* AliHBTWeightsPID::Instance() { //Creates an instance of the class //or returns pointer to already existing one if (fgWeightsPID) { return fgWeightsPID; } else { - fgWeightsPID = new AliHBTLLWeightsPID(); + fgWeightsPID = new AliHBTWeightsPID(); return fgWeightsPID; } } -Double_t AliHBTLLWeightsPID::GetWeightPID(const AliHBTPair* trackpair) +Double_t AliHBTWeightsPID::GetWeightPID(const AliHBTPair* trackpair) { //Calculates the weight of "trackpair" AliHBTParticle *track1 = trackpair->Particle1(); diff --git a/HBTAN/AliHBTLLWeightsPID.h b/HBTAN/AliHBTWeightsPID.h similarity index 66% rename from HBTAN/AliHBTLLWeightsPID.h rename to HBTAN/AliHBTWeightsPID.h index ee6aa71f3c7..8dbd1bfee50 100644 --- a/HBTAN/AliHBTLLWeightsPID.h +++ b/HBTAN/AliHBTWeightsPID.h @@ -1,5 +1,5 @@ -#ifndef ALIHBTLLWEIGHTSPID_H -#define ALIHBTLLWEIGHTSPID_H +#ifndef ALIHBTWeightSPID_H +#define ALIHBTWeightSPID_H ///////////////////////////////////////////////////////////// // //This class introduces the weights calculated according @@ -14,21 +14,21 @@ class TF1; class TH1; class AliHBTPair; -class AliHBTLLWeightsPID: public TObject +class AliHBTWeightsPID: public TObject { public: - AliHBTLLWeightsPID(); - AliHBTLLWeightsPID(const AliHBTLLWeightsPID &source):TObject(source) { + AliHBTWeightsPID(); + AliHBTWeightsPID(const AliHBTWeightsPID &source):TObject(source) { //Copy ctor needed by the coding conventions but not used - Fatal("AliHBTLLWeightsPID","copy ctor not implemented"); + Fatal("AliHBTWeightsPID","copy ctor not implemented"); } - AliHBTLLWeightsPID & operator=(const AliHBTLLWeightsPID &/*source*/) { + AliHBTWeightsPID & operator=(const AliHBTWeightsPID &/*source*/) { //Assignment operator needed by the coding conventions but not used - Fatal("AliHBTLLWeightsPID","assignment operator not implemented"); + Fatal("AliHBTWeightsPID","assignment operator not implemented"); return * this; } - virtual ~AliHBTLLWeightsPID(){;} - static AliHBTLLWeightsPID* Instance(); + virtual ~AliHBTWeightsPID(){;} + static AliHBTWeightsPID* Instance(); Double_t GetWeightPID(const AliHBTPair* trackpair); //get weight calculated Batyunia's algorithm @@ -38,7 +38,7 @@ class AliHBTLLWeightsPID: public TObject Float_t fEfficTOF1; // ...? Float_t fEfficTOF2; // ...? - static AliHBTLLWeightsPID *fgWeightsPID;// pointer to wrapper of Fortran Lednicky code + static AliHBTWeightsPID *fgWeightsPID;// pointer to wrapper of Fortran Lednicky code TH1 *fPtK; //comment? TH1 *fPtKefftpc;//comment? TH1 *fPtKefftpcboth;//comment? @@ -52,7 +52,7 @@ class AliHBTLLWeightsPID: public TObject TF1 *fEffic3polTOF;//comment? TF1 *fEffic4polTOF;//comment? - ClassDef(AliHBTLLWeightsPID,1) + ClassDef(AliHBTWeightsPID,1) }; #endif diff --git a/HBTAN/HBTAnalysisLinkDef.h b/HBTAN/HBTAnalysisLinkDef.h index c0efe59b2d8..363e119637a 100644 --- a/HBTAN/HBTAnalysisLinkDef.h +++ b/HBTAN/HBTAnalysisLinkDef.h @@ -156,21 +156,26 @@ #pragma link C++ class AliHBTTwoTrackEffFctn+; #pragma link C++ class AliHBTTwoTrackEffFctn3D+; +#pragma link C++ class AliHBTWeights+; +#pragma link C++ class AliHBTCrab+; #pragma link C++ class AliHBTLLWeights+; -#pragma link C++ class AliHBTLLWeightQInvFctn+; -#pragma link C++ class AliHBTLLWeightQOutFctn+; -#pragma link C++ class AliHBTLLWeightQSideFctn+; -#pragma link C++ class AliHBTLLWeightQLongFctn+; -#pragma link C++ class AliHBTLLWeightTwoKStarFctn+; -#pragma link C++ class AliHBTLLWeightQOutQSideFctn+; -#pragma link C++ class AliHBTLLWeightQOutQLongFctn+; -#pragma link C++ class AliHBTLLWeightQSideQLongFctn+; - -#pragma link C++ class AliHBTLLWeightTheorQInvFctn+; -#pragma link C++ class AliHBTLLWeightTheorQOutFctn+; -#pragma link C++ class AliHBTLLWeightTheorQSideFctn+; -#pragma link C++ class AliHBTLLWeightTheorQLongFctn+; -#pragma link C++ class AliHBTLLWeightsPID+; + +#pragma link C++ class AliHBTWeightQInvFctn+; +#pragma link C++ class AliHBTWeightQOutFctn+; +#pragma link C++ class AliHBTWeightQSideFctn+; +#pragma link C++ class AliHBTWeightQLongFctn+; +#pragma link C++ class AliHBTWeightTwoKStarFctn+; +#pragma link C++ class AliHBTWeightQOutQSideFctn+; +#pragma link C++ class AliHBTWeightQOutQLongFctn+; +#pragma link C++ class AliHBTWeightQSideQLongFctn+; + +#pragma link C++ class AliHBTWeightTheorQInvFctn+; +#pragma link C++ class AliHBTWeightTheorQOutFctn+; +#pragma link C++ class AliHBTWeightTheorQSideFctn+; +#pragma link C++ class AliHBTWeightTheorQLongFctn+; +#pragma link C++ class AliHBTWeightTheorOSLFctn+; +#pragma link C++ class AliHBTWeightsPID+; + #pragma link C++ class AliHBTPositionRandomizer+; #pragma link C++ class AliHBTRndm+; #pragma link C++ class AliHBTRndmGaussBall+; @@ -178,6 +183,4 @@ #pragma link C++ class AliHBTMonVyDistributionVsVxFctn+; #pragma link C++ class AliHBTMonRtDistributionVsVzFctn+; //#pragma link C++ class AliHBTGoComPair+; -//#pragma link C++ class AliHBTCorrectQ3DCorrelFctn+; -//#pragma link C++ class AliHBTCorrectQInvCorrelFctn+; #endif diff --git a/HBTAN/libHBTAN.pkg b/HBTAN/libHBTAN.pkg index 6dda49eb980..da9bbe146f5 100644 --- a/HBTAN/libHBTAN.pkg +++ b/HBTAN/libHBTAN.pkg @@ -1,4 +1,4 @@ -SRCS = AliHBTAnalysis.cxx \ +SRCS = AliHBTWeights.cxx AliHBTCrab.cxx AliHBTAnalysis.cxx \ AliHBTReader.cxx AliHBTReaderKineTree.cxx \ AliHBTReaderESD.cxx AliHBTReaderInternal.cxx \ AliHBTReaderTPC.cxx AliHBTReaderPPprod.cxx \ @@ -10,8 +10,8 @@ AliHBTFunction.cxx AliHBTCorrelFctn.cxx \ AliHBTMonitorFunction.cxx AliHBTTwoTrackEffFctn.cxx \ AliHBTQResolutionFctns.cxx AliHBTQDistributionFctns.cxx \ AliHBTMonDistributionFctns.cxx AliHBTMonResolutionFctns.cxx \ -AliHBTLLWeights.cxx AliHBTLLWeightFctn.cxx \ -AliHBTLLWeightsPID.cxx AliHBTLLWeightTheorFctn.cxx \ +AliHBTLLWeights.cxx AliHBTWeightFctn.cxx \ +AliHBTWeightsPID.cxx AliHBTWeightTheorFctn.cxx \ AliHBTPositionRandomizer.cxx AliHBTEventBuffer.cxx \ AliHBTCorrFitFctn.cxx AliHBTCorrectCorrelFctn.cxx diff --git a/HBTAN/volya_complex.h b/HBTAN/volya_complex.h new file mode 100644 index 00000000000..76603a716c4 --- /dev/null +++ b/HBTAN/volya_complex.h @@ -0,0 +1,799 @@ +//This file is part of CRAB +//http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html + +/*-------------------------------------------------- +complex.cxx +Version 1.0.1 + +-- Some modification for speed +-- Some protections from overflow and underflow +-- Solving of equations of order 2,3,4: + Solve2, Solve3, Solve4; +--------------------------------------------------*/ + +#ifndef __COMPLEX_CXX__ +#define __COMPLEX_CXX__ + +#ifndef small_epsilon +#define small_epsilon 1e-12 +#endif + +#ifndef PI +#define PI 3.141592653589793238462643 +#endif + +#include +#include +#include + +/////////////////////////// Helper Error Function + +void ComplexError (char *msg) +{ + cout << endl << msg << endl; + exit (1); +} + +void ComplexWarning (char *msg) +{ + cout << endl << msg << endl; +} + +class Complex +{ +public: + double Re,Im; + + Complex (double a=0, double b=0): Re(a),Im(b) {} + + double GetRe () { return Re; } + double GetIm () { return Im; } + void SetRe (double a) { Re = a; } + void SetIm (double a) { Im = a; } + void Set (double a, double b) { Re = a; Im = b; } + void SetPolar (double, double); + void SetAbs (double); + void SetArg (double); + Complex operator ! () const; + Complex operator ! (); + Complex operator + () const; + Complex operator + (); + Complex operator - () const; + Complex operator - (); + Complex operator += (const Complex &); + Complex operator += (double); + Complex operator -= (const Complex &); + Complex operator -= (double); + Complex operator *= (const Complex &); + Complex operator *= (double); + Complex operator /= (const Complex &); + Complex operator /= (double); + + friend ostream& operator << (ostream&, const Complex &); + friend istream& operator >> (istream&, Complex &); +}; + + const Complex ImUnit (0,1); + + void Zero (Complex &); + Complex Polar (double, double); + double real (const Complex &); + double imag (const Complex &); + double Arg (const Complex &); + double phase(const Complex &); + double abs (const Complex &); + void Solve2 (Complex*, const Complex&, const Complex&); + void Solve3 (Complex*, const Complex&, const Complex&, const Complex&); + void Solve4 (Complex*, const Complex&, const Complex&, const Complex&, const Complex&); + Complex Solve2 (const Complex&, const Complex&, int RootNumber); + Complex Solve3 (const Complex&, const Complex&, const Complex&, int RootNumber); + Complex Solve4 (const Complex&, const Complex&, const Complex&, const Complex&, int RootNumber); + Complex conj (const Complex &); + Complex sqr (const Complex &); + Complex sqrt (const Complex &, int=0); + Complex pow (Complex, int); + Complex pow (Complex, double); + Complex pow (Complex, Complex); + Complex root (const Complex &, int, int=0); + Complex exp (const Complex &); + Complex log (const Complex &); + Complex sin (const Complex &); + Complex cos (const Complex &); + Complex tan (const Complex &); + Complex cot (const Complex &); + Complex sec (const Complex &); + Complex csc (const Complex &); + Complex sinh (const Complex &); + Complex cosh (const Complex &); + Complex tanh (const Complex &); + Complex coth (const Complex &); + Complex sech (const Complex &); + Complex csch (const Complex &); + Complex asin (const Complex &, int=0); + Complex acos (const Complex &, int=0); + Complex atan (const Complex &); + Complex acot (const Complex &); + Complex asec (const Complex &, int=0); + Complex acsc (const Complex &, int=0); + Complex asinh(const Complex &, int=0); + Complex acosh(const Complex &, int=0); + Complex atanh(const Complex &); + Complex acoth(const Complex &); + Complex asech(const Complex &, int=0); + Complex acsch(const Complex &, int=0); + + Complex operator ^ (const Complex &, int); + Complex operator ^ (const Complex &, double); + Complex operator ^ (const Complex &, const Complex &); + + double operator += (double&, const Complex &); + double operator -= (double&, const Complex &); + double operator *= (double&, const Complex &); + double operator /= (double&, const Complex &); + + Complex operator + (const Complex &, const Complex &); + Complex operator + (const Complex &, double); + Complex operator + (double, const Complex &); + Complex operator - (const Complex &, const Complex &); + Complex operator - (const Complex &, double); + Complex operator - (double, const Complex &); + Complex operator * (const Complex &, const Complex &); + Complex operator * (const Complex &, double); + Complex operator * (double, const Complex &); + Complex operator / (const Complex &, const Complex &); + Complex operator / (const Complex &, double); + Complex operator / (double, const Complex &); + + int operator == (const Complex &, double); + int operator == (double, const Complex &); + int operator == (const Complex &, const Complex &); + int operator != (const Complex &, double); + int operator != (double, const Complex &); + int operator != (const Complex &, const Complex &); + int operator >= (const Complex &, double); + int operator >= (double, const Complex &); + int operator >= (const Complex &, const Complex &); + int operator <= (const Complex &, double); + int operator <= (double, const Complex &); + int operator <= (const Complex &, const Complex &); + int operator > (const Complex &, double); + int operator > (double, const Complex &); + int operator > (const Complex &, const Complex &); + int operator < (const Complex &, double); + int operator < (double, const Complex &); + int operator < (const Complex &, const Complex &); + +//////////////////////////////////////////////// + +inline double sqr(double a) { return a*a; } + +inline Complex Complex::operator + () const +{ + return *this; +} + +inline Complex Complex::operator + () +{ + return *this; +} + +inline Complex Complex::operator - () const +{ + return Complex (-Re, -Im); +} + +inline Complex Complex::operator - () +{ + return Complex (-Re, -Im); +} + +inline Complex Complex::operator ! () const // complex conjugate // +{ + return Complex (Re, -Im); +} + +inline Complex Complex::operator ! () // complex conjugate // +{ + return Complex (Re, -Im); +} + +/////////////////////// +=,-=,*=,/= + +double operator += (double &a, const Complex &b) +{ + if (fabs(b.Im) < small_epsilon) a+=b.Re; + else + ComplexError ("Error in double+=Complex: Complex is not double number"); + return a; +} + +double operator -= (double &a, const Complex &b) +{ + if (fabs(b.Im) < small_epsilon) a-=b.Re; + else + ComplexError ("Error in double -=Complex: Complex is not double number"); + return a; +} + +double operator *= (double &a, const Complex &b) +{ + if (fabs(b.Im) < small_epsilon) a*=b.Re; + else + ComplexError ("Error in double *=Complex: Complex is not double number"); + return a; +} + +double operator /= (double &a, const Complex &b) +{ + if (fabs(b.Im) < small_epsilon) a/=b.Re; + else + ComplexError ("Error in double /=Complex: Complex is not double number"); + return a; +} + +inline Complex Complex::operator += (const Complex &a) +{ + Re+=a.Re; + Im+=a.Im; + return *this; +} + +inline Complex Complex::operator += (double a) +{ + Re+=a; + return *this; +} + +inline Complex Complex::operator -= (const Complex &a) +{ + Re-=a.Re; + Im-=a.Im; + return *this; +} + +inline Complex Complex::operator -= (double a) +{ + Re-=a; + return *this; +} + +Complex Complex::operator *= (const Complex &a) +{ + double t1=Re*a.Re, t2=Im*a.Im; + Im = (Re+Im) * (a.Re+a.Im) - t1 - t2; + Re = t1 - t2; + return *this; +} + +inline Complex Complex::operator *= (double a) +{ + Re*=a; + Im*=a; + return *this; +} + +Complex Complex::operator /= (const Complex &a) +{ + double t1, t2, temp; + if (fabs(a.Re) >= fabs(a.Im)) + { + t1 = a.Im / a.Re; + t2 = a.Re + a.Im * t1; + temp = (Re + Im * t1) / t2; + Im = (Im - Re * t1) / t2; + Re = temp; + } + else + { + t1 = a.Re / a.Im; + t2 = a.Re * t1 + a.Im; + temp = (Re * t1 + Im) / t2; + Im = (Im * t1 - Re) / t2; + Re = temp; + } + return *this; +} + +inline Complex Complex::operator /= (double a) +{ + Re/=a; + Im/=a; + return *this; +} + +/////////////////////////// +,-,*,/ + +inline Complex operator + (const Complex &a, const Complex &b) +{ + return Complex (a.Re+b.Re, a.Im+b.Im); +} + +inline Complex operator + (const Complex &a, double b) +{ + return Complex (a.Re+b, a.Im); +} + +inline Complex operator + (double b, const Complex &a) +{ + return Complex (a.Re+b, a.Im); +} + +inline Complex operator - (const Complex &a, const Complex &b) +{ + return Complex (a.Re-b.Re, a.Im-b.Im); +} + +inline Complex operator - (const Complex &a, double b) +{ + return Complex (a.Re-b, a.Im); +} + +inline Complex operator - (double b, const Complex &a) +{ + return Complex (b-a.Re, -a.Im); +} + +Complex operator * (const Complex &a, const Complex &b) +{ + double t1 = a.Re * b.Re; + double t2 = a.Im * b.Im; + return Complex (t1 - t2, (a.Re+a.Im) * (b.Re+b.Im) - t1 - t2); +} + +inline Complex operator * (const Complex &a, double b) +{ + return Complex (a.Re*b, a.Im*b); +} + +inline Complex operator * (double b, const Complex &a) +{ + return Complex (a.Re*b, a.Im*b); +} + +inline Complex operator / (const Complex &a, double b) +{ + return Complex (a.Re/b, a.Im/b); +} + +inline Complex operator / (const Complex &a, const Complex &b) +{ + double t1, t2; + if (fabs(b.Re) >= fabs(b.Im)) + { + t1 = b.Im / b.Re; + t2 = b.Re + b.Im * t1; + return Complex ((a.Re + a.Im * t1) / t2, (a.Im - a.Re * t1) / t2); + } + else + { + t1 = b.Re / b.Im; + t2 = b.Re * t1 + b.Im; + return Complex ((a.Re * t1 + a.Im) / t2, (a.Im * t1 - a.Re) / t2); + } +} + +inline Complex operator / (double b, const Complex &a) +{ + return (Complex (b,0)) / a; +} + +//////////////////////// <<,>> + +ostream& operator << (ostream &stream, const Complex &a) +{ + stream<<" "<> (istream &stream, Complex &a) +{ + stream>>a.Re>>a.Im; + return stream; +} + +/////////////////////// ==,!= + +inline int operator == (const Complex &a, double b) +{ + return (a.Re==b) && (fabs(a.Im)small_epsilon); +} + +inline int operator != (double b, const Complex &a) +{ + return (a.Re!=b) || (fabs(a.Im)>small_epsilon); +} + +inline int operator != (const Complex &a, const Complex &b) +{ + return (a.Re!=b.Re) || (a.Im!=b.Im); +} + +/////////////////////// >=,<= + +inline int operator >= (const Complex &a, double b) +{ + return abs(a) >= b; +} + +inline int operator >= (double b, const Complex &a) +{ + return b >= abs(a); +} + +inline int operator >= (const Complex &a, const Complex &b) +{ + return abs(a) >= abs(b); +} + +inline int operator <= (const Complex &a, double b) +{ + return abs(a) <= b; +} + +inline int operator <= (double b, const Complex &a) +{ + return b <= abs(a); +} + +inline int operator <= (const Complex &a, const Complex &b) +{ + return abs(a) <= abs(b); +} + +/////////////////////// >,< + +inline int operator > (const Complex &a, double b) +{ + return abs(a) > b; +} + +inline int operator > (double b, const Complex &a) +{ + return b > abs(a); +} + +inline int operator > (const Complex &a, const Complex &b) +{ + return abs(a) > abs(b); +} + +inline int operator < (const Complex &a, double b) +{ + return abs(a) < b; +} + +inline int operator < (double b, const Complex &a) +{ + return b < abs(a); +} + +inline int operator < (const Complex &a, const Complex &b) +{ + return abs(a) < abs(b); +} + +///////////////////////// Functions +double real (const Complex &a) +{ + return a.Re; +} +double imag (const Complex &a) +{ + return a.Im; +} +double abs (const Complex &a) +{ + if (a.Im == 0) return fabs(a.Re); + if (a.Re == 0) return fabs(a.Im); + double R = fabs(a.Re), I = fabs(a.Im); + return (R >= I) ? + (R * sqrt (1 + sqr(a.Im/a.Re))) : + (I * sqrt (1 + sqr(a.Re/a.Im))) ; +} + +double Arg (const Complex &a) +{ + return atan2(a.Im,a.Re); +} + +double phase (const Complex &a) +{ + return atan2(a.Im,a.Re); +} + +inline void Zero (Complex &a) { a.Re=a.Im=0; } + +Complex conj (const Complex &a) { return !a; } +Complex pow (Complex a, int n) { return a^n; } +Complex pow (Complex a, double n) { return a^n; } +Complex pow (Complex a, Complex b) { return a^b; } + +inline Complex sqr (const Complex &a) { return a*a; } + +Complex sqrt (const Complex &a, int flag) +{ + if ((a.Re>=0) && (fabs(a.Im)= I) ? + sqrt (R/2 * ( 1 + sqrt (1 + sqr(a.Im/a.Re)))): + sqrt (I/2 * (R/I + sqrt (1 + sqr(a.Re/a.Im)))); + Complex c; + if (a.Re >= 0) + { + c.Re = w; + c.Im = a.Im / (2*w); + } + else + { + c.Re = I / (2*w); + c.Im = (a.Im >= 0) ? w : -w; + } + return ((flag && (c.Re<0)) || (!flag && (c.Re>=0))) ? c : -c; +} + +Complex operator ^ (const Complex &a, int n) +{ + Complex c(1,0); + if (n==0) return c; + if (n>0) for (int i=0;in;j--) c*=a; + if (n>0) return c; + else return 1/c; +} + +Complex operator ^ (const Complex &a, double n) +{ + return exp(n*log(a)); +} + +Complex operator ^ (const Complex &a, const Complex &b) +{ + return exp(b*log(a)); +} + +Complex root (const Complex &z, int n, int k) +{ + double c=exp(log(abs(z))/n); + double t=(Arg(z)+2*PI*k)/n; + return Complex (c*cos(t), c*sin(t)); +} + +Complex exp (const Complex &a) +{ + double t=exp(a.Re); + return Complex (t*cos(a.Im), t*sin(a.Im)); +} + +Complex log (const Complex &a) +{ + if (a==0) + ComplexError("Error in function log(Complex): argument is 0"); + return Complex (log(abs(a)), Arg(a)); +} + +Complex sin (const Complex &a) +{ + return Complex (sin(a.Re)*cosh(a.Im), cos(a.Re)*sinh(a.Im)); +} + +Complex cos (const Complex &a) +{ + return Complex (cos(a.Re)*cosh(a.Im), -sin(a.Re)*sinh(a.Im)); +} + +Complex tan (const Complex &a) +{ + return sin(a)/cos(a); +} + +Complex cot (const Complex &a) +{ + return cos(a)/sin(a); +} + +Complex sec (const Complex &a) +{ + return 1/cos(a); +} + +Complex csc (const Complex &a) +{ + return 1/sin(a); +} + +Complex sinh (const Complex &a) +{ + return Complex (sinh(a.Re)*cos(a.Im), cosh(a.Re)*sin(a.Im)); +} + +Complex cosh (const Complex &a) +{ + return Complex (cosh(a.Re)*cos(a.Im), sinh(a.Re)*sin(a.Im)); +} + +Complex tanh (const Complex &a) +{ + return sinh(a)/cosh(a); +} + +Complex coth (const Complex &a) +{ + return cosh(a)/sinh(a); +} + +Complex sech (const Complex &a) +{ + return 1/cosh(a); +} + +Complex csch (const Complex &a) +{ + return 1/sinh(a); +} +//////////////////////// Inverce trigonometric functions + +Complex asin (const Complex &a, int flag) +{ + return -ImUnit * log(ImUnit*a + sqrt(1-sqr(a), flag)); +} + +Complex acos (const Complex &a, int flag) +{ + return -ImUnit * log(a + ImUnit*sqrt(1-sqr(a), flag)); +} + +Complex atan (const Complex &a) +{ + return ImUnit/2 * log((ImUnit+a)/(ImUnit-a)); +} + +Complex acot (const Complex &a) +{ + return ImUnit/2 * log((a-ImUnit)/(a+ImUnit)); +} + +Complex asec (const Complex &a, int flag) +{ + return acos(1/a, flag); +} + +Complex acsc (const Complex &a, int flag) +{ + return asin(1/a, flag); +} + +Complex asinh (const Complex &a, int flag) +{ + return log(a + sqrt(sqr(a)+1, flag)); +} + +Complex acosh (const Complex &a, int flag) +{ + return log(a + sqrt(sqr(a)-1, flag)); +} + +Complex atanh (const Complex &a) +{ + return log((1+a)/(1-a)) / 2; +} + +Complex acoth (const Complex &a) +{ + return log((a+1)/(a-1)) / 2; +} + +Complex asech (const Complex &a, int flag) +{ + return acosh(1/a, flag); +} + +Complex acsch (const Complex &a, int flag) +{ + return asinh(1/a, flag); +} + +Complex Polar (double a, double b) +{ + return Complex (a*cos(b), a*sin(b)); +} + +void Complex::SetPolar (double a, double b) +{ + Re=a*cos(b); + Im=a*sin(b); +} + +void Complex::SetAbs (double a) +{ + double b=Arg(*this); + Re=a*cos(b); + Im=a*sin(b); +} + +void Complex::SetArg (double b) +{ + double a=abs(*this); + Re=a*cos(b); + Im=a*sin(b); +} + +void Solve2 (Complex* z, const Complex &b, const Complex &c) + // finding z of equation: z^2 + b*z + c = 0 +{ + Complex t = sqrt(sqr(b)-4*c); + Complex q = ((!b * t).Re >= 0) ? (-(b + t) / 2) : (-(b - t) / 2); + z[0] = q; + z[1] = c/q; +} + +Complex Solve2 (const Complex &b, const Complex &c, int RootNumber) +{ + if ((RootNumber < 0) || (RootNumber > 1)) + ComplexError ("Error in Solve2: wrong root number"); + Complex z[2]; + Solve2 (z,b,c); + return z[RootNumber]; +} + +void Solve3 (Complex* z, const Complex &a2, const Complex &a1, const Complex &a0) + // finding z of equation: z^3 + a2*z^2 + a1*z + a0 = 0 +{ + Complex q, r, t, a, b, zero(0,0); + q = (sqr(a2) - 3*a1) / 9; + r = (2*(a2^3) - 9*a2*a1 + 27*a0) / 54; + t = sqrt((r^2) - (q^3)); + a = ((!r * t) >=0.0) ? -((r+t)^(1.0/3)) : -((r-t)^(1.0/3)); + b = ((a == zero) ? zero : (q/a)); + z[0] = -(a+b)/2 - a2/3 + ImUnit*sqrt(3)*(a-b)/2; + z[1] = -(a+b)/2 - a2/3 - ImUnit*sqrt(3)*(a-b)/2; + z[2] = a + b - a2/3; +} + +Complex Solve3 (const Complex &a2, const Complex &a1, const Complex &a0, int RootNumber) +{ + if ((RootNumber < 0) || (RootNumber > 2)) + ComplexError ("Error in Solve3: wrong root number"); + Complex z[3]; + Solve3 (z,a2,a1,a0); + return z[RootNumber]; +} + +void Solve4 (Complex* z, const Complex &a3, const Complex &a2, const Complex &a1, const Complex &a0) + // finding z of equation: z^4 + a3*z^3 + a2*z^2 + a1*z + a0 = 0 +{ + Complex u[3], t1, t2, t; + Solve3 (u, -a2, (a1*a3 - 4*a0), -((a1^2) + a0*(a3^2) - 4*a0*a2)); + t = *u; + t1 = sqrt((a3^2)/4 - a2 + t); + t2 = sqrt((t^2)/4 - a0); + Solve2 ( z, (a3/2 - t1), (t/2 + t2)); + Solve2 (&(z[2]), (a3/2 + t1), (t/2 - t2)); +} + +Complex Solve4 (const Complex &a3, const Complex &a2, const Complex &a1, const Complex &a0, int RootNumber) +{ + if ((RootNumber < 0) || (RootNumber > 3)) + ComplexError ("Error in Solve4: wrong root number"); + Complex z[4]; + Solve4 (z,a3,a2,a1,a0); + return z[RootNumber]; +} + + +#endif // __COMPLEX_CXX__ // -- 2.39.3