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))
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)
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)
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)
--- /dev/null
+#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 <TMath.h>
+#include <TPDGCode.h>
+
+#include "volya_complex.h"
+
+//#include <complex>
+//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;
+}
+//===================================================================
+
--- /dev/null
+/* $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
}
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;
{
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;
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;
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();
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; k<nbinsZ; k++)
for (UInt_t j = offsetY; j<nbinsY; j++)
{
if ( fNumerator->GetBinContent(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;
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
+++ /dev/null
-#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
/**************************************************************/
AliHBTLLWeights::AliHBTLLWeights(const AliHBTLLWeights &/*source*/):
- TObject(),
+ AliHBTWeights(),
fTest(kTRUE),
fColoumbSwitch(kTRUE),
fQuantStatSwitch(kTRUE),
}
/************************************************************/
+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
#ifndef ALIHBTLLWEIGHTS_H
#define ALIHBTLLWEIGHTS_H
-#include <TObject.h>
+#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
#include "AliHBTPair.h"
#include "AliHBTParticle.h"
-#include "AliHBTLLWeights.h"
+#include "AliHBTWeights.h"
ClassImp(AliHBTPair)
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),
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),
}
/************************************************************************/
-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;
}
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
Bool_t fQInvLNotCalc;
void CalculateQInvL();
- Double_t fLLWeight;
- Bool_t ffLLWeightNotCalc;
+ Double_t fWeight;
+ Bool_t fWeightNotCalc;
Double_t fPxSum;
Double_t fPySum;
fKStarNotCalc = kTRUE;
fQInvLNotCalc = kTRUE;
fGammaCMSLCNotCalc = kTRUE;
- ffLLWeightNotCalc = kTRUE;
+ fWeightNotCalc = kTRUE;
}
/****************************************************************/
inline
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;
for (Int_t s = 0; s<kNSpecies; s++)
{
- if (w[s] == 0.0) continue;
+ Float_t pp = w[s];
+ if (pp == 0.0) continue;
Int_t pdgcode = charge*GetSpeciesPdgCode((ESpecies)s);
if(Pass(pdgcode)) continue; //check if we are intersted with particles of this type
}//for (Int_t i = 0;i<ntr; i++) -- loop over tracks
Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
- fParticlesEvent->GetNumberOfParticles(), fTracksEvent->GetNumberOfParticles(),
+ fTracksEvent->GetNumberOfParticles(), fParticlesEvent->GetNumberOfParticles(),
fNEventsRead,fCurrentEvent,fCurrentDir);
fCurrentEvent++;
-#include "AliHBTLLWeightFctn.h"
+#include "AliHBTWeightFctn.h"
/* $Id$ */
//_________________________________________________________________________
//
//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
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);
}
}
/**************************************************************/
-TH1* AliHBTLLWeightQInvFctn::GetResult()
+TH1* AliHBTWeightQInvFctn::GetResult()
{
//returns ratio of numerator and denominator
return GetRatio(Scale());
/**************************************************************************************/
/**************************************************************************************/
-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
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);
}
}
/**************************************************************/
-TH1* AliHBTLLWeightQOutFctn::GetResult()
+TH1* AliHBTWeightQOutFctn::GetResult()
{
//returns ratio of numerator and denominator
/*************************************************************************************/
/*************************************************************************************/
-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
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);
}
}
/**************************************************************/
-TH1* AliHBTLLWeightQLongFctn::GetResult()
+TH1* AliHBTWeightQLongFctn::GetResult()
{
//returns ratio of numerator and denominator
/*************************************************************************************/
/*************************************************************************************/
-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
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);
}
}
/**************************************************************/
-TH1* AliHBTLLWeightQSideFctn::GetResult()
+TH1* AliHBTWeightQSideFctn::GetResult()
{
//returns ratio of numerator and denominator
/*************************************************************************************/
/*************************************************************************************/
-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
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);
}
}
/**************************************************************/
-TH1* AliHBTLLWeightTwoKStarFctn::GetResult()
+TH1* AliHBTWeightTwoKStarFctn::GetResult()
{
//returns ratio of numerator and denominator
/*************************************************************************************/
/*************************************************************************************/
-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)
{
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);
}
}
/**************************************************************/
-TH1* AliHBTLLWeightQOutQSideFctn::GetResult()
+TH1* AliHBTWeightQOutQSideFctn::GetResult()
{
//returns result
return GetRatio(Scale());
/*************************************************************************************/
/*************************************************************************************/
-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)
{
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);
}
/**************************************************************/
-TH1* AliHBTLLWeightQOutQLongFctn::GetResult()
+TH1* AliHBTWeightQOutQLongFctn::GetResult()
{
//returns result
return GetRatio(Scale());
/*************************************************************************************/
/*************************************************************************************/
-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)
{
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);
}
}
/**************************************************************/
-TH1* AliHBTLLWeightQSideQLongFctn::GetResult()
+TH1* AliHBTWeightQSideQLongFctn::GetResult()
{
//returns result
return GetRatio(Scale());
-#ifndef ALIHBTLLWEIGHTQINVFCTN_H
-#define ALIHBTLLWEIGHTQINVFCTN_H
+#ifndef ALIHBTWeightQINVFCTN_H
+#define ALIHBTWeightQINVFCTN_H
/* $Id$ */
#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);
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);
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);
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);
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);
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)
};
// 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
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());
/*************************************************************/
/*************************************************************/
-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
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());
/*************************************************************/
/*************************************************************/
-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
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());
/*************************************************************/
/*************************************************************/
-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
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());
--- /dev/null
+#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
--- /dev/null
+#include "AliHBTWeights.h"
+
+ClassImp(AliHBTWeights)
+
+AliHBTWeights* AliHBTWeights::fgWeights = 0x0;
+
+AliHBTWeights::~AliHBTWeights()
+ {
+ delete fgWeights;
+ fgWeights = 0x0;
+ }
--- /dev/null
+#ifndef ALIHBTWEIGHTS_H
+#define ALIHBTWEIGHTS_H
+
+#include <TObject.h>
+
+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
//Author: Ludmila Malinina, JINR (malinina@sunhe.jinr.ru)
//-----------------------------------------------------------
-#include "AliHBTLLWeightsPID.h"
+#include "AliHBTWeightsPID.h"
#include "AliHBTPair.h"
#include "AliHBTParticle.h"
#include <TRandom.h>
#include <TH1.h>
#include <TF1.h>
-ClassImp(AliHBTLLWeightsPID)
+ClassImp(AliHBTWeightsPID)
-AliHBTLLWeightsPID* AliHBTLLWeightsPID::fgWeightsPID=NULL;
+AliHBTWeightsPID* AliHBTWeightsPID::fgWeightsPID=NULL;
-AliHBTLLWeightsPID::AliHBTLLWeightsPID()
+AliHBTWeightsPID::AliHBTWeightsPID()
{
//ctor
//initial parameters of model
}
-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();
-#ifndef ALIHBTLLWEIGHTSPID_H
-#define ALIHBTLLWEIGHTSPID_H
+#ifndef ALIHBTWeightSPID_H
+#define ALIHBTWeightSPID_H
/////////////////////////////////////////////////////////////
//
//This class introduces the weights calculated according
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
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?
TF1 *fEffic3polTOF;//comment?
TF1 *fEffic4polTOF;//comment?
- ClassDef(AliHBTLLWeightsPID,1)
+ ClassDef(AliHBTWeightsPID,1)
};
#endif
#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+;
#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
-SRCS = AliHBTAnalysis.cxx \
+SRCS = AliHBTWeights.cxx AliHBTCrab.cxx AliHBTAnalysis.cxx \
AliHBTReader.cxx AliHBTReaderKineTree.cxx \
AliHBTReaderESD.cxx AliHBTReaderInternal.cxx \
AliHBTReaderTPC.cxx AliHBTReaderPPprod.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
--- /dev/null
+//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 <iostream.h>
+#include <math.h>
+#include <stdlib.h>
+
+/////////////////////////// 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<<" "<<a.Re<<" "<<a.Im<<" ";
+ return stream;
+}
+
+istream& operator >> (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 (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)<small_epsilon))
+ return flag ? -sqrt(a.Re) : sqrt(a.Re);
+ double R = fabs(a.Re), I = fabs(a.Im);
+ double w = (R >= 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;i<n;i++) c*=a;
+ else for (int j=0;j>n;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__ //
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff