From 7f92929efeae1feda0e7d9fdf4f6ce4632d59444 Mon Sep 17 00:00:00 2001 From: skowron Date: Tue, 25 Jun 2002 17:50:25 +0000 Subject: [PATCH 1/1] C++ interface to Lednicky's algorithm. Initial revision --- HBTAN/AliHBTLLWeightFctn.cxx | 54 +++++ HBTAN/AliHBTLLWeightFctn.h | 39 ++++ HBTAN/AliHBTLLWeights.cxx | 384 +++++++++++++++++++++++++++++++++++ HBTAN/AliHBTLLWeights.h | 111 ++++++++++ 4 files changed, 588 insertions(+) create mode 100644 HBTAN/AliHBTLLWeightFctn.cxx create mode 100644 HBTAN/AliHBTLLWeightFctn.h create mode 100644 HBTAN/AliHBTLLWeights.cxx create mode 100644 HBTAN/AliHBTLLWeights.h diff --git a/HBTAN/AliHBTLLWeightFctn.cxx b/HBTAN/AliHBTLLWeightFctn.cxx new file mode 100644 index 00000000000..30e2d65a4fc --- /dev/null +++ b/HBTAN/AliHBTLLWeightFctn.cxx @@ -0,0 +1,54 @@ +#include +#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 index 00000000000..f23a746eb83 --- /dev/null +++ b/HBTAN/AliHBTLLWeightFctn.h @@ -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 index 00000000000..0593895a8df --- /dev/null +++ b/HBTAN/AliHBTLLWeights.cxx @@ -0,0 +1,384 @@ +#include +#include "AliHBTLLWeights.h" +#include "AliPDG.h" +#include "AliHBTPair.h" +#include "AliHBTParticle.h" +#include +#include +#include + +/*******************************************************************/ +/****** 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 index 00000000000..d9495022742 --- /dev/null +++ b/HBTAN/AliHBTLLWeights.h @@ -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 +#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 -- 2.43.0