CRAB added
authorskowron <skowron@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 28 Oct 2003 17:59:19 +0000 (17:59 +0000)
committerskowron <skowron@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 28 Oct 2003 17:59:19 +0000 (17:59 +0000)
24 files changed:
HBTAN/AliHBTCorrFitFctn.cxx
HBTAN/AliHBTCorrelFctn.h
HBTAN/AliHBTCrab.cxx [new file with mode: 0644]
HBTAN/AliHBTCrab.h [new file with mode: 0644]
HBTAN/AliHBTFunction.cxx
HBTAN/AliHBTFunction.h
HBTAN/AliHBTLLWeightTheorFctn.h [deleted file]
HBTAN/AliHBTLLWeights.cxx
HBTAN/AliHBTLLWeights.h
HBTAN/AliHBTPair.cxx
HBTAN/AliHBTPair.h
HBTAN/AliHBTParticle.h
HBTAN/AliHBTReaderESD.cxx
HBTAN/AliHBTWeightFctn.cxx [moved from HBTAN/AliHBTLLWeightFctn.cxx with 64% similarity]
HBTAN/AliHBTWeightFctn.h [moved from HBTAN/AliHBTLLWeightFctn.h with 67% similarity]
HBTAN/AliHBTWeightTheorFctn.cxx [moved from HBTAN/AliHBTLLWeightTheorFctn.cxx with 57% similarity]
HBTAN/AliHBTWeightTheorFctn.h [new file with mode: 0644]
HBTAN/AliHBTWeights.cxx [new file with mode: 0644]
HBTAN/AliHBTWeights.h [new file with mode: 0644]
HBTAN/AliHBTWeightsPID.cxx [moved from HBTAN/AliHBTLLWeightsPID.cxx with 95% similarity]
HBTAN/AliHBTWeightsPID.h [moved from HBTAN/AliHBTLLWeightsPID.h with 66% similarity]
HBTAN/HBTAnalysisLinkDef.h
HBTAN/libHBTAN.pkg
HBTAN/volya_complex.h [new file with mode: 0644]

index 8c44d77..3661402 100644 (file)
@@ -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))
index bf9ea38..88dce65 100644 (file)
@@ -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 (file)
index 0000000..db1070c
--- /dev/null
@@ -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 <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;
+}
+//===================================================================
+
diff --git a/HBTAN/AliHBTCrab.h b/HBTAN/AliHBTCrab.h
new file mode 100644 (file)
index 0000000..35ccbca
--- /dev/null
@@ -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
index 767834a..aceda9f 100644 (file)
@@ -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; k<nbinsZ; k++)
     for (UInt_t j = offsetY; j<nbinsY; j++)
@@ -763,18 +760,18 @@ Double_t AliHBTFunction3D::Scale()
        {
         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;
index 1132c40..c5f91b3 100644 (file)
@@ -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 (file)
index 7acd78e..0000000
+++ /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
index 49e1868..32c19c7 100644 (file)
@@ -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
index 6a6d6f6..b625d54 100644 (file)
@@ -9,14 +9,15 @@
 #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
 
index a385fb6..fe5c020 100644 (file)
@@ -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; 
 }
index f01c592..7c79d67 100644 (file)
@@ -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 
index 7f3f9d9..3fa2047 100644 (file)
@@ -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;
 
index 9ea9657..232b67e 100644 (file)
@@ -198,7 +198,8 @@ Int_t AliHBTReaderESD::ReadNext()
 
         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 
@@ -239,7 +240,7 @@ Int_t AliHBTReaderESD::ReadNext()
       }//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++;
similarity index 64%
rename from HBTAN/AliHBTLLWeightFctn.cxx
rename to HBTAN/AliHBTWeightFctn.cxx
index 0991d5d..2d2b91f 100644 (file)
@@ -1,4 +1,4 @@
-#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
@@ -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());
similarity index 67%
rename from HBTAN/AliHBTLLWeightFctn.h
rename to HBTAN/AliHBTWeightFctn.h
index 97c00d0..278eb5c 100644 (file)
@@ -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)
  
 };
 
similarity index 57%
rename from HBTAN/AliHBTLLWeightTheorFctn.cxx
rename to HBTAN/AliHBTWeightTheorFctn.cxx
index 515ab87..5bd7a34 100644 (file)
@@ -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 (file)
index 0000000..6ec3184
--- /dev/null
@@ -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 (file)
index 0000000..9952c11
--- /dev/null
@@ -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 (file)
index 0000000..a9c7dd5
--- /dev/null
@@ -0,0 +1,23 @@
+#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
similarity index 95%
rename from HBTAN/AliHBTLLWeightsPID.cxx
rename to HBTAN/AliHBTWeightsPID.cxx
index 0a5d4e6..9b439b4 100644 (file)
@@ -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 <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
@@ -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();
similarity index 66%
rename from HBTAN/AliHBTLLWeightsPID.h
rename to HBTAN/AliHBTWeightsPID.h
index ee6aa71..8dbd1bf 100644 (file)
@@ -1,5 +1,5 @@
-#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
    
@@ -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  
index c0efe59..363e119 100644 (file)
 #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
index 6dda49e..da9bbe1 100644 (file)
@@ -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 (file)
index 0000000..76603a7
--- /dev/null
@@ -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 <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__ //