]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
C++ interface to Lednicky's algorithm. Initial revision
authorskowron <skowron@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 25 Jun 2002 17:50:25 +0000 (17:50 +0000)
committerskowron <skowron@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 25 Jun 2002 17:50:25 +0000 (17:50 +0000)
HBTAN/AliHBTLLWeightFctn.cxx [new file with mode: 0644]
HBTAN/AliHBTLLWeightFctn.h [new file with mode: 0644]
HBTAN/AliHBTLLWeights.cxx [new file with mode: 0644]
HBTAN/AliHBTLLWeights.h [new file with mode: 0644]

diff --git a/HBTAN/AliHBTLLWeightFctn.cxx b/HBTAN/AliHBTLLWeightFctn.cxx
new file mode 100644 (file)
index 0000000..30e2d65
--- /dev/null
@@ -0,0 +1,54 @@
+#include <iostream.h>
+#include "AliHBTLLWeightFctn.h"
+#include "AliHBTLLWeights.h"
+
+//--for test--AliHBTLLWeightQInvFctn* yyy= new AliHBTLLWeightQInvFctn();
+
+ClassImp( AliHBTLLWeightQInvFctn )  
+/****************************************************************/
+AliHBTLLWeightQInvFctn::AliHBTLLWeightQInvFctn(Int_t nbins, Double_t maxXval, Double_t minXval):
+           AliHBTTwoPairFctn1D(nbins,maxXval,minXval)
+{
+  Rename("Correlation function, method of weights(Lednicky's algorithm)");
+}
+/****************************************************************/
+void  AliHBTLLWeightQInvFctn::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
+{
+
+  trackpair = CheckPair(trackpair);
+  partpair  = CheckPair(partpair);
+  if ( trackpair && partpair)     
+  {
+     fNumerator->Fill(trackpair->GetQInv(),
+          AliHBTLLWeights::Instance()->GetWeight(partpair));
+  }
+  
+} 
+                
+
+/****************************************************************/
+ void  AliHBTLLWeightQInvFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
+{
+  trackpair = CheckPair(trackpair);
+  partpair  = CheckPair(partpair);
+  if ( trackpair && partpair)  
+  {
+     fDenominator->Fill(trackpair->GetQInv());
+  }
+}
+/**************************************************************/
+TH1* AliHBTLLWeightQInvFctn::GetResult() 
+                                                                               
+{ 
+//returns ratio of numerator and denominator                                    
+ TH1* res = GetRatio(Scale());                                                  
+  
+ if(res)
+   {                                                                             
+    res->GetXaxis()->SetTitle("Qinv [GeV/c]");                                       
+    res->GetYaxis()->SetTitle("C(Qinv)");                                          
+    res->SetTitle(GetTitle());
+   }                                                                             
+  return res;                                                                     
+}                    
+                                                              
diff --git a/HBTAN/AliHBTLLWeightFctn.h b/HBTAN/AliHBTLLWeightFctn.h
new file mode 100644 (file)
index 0000000..f23a746
--- /dev/null
@@ -0,0 +1,39 @@
+//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 reconstructed tracks.
+//Weights are calculated with corresponding simulated particles momenta.
+//Denominator is filled with mixing unweighted reconstructed tracks.
+//One needs both pairs 
+//(simulated and recontructed), thus function is of class AliHBTTwoPairFctn1D.
+
+#ifndef ALIHBTLLWEIGHTFCTN_H
+#define ALIHBTLLWEIGHTFCTN_H
+#include "AliHBTFunction.h"
+
+
+class AliHBTLLWeights;
+
+class AliHBTLLWeightQInvFctn: public AliHBTTwoPairFctn1D
+{
+  friend class AliHBTOnePairFctn1D;
+
+  public:
+      AliHBTLLWeightQInvFctn(Int_t nbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0);
+      virtual  ~AliHBTLLWeightQInvFctn(){};
+      TH1* GetResult(); 
+
+      void   ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair);
+      void   ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair);
+      
+      Double_t GetValue(AliHBTPair* trackpair, AliHBTPair* partpair)         
+           { return trackpair->GetQInv()-partpair->GetQInv();} //isn't use                                                                    
+       
+
+  protected:
+
+  private:
+  public:
+     ClassDef(AliHBTLLWeightQInvFctn,1)
+};
+  
+#endif
diff --git a/HBTAN/AliHBTLLWeights.cxx b/HBTAN/AliHBTLLWeights.cxx
new file mode 100644 (file)
index 0000000..0593895
--- /dev/null
@@ -0,0 +1,384 @@
+#include <iostream.h>
+#include "AliHBTLLWeights.h"
+#include "AliPDG.h"
+#include "AliHBTPair.h"
+#include "AliHBTParticle.h"
+#include <TList.h>
+#include <TRandom.h>                                                                     
+#include <TMath.h>                                                                       
+
+/*******************************************************************/
+/******      ROUTINES    USED    FOR     COMMUNUCATION      ********/
+/********************     WITH      FORTRAN     ********************/
+/*******************************************************************/
+#ifndef WIN32
+# define led_bldata led_bldata_
+# define fsiini fsiini_
+# define ltran12 ltran12_
+# define fsiw fsiw_
+# define type_of_call
+#else
+# define led_bldata LED_BLDATA
+# define fsiini FSIINI
+# define ltran12 LTRAN12
+# define fsiw FSIW
+# define type_of_call _stdcall
+#endif
+/****************************************************************/
+extern "C" void type_of_call led_bldata(); 
+extern "C" void type_of_call fsiini();
+extern "C" void type_of_call ltran12();
+extern "C" void type_of_call fsiw();
+/**************************************************************/
+
+ClassImp(AliHBTLLWeights)  
+
+//Interface to Richard Lednicky's weight alghorithm (Fortran implementation)
+//Ludmila Malinina JINR
+//Piotr Krzysztof Skowronski CERN/WUT
+AliHBTLLWeights* AliHBTLLWeights::fgLLWeights=NULL; 
+
+AliHBTLLWeights::AliHBTLLWeights()
+{
+// Default Constructor 
+    fPID1 = 0;
+    fPID2 = 0;
+    SetRandomPosition();
+    SetColWithResidNuclSwitch();
+    SetStrongInterSwitch();
+    SetQuantumStatistics();
+    SetColoumb();
+    SetTest();
+}
+
+
+ AliHBTLLWeights* AliHBTLLWeights::Instance()
+{                                                                                             
+  if (fgLLWeights) 
+   { 
+     return fgLLWeights;                                                                   
+   } 
+  else 
+   {                                                                                  
+    fgLLWeights = new AliHBTLLWeights();                                                        
+    return fgLLWeights;                                                                   
+  }                                                                                         
+}                                                                                             
+                                            
+
+Double_t AliHBTLLWeights::GetWeight(const AliHBTPair* partpair)
+{
+    AliHBTParticle *part1 = partpair->Particle1();
+    AliHBTParticle *part2 = partpair->Particle2();
+   
+    if ( (part1 == 0x0) || (part2 == 0x0))
+     {
+      Error("GetWeight","Null particle pointer");
+      return 0.0;
+     }
+   
+    Double_t part1Momentum[]={part1->Px(),part1->Py(),part1->Pz()};
+    Double_t part2Momentum[]={part2->Px(),part2->Py(),part2->Pz()};
+              
+   
+    if ( (part1->Px() == part2->Px()) && (part1->Py() == part2->Py()) 
+       && (part1->Pz() == part2->Pz()) )
+      {//if particles have the same momentum 
+        return 0.0;
+      }
+
+             
+    if ((!fRandomPosition) && (part1->Vx()  == part2->Vx()) && (part1->Vy()  == part2->Vy())
+          && (part1->Vz()  == part2->Vz()) )
+      { //if particles have the same position       
+       return 0.0;
+      }
+
+//put momenta of particles in LAB into Fortran commons
+
+    FSI_MOM.P1X=part1Momentum[0];
+    FSI_MOM.P1Y=part1Momentum[1];
+    FSI_MOM.P1Z=part1Momentum[2];
+      
+    FSI_MOM.P2X=part2Momentum[0];
+    FSI_MOM.P2Y=part2Momentum[1];
+    FSI_MOM.P2Z=part2Momentum[2];
+        
+    if (fRandomPosition) 
+     {
+       Double_t rxcm = fsigma*gRandom->Gaus();
+       Double_t rycm = fsigma*gRandom->Gaus();
+       Double_t rzcm = fsigma*gRandom->Gaus();   
+    
+       FSI_PRF.X=rxcm*fwcons;
+       FSI_PRF.Y=rycm*fwcons;
+       FSI_PRF.Z=rzcm*fwcons;
+       FSI_PRF.T=0.;
+
+       Double_t rps=rxcm*rxcm+rycm*rycm+rzcm*rzcm;
+       Double_t rp=TMath::Sqrt(rps);
+       FSI_PRF.RP=rp;
+       FSI_PRF.RPS=rps;
+     }    
+          
+    ltran12();
+    fsiw();
+    return LEDWEIGHT.WEIN;
+ }
+/************************************************************/
+void AliHBTLLWeights::Init()
+ {
+//---------------------------------------------------------------------
+
+//initial parameters of model
+
+  FSI_NS.NS = approximationModel;  
+  
+  if(!ftest) LEDWEIGHT.ITEST=0;
+    
+  if(ftest)
+   {
+     LEDWEIGHT.ITEST=1;
+     if(fColoumbSwitch) FSI_NS.ICH =1;
+     else FSI_NS.ICH=0;
+     
+     if(fStrongInterSwitch) FSI_NS.ISI=1;
+     else FSI_NS.ISI=0;
+   
+     if(fQuantStatSwitch) FSI_NS.IQS=1;
+     else FSI_NS.IQS=0;
+   
+     if(fColWithResidNuclSwitch) FSI_NS.I3C=1;
+     else FSI_NS.I3C=0;
+   }
+  
+  if(fRandomPosition) LEDWEIGHT.IRANPOS=1;
+  else LEDWEIGHT.IRANPOS=0;
+
+  if ( (fPID1 == 0) || (fPID2 == 0) )
+   {
+    Fatal("Init","Particles types are not set");
+    return;//pro forma
+   }
+
+  FSI_NS.LL = GetPairCode(fPID1,fPID2);
+   
+  if (FSI_NS.LL == 0) 
+   {
+     Fatal("Init","Particles types are not supported");
+     return;//pro forma
+   }
+
+  TParticlePDG* tpart1 = TDatabasePDG::Instance()->GetParticle(fPID1);
+  if (tpart1 == 0x0)
+   {
+     Fatal("init","We can not find particle with ID=%d in our DataBase",fPID1);
+     return;
+   }
+  
+  FSI_POC.AM1=tpart1->Mass();
+  FSI_POC.C1=tpart1->Charge(); 
+
+  TParticlePDG* tpart2 = TDatabasePDG::Instance()->GetParticle(fPID2);
+  if (tpart2 == 0x0)
+   {
+     Fatal("init","We can not find particle with ID=%d in our DataBase",fPID2);
+     return;
+   }
+
+  FSI_POC.AM2=tpart2->Mass();
+  FSI_POC.C1=tpart2->Charge();
+
+  led_bldata();
+  fsiini();
+
+
+//constants for radii simulation 
+
+  if(fRandomPosition)
+   {
+    fsigma =TMath::Sqrt(2.)*fRadius;     
+    fwcons =FSI_CONS.W;
+   } 
+} 
+
+Int_t AliHBTLLWeights::GetPairCode(const AliHBTPair* partpair)
+{
+ return GetPairCode(partpair->Particle1()->GetPdgCode(),partpair->Particle2()->GetPdgCode());
+}
+
+Int_t AliHBTLLWeights::GetPairCode(Int_t pid1,Int_t pid2)
+{
+//   pairCode   1  2  3  4   5    6   7  8  9 10  11  12  13  14 15 16 17 18  19  20   21   22  23 24 25 26    27     28
+//   hpid:      n  p  n alfa pi+ pi0 pi+ n  p pi+ pi+ pi+ pi- K+ K+ K+ K-  d  d    t   t    K0  K0  d p  p      p      n
+//   lpid:      n  p  p alfa pi- pi0 pi+ d  d  K-  K+  p   p  K- K+ p  p   d alfa  t  alfa  K0  K0b t t alfa lambda lambda
+//   NS=1 y/n:  +  +  +  +   +    -   -  -  -  -   -   -   -  -  -  -  -   -  -    -    -    -  -   - -  -      -      -
+
+//alphas, deuterons and tyts are NOT supported here
+
+  Int_t chargefactor = 1;
+  Int_t hpid; //pid in higher row
+  Int_t lpid; //pid in lower row
+  Int_t code; //pairCode
+  
+  Bool_t swap;
+  
+//determine the order of selcetion in switch  
+  if (TMath::Abs(pid1) < TMath::Abs(pid2) ) 
+   {
+    if (pid1<0) chargefactor=-1;
+    hpid=pid2*chargefactor;
+    lpid=pid1*chargefactor;
+    swap = kFALSE;
+   } 
+  else 
+   {
+    if (pid2<0) chargefactor=-1;
+    hpid=pid1*chargefactor;
+    lpid=pid2*chargefactor;
+    swap = kTRUE;
+   }
+
+//Determine the pair code
+  switch (hpid) //switch on first  particle id
+   {
+     case kNeutron:
+      switch (lpid)
+       {
+         case kNeutron: 
+           code = 1;  //neutron neutron
+           break;
+        
+         case kProton: 
+           code = 3;  //neutron proton
+           break;
+           
+         case kLambda0: 
+           code = 28;  //neutron lambda
+           break;
+           
+         default: 
+           return 0; //given pair not supported
+           break;
+       }
+      break;
+
+     case kProton:
+      switch (lpid)
+       {
+         case kProton:
+           code = 2; //proton proton
+           break;
+           
+         case kLambda0: 
+           code = 27;//proton lambda
+           break;
+           
+         default: 
+           return 0; //given pair not supported
+           break;
+           
+       }
+      break;
+
+     case kPiPlus:
+      switch (lpid)
+       {
+         case kPiPlus:
+           code = 7; //piplus piplus
+           break;
+
+         case kPiMinus:
+           code = 5; //piplus piminus
+           break;
+        
+         case kKMinus:
+           code = 10; //piplus Kminus
+           break;
+
+         case kKPlus:
+           code = 11; //piplus Kplus
+           break;
+
+         case kProton:
+           code = 12; //piplus proton
+           chargefactor*=-1;
+           break;
+
+         default: 
+           return 0; //given pair not supported
+           break;
+       }
+      break;
+     case kPi0:
+      switch (lpid)
+       {
+         case kPi0:
+           code = 6;
+           break;
+           
+         default: 
+           return 0; //given pair not supported
+           break;
+       }
+      break;
+      
+     case kKPlus:
+      switch (lpid)
+       {
+         case kKMinus:
+           code = 14; //Kplus Kminus
+           break;
+
+         case kKPlus:
+           code = 15; //Kplus Kplus
+           break;
+
+         case kProton:
+           code = 16; //Kplus proton
+           break;
+           
+         default: 
+           return 0; //given pair not supported
+           break;
+       }
+      break;
+      
+     case kKMinus:
+      switch (lpid)
+       {
+         case kProton:
+           code = 17; //Kminus proton
+           chargefactor*=1;
+           break;
+           
+         default: 
+           return 0; //given pair not supported
+           break;
+       }
+      break;
+      
+     case kK0:
+      switch (lpid)
+       {
+         case kK0:
+           code = 2; //Kzero Kzero
+           break;
+         
+         case kK0Bar:
+           code = 17; //Kzero KzeroBar
+           break;
+
+         default: 
+           return 0; //given pair not supported
+           break;
+       }
+      break;
+
+     default: return 0;
+   }
+  return code;
+}
+
diff --git a/HBTAN/AliHBTLLWeights.h b/HBTAN/AliHBTLLWeights.h
new file mode 100644 (file)
index 0000000..d949502
--- /dev/null
@@ -0,0 +1,111 @@
+//This class introduce the weights calculation according with Lednicky's algorithm.
+//The detailed description of the algorithm can be found in comments to fortran code:
+//fsiw.f, fsiini.f  
+#ifndef ALIHBTLLWEIGHTS_H
+#define ALIHBTLLWEIGHTS_H
+
+#include <TObject.h>
+#include "WLedCOMMONS.h"
+
+class AliHBTPair;
+class AliHBTLLWeights: public TObject
+ {
+   public:
+     virtual ~AliHBTLLWeights(){;}
+     static AliHBTLLWeights* Instance();
+     
+     void Init(); //put the initial values in fortran commons fsiini, led_bldata
+     Double_t GetWeight(const AliHBTPair* partpair); //get weight calculated by Lednicky's algorithm
+
+     void SetTest(Bool_t rtest = kTRUE){ftest = rtest;} //if ftest=0: 
+     //physical values of the following  parameters are put automatically                       
+     //            in FSIINI (their values are not required)          
+     // ftest=1: any values of the following parameters are allowed,    
+     //the following parameters are required:                           
+
+     void SetColoumb(Bool_t col = kTRUE){fColoumbSwitch = col;}//: (ICH in fortran code) Coulomb interaction between the two particles ON (OFF)
+     void SetQuantumStatistics(Bool_t qss = kTRUE){fQuantStatSwitch = qss;}//IQS: quantum statistics for the two particles ON (OFF) //if non-identical particles automatically off
+     void SetStrongInterSwitch(Bool_t sis = kTRUE){fStrongInterSwitch = sis;}//ISI: strong interaction between the two particles ON (OFF)
+     void SetColWithResidNuclSwitch(Bool_t crn = kTRUE){fColWithResidNuclSwitch = crn;}//I3C: Coulomb interaction with residual nucleus ON (OFF)  
+     void SetApproxModel(Int_t ap){approximationModel=ap;}//NS in Fortran code, 
+     //   NS=1  Square well potential,                                                             
+     //   NS=3  not used                                                                           
+     //   NS=4  scattered wave approximated by the spherical wave,                                 
+     //   NS=2  same as NS=4 but the approx. of equal emission times in PRF                        
+     //         not required (t=0 approx. used in all other cases).      
+
+     
+     void SetRandomPosition(Bool_t rp = kTRUE){fRandomPosition = rp;} //ON=kTRUE(OFF=kFALSE)
+     // ON -- calculation of the Gauss source radii if the generator don't allows the source generation (for example MeVSim)
+     //if ON the following parameters are requested:
+     void SetR1dw(Double_t R){fRadius=R;}   //spherical source model radii                                                                                                                                                                                             
+     void SetLambdaw(Double_t la){flambda=la;}  //lambda=haoticity                                                   
+
+     
+     void SetParticlesTypes(Int_t pid1, Int_t pid2){fPID1 = pid1; fPID2 = pid2;} //set AliRoot particles types   
+    
+     void SetNucleusCharge(Double_t ch){fNuclCharge=ch;} // not used now  (see comments in fortran code)
+     void SetNucleusMass(Double_t mass){fNuclMass=mass;} // (see comments in fortran code)
+
+
+   protected:
+     
+     Bool_t ftest; 
+     Bool_t fColoumbSwitch; 
+     Bool_t fQuantStatSwitch; 
+     Bool_t fStrongInterSwitch;//Switches strong interactions TRUE=ON
+     Bool_t fColWithResidNuclSwitch;//Switches couloumb interaction 
+                                    //with residual nucleus TRUE=ON          
+     Double_t fNuclMass; //mass 
+     Double_t fNuclCharge; //charge    
+     
+     Bool_t  fRandomPosition;
+     Double_t fRadius;
+     Double_t flambda;
+
+     
+     Double_t wein;
+
+     Int_t approximationModel; //approximation used to calculate Bethe-Salpeter amplitude
+                             //   ==1  Square well potential,
+                             //   ==3  not used
+                             //   ==4  scattered wave approximated by the spherical wave,
+                             //   ==2  same as NS=4 but the approx. of equal emission times in PRF
+                             //         not required (t=0 approx. used in all other cases).
+                             //  Note: if ==2,4, the B-S amplitude diverges at zero distance r* in
+                             //         the two-particle c.m.s.; user can specify a cutoff AA in
+                             //         SUBROUTINE FSIINI, for example:
+                             //         IF(NS.EQ.2.OR.NS.EQ.4)AA=5.D0 !! in 1/GeV --> AA=1. fm
+     
+     Int_t fPID1;
+     Int_t fPID2;
+
+     static  AliHBTLLWeights *fgLLWeights;// pointer to wrapper of Fortran Lednicky code
+
+
+     static Int_t GetPairCode(Int_t pid1,Int_t pid2);
+     static Int_t GetPairCode(const AliHBTPair* partpair);//calculate automatically internal FSIW 
+         //     C----------------------------------------------------------------------               
+         //     C-   LL       1  2  3  4   5    6   7  8  9 10  11  12  13  14 15 16 17               
+         //     C-   part. 1: n  p  n alfa pi+ pi0 pi+ n  p pi+ pi+ pi+ pi- K+ K+ K+ K-               
+         //     C-   part. 2: n  p  p alfa pi- pi0 pi+ d  d  K-  K+  p   p  K- K+ p  p                
+         //     C   NS=1 y/n: +  +  +  +   +    -   -  -  -  -   -   -   -  -  -  -  -                
+         //     C----------------------------------------------------------------------               
+         //     C-   LL       18  19 20  21  22 23 24 25 26    27     28                              
+         //     C-   part. 1:  d  d   t  t   K0 K0  d p  p      p      n                              
+         //     C-   part. 2:  d alfa t alfa K0 K0b t t alfa lambda lambda                            
+         //     C   NS=1 y/n:  -  -   -  -   -  -   - -  -      +      +                              
+         //     C----------------------------------------------------------------------               
+                                                                                  
+
+     Double_t fsigma; //constants for spherical source model 
+     Double_t fwcons; //
+     
+   private:
+   AliHBTLLWeights();
+
+   public:
+     ClassDef(AliHBTLLWeights,1)
+ };
+
+#endif