/*
$Log$
+Revision 1.6 2000/11/30 17:16:14 coppedis
+Changes suggested by fca
+
Revision 1.5 2000/11/22 11:30:12 coppedis
Major code revision
fCosy = 0.;
fCosz = 1.;
fPseudoRapidity = 0.;
+
fFermiflag = 1;
// LHC values for beam divergence and crossing angle
fBeamDiv = 0.000032;
fBeamCrossAngle = 0.0001;
fBeamCrossPlane = 2;
+
+ Int_t i, j;
+ for(i=0; i<201; i++){
+ fProbintp[i] = 0;
+ fProbintn[i] = 0;
+ }
+ for(j=0; j<3; j++){
+ fPp[i] = 0;
+ }
+ fDebugOpt = 0;
}
//_____________________________________________________________________________
void AliGenZDC::Init()
{
- printf(" AliGenZDC initialized with:\n");
+ printf("\n\n AliGenZDC initialized with:\n");
printf(" Fermi flag = %d, Beam Divergence = %f, Crossing Angle "
"= %f, Crossing Plane = %d\n\n", fFermiflag, fBeamDiv, fBeamCrossAngle,
fBeamCrossPlane);
+
//Initialize Fermi momentum distributions for Pb-Pb
FermiTwoGaussian(207.,82.,fPp,fProbintp,fProbintn);
}
//
Int_t i;
- Double_t mass, pLab[3], fP0, ddp[3], dddp0, dddp[3];
- Float_t ptot = fPMin;
+ Double_t Mass, pLab[3], fP0, fP[3], fBoostP[3], ddp[3], dddp0, dddp[3];
+ Float_t fPTrack[3], ptot = fPMin;
Int_t nt;
if(fPseudoRapidity==0.){
pLab[2] = ptot*TMath::Cos(scang);
}
for(i=0; i<=2; i++){
- fPInit[i] = pLab[i];
fP[i] = pLab[i];
}
+ if(fDebugOpt == 1){
+ printf("\n\n Initial momentum : \n");
+ for(i=0; i<=2; i++){
+ printf(" p[%d] = %f\n",i,pLab[i]);
+ }
+ }
+
// Beam divergence and crossing angle
if(fBeamCrossAngle!=0.) {
BeamDivCross(1,fBeamDiv,fBeamCrossAngle,fBeamCrossPlane,pLab);
if((fIpart==kProton) || (fIpart==kNeutron)){
ExtractFermi(fIpart,fPp,fProbintp,fProbintn,ddp);
}
- mass=gAlice->PDGDB()->GetParticle(fIpart)->Mass();
- fP0 = TMath::Sqrt(fP[0]*fP[0]+fP[1]*fP[1]+fP[2]*fP[2]+mass*mass);
+ Mass=gAlice->PDGDB()->GetParticle(fIpart)->Mass();
+ fP0 = TMath::Sqrt(fP[0]*fP[0]+fP[1]*fP[1]+fP[2]*fP[2]+Mass*Mass);
for(i=0; i<=2; i++){
dddp[i] = ddp[i];
}
- dddp0 = TMath::Sqrt(dddp[0]*dddp[0]+dddp[1]*dddp[1]+dddp[2]*dddp[2]+mass*mass);
+ dddp0 = TMath::Sqrt(dddp[0]*dddp[0]+dddp[1]*dddp[1]+dddp[2]*dddp[2]+Mass*Mass);
TVector3 b(fP[0]/fP0, fP[1]/fP0, fP[2]/fP0);
TLorentzVector pFermi(dddp[0], dddp[1], dddp[2], dddp0);
}
Float_t polar[3] = {0,0,0};
-// printf("fPTrack = %f, %f, %f \n",fPTrack[0],fPTrack[1],fPTrack[2]);
gAlice->SetTrack(fTrackIt,-1,fIpart,fPTrack,fOrigin.GetArray(),polar,0,
kPPrimary,nt);
+ if(fDebugOpt == 1){
+ printf("\n\n Track momentum:\n");
+ printf("\n fPTrack = %f, %f, %f \n",fPTrack[0],fPTrack[1],fPTrack[2]);
+ }
}
//_____________________________________________________________________________
-void AliGenZDC::FermiTwoGaussian(Double_t A, Float_t Z, Double_t* fPp, Double_t*
- fProbintp, Double_t* fProbintn)
+void AliGenZDC::FermiTwoGaussian(Float_t A, Float_t Z, Double_t *fPp,
+ Double_t *fProbintp, Double_t *fProbintn)
{
//
// Momenta distributions according to the "double-gaussian"
// distribution (Ilinov) - equal for protons and neutrons
//
-// printf(" Initialization of Fermi momenta distribution\n");
+
fProbintp[0] = 0;
fProbintn[0] = 0;
Double_t sig1 = 0.113;
alfa*f2/(TMath::Power(sig2,3.)))*0.005;
fProbintp[i] = fProbintp[i-1] + probp;
fProbintn[i] = fProbintp[i];
-// printf(" fProbintp[%d] = %f, fProbintp[%d] = %f\n",i,fProbintp[i],i,fProbintn[i]);
+ }
+ if(fDebugOpt == 1){
+ printf("\n\n Initialization of Fermi momenta distribution \n");
+ for(Int_t i=0; i<=200; i++){
+ printf(" fProbintp[%d] = %f, fProbintn[%d] = %f\n",i,fProbintp[i],i,fProbintn[i]);
+ }
}
}
//_____________________________________________________________________________
-void AliGenZDC::ExtractFermi(Int_t id, Double_t* fPp, Double_t* fProbintp,
- Double_t* fProbintn, Double_t* ddp)
+void AliGenZDC::ExtractFermi(Int_t id, Double_t *fPp, Double_t *fProbintp,
+ Double_t *fProbintn, Double_t *ddp)
{
//
// Compute Fermi momentum for spectator nucleons
//
-// printf(" Extraction of Fermi momentum\n");
-
+
Int_t i;
Float_t xx = gRandom->Rndm();
assert ( id==kProton || id==kNeutron );
ddp[0] = pext*TMath::Sin(tet)*TMath::Cos(phi);
ddp[1] = pext*TMath::Sin(tet)*TMath::Sin(phi);
ddp[2] = pext*cost;
-// printf(" pFx = %f pFy = %f pFz = %f \n",ddp[0],ddp[1],ddp[2]);
+
+ if(fDebugOpt == 1){
+ printf("\n\n Extraction of Fermi momentum\n");
+ printf("\n pxFermi = %f pyFermi = %f pzFermi = %f \n",ddp[0],ddp[1],ddp[2]);
+ }
}
//_____________________________________________________________________________
void AliGenZDC::BeamDivCross(Int_t icross, Float_t fBeamDiv, Float_t fBeamCrossAngle,
- Int_t fBeamCrossPlane, Double_t* pLab)
+ Int_t fBeamCrossPlane, Double_t *pLab)
{
Double_t tetpart, fipart, tetdiv=0, fidiv=0, angleSum[2], tetsum, fisum;
Double_t rvec;
-// printf(" Beam divergence and crossing angle\n");
Int_t i;
Double_t pmq = 0.;
for(i=0; i<=2; i++){
pLab[0] = pmod*TMath::Sin(tetsum)*TMath::Cos(fisum);
pLab[1] = pmod*TMath::Sin(tetsum)*TMath::Sin(fisum);
pLab[2] = pmod*TMath::Cos(tetsum);
- for(i=0; i<=2; i++){
- fDivP[i] = pLab[i];
+ if(fDebugOpt == 1){
+ printf("\n\n Beam divergence and crossing angle\n");
+ for(i=0; i<=2; i++){
+ printf(" pLab[%d] = %f\n",i,pLab[i]);
+ }
}
}
//_____________________________________________________________________________
void AliGenZDC::AddAngle(Double_t theta1, Double_t phi1, Double_t theta2,
- Double_t phi2, Double_t* angleSum)
+ Double_t phi2, Double_t *angleSum)
{
Double_t temp, conv, cx, cy, cz, ct1, st1, ct2, st2, cp1, sp1, cp2, sp2;
Double_t rtetsum, tetsum, fisum;
AliGenZDC();
AliGenZDC(Int_t npart);
virtual ~AliGenZDC() {}
- virtual void Init();
- virtual void Generate();
+ void Init();
+ void Generate();
+
+ // Fermi smearing, beam divergence and crossing angle
+ void FermiTwoGaussian(Float_t A, Float_t Z, Double_t *pp,
+ Double_t *probintp, Double_t *probintn);
+ void ExtractFermi(Int_t id, Double_t *pp, Double_t *probintp,
+ Double_t *probintn, Double_t *pFermi);
+ void BeamDivCross(Int_t icross, Float_t divergence, Float_t crossangle,
+ Int_t crossplane, Double_t *pLab);
+ void AddAngle(Double_t theta1, Double_t phi1, Double_t theta2,
+ Double_t phi2, Double_t *angle);
+
// Parameters that could be set for generation
- virtual void SetParticle(Int_t ipart) {fIpart=ipart;};
- virtual void SetMomentum(Float_t ptot) {fPMin=ptot; fPMax=ptot;};
- virtual void SetDirection(Float_t zpsrp, Float_t cosx, Float_t cosy, Float_t cosz)
- {fPseudoRapidity=zpsrp; fCosx=cosx; fCosy=cosy; fCosz=cosz;};
- virtual void SetFermi(Int_t Fflag) {fFermiflag=Fflag;};
- virtual void SetDiv(Float_t bmdiv, Float_t bmcra, Int_t iflcr)
- {fBeamDiv=bmdiv; fBeamCrossAngle=bmcra; fBeamCrossPlane=iflcr;};
+ void SetParticle(Int_t ipart) {fIpart=ipart;};
+ void SetMomentum(Float_t ptot) {fPMin=ptot; fPMax=ptot;};
+ void SetDirection(Float_t zpsrp, Float_t cosx, Float_t cosy, Float_t cosz)
+ {fPseudoRapidity=zpsrp; fCosx=cosx; fCosy=cosy; fCosz=cosz;};
+ void SetFermi(Int_t Fflag) {fFermiflag=Fflag;};
+ void SetDiv(Float_t bmdiv, Float_t bmcra, Int_t iflcr)
+ {fBeamDiv=bmdiv; fBeamCrossAngle=bmcra; fBeamCrossPlane=iflcr;};
+ void SetDebug() {fDebugOpt = 1;};
// Getters
Double_t GetFermi2p(Int_t key) {return fProbintp[key];}
Double_t GetFermi2n(Int_t key) {return fProbintn[key];}
- Float_t GetInMomentum(Int_t key) {return fPInit[key];};
- Float_t GetBoostMomentum(Int_t key) {return fBoostP[key];};
- Float_t GetDivMomentum(Int_t key) {return fDivP[key];};
- Float_t GetTrackMomentum(Int_t key) {return fPTrack[key];};
-
- // Fermi smearing, beam divergence and crossing angle
- virtual void FermiTwoGaussian(Double_t A, Float_t Z, Double_t* pp,
- Double_t* probintp, Double_t* probintn);
- virtual void ExtractFermi(Int_t id, Double_t* pp, Double_t* probintp,
- Double_t* probintn, Double_t* pFermi);
- virtual void BeamDivCross(Int_t icross, Float_t divergence, Float_t crossangle,
- Int_t crossplane, Double_t* pLab);
- virtual void AddAngle(Double_t theta1, Double_t phi1, Double_t theta2,
- Double_t phi2, Double_t* angle);
-
+
protected:
- Int_t fIpart; // Particle to generate
- Float_t fCosx; // Cos x of particle
- Float_t fCosy; // Cos y of particle
- Float_t fCosz; // Cos z of particle
- Float_t fPseudoRapidity; // Pseudo Rapidity of particle
- Int_t fFermiflag; // Fermi momentum flag
- Float_t fBeamDiv; // Beam divergence
- Float_t fBeamCrossAngle; // Beam crossing angle
- Int_t fBeamCrossPlane; // Beam crossing plane
- Double_t fProbintp[201]; // for protons
- Double_t fProbintn[201]; // for neutrons
- Double_t fPp[201]; // for protons
- Double_t fP[3]; // temporary momentum
- Float_t fPInit[3]; // initial momentum
- Float_t fBoostP[3]; // boosted momentum
- Float_t fDivP[3]; // divergence
- Float_t fPTrack[3]; // track momentum
+ Int_t fIpart; // Particle to be generated
+ Float_t fCosx; // Director cos of the track - x direction
+ Float_t fCosy; // Director cos of the track - y direction
+ Float_t fCosz; // Director cos of the track - z direction
+ Float_t fPseudoRapidity; // Pseudorapidity (!=0 -> eta of the particle)
+ // (=0 -> director cos of the track)
+ Int_t fFermiflag; // Fermi momentum flag (=1 -> Fermi smearing)
+ Float_t fBeamDiv; // Beam divergence (angle in rad)
+ Float_t fBeamCrossAngle; // Beam crossing angle (angle in rad)
+ Int_t fBeamCrossPlane; // Beam crossing plane
+ // (=1 -> horizontal, =2 -> vertical plane)
+ Double_t fProbintp[201]; // Protons momentum distribution due to Fermi
+ Double_t fProbintn[201]; // Neutrons momentum distribution due to Fermi
+ Double_t fPp[201]; //
+ Int_t fDebugOpt; // Option for debugging
- ClassDef(AliGenZDC,1) // Generator for AliZDC class
+ ClassDef(AliGenZDC,1) // Generator for AliZDC class
};
#endif
{
// Default constructor
+ fDetector = 0;
+ fQuadrant = 0;
+ fADCValue = 0;
}
//____________________________________________________________________________
AliZDCDigit::AliZDCDigit(Int_t Det, Int_t Quad, Float_t ADCValue)
{
- // constructor with all data
-
-
+ // Constructor
+
fDetector = Det;
fQuadrant = Quad;
fADCValue = ADCValue;
//____________________________________________________________________________
AliZDCDigit::AliZDCDigit(const AliZDCDigit & digit)
{
- // copy constructor
+ // Copy constructor
fDetector = digit.fDetector;
fQuadrant = digit.fQuadrant;
return *this ;
}
-// protected:
+ protected:
-// Int_t fNprimary; // Number of primaries
Int_t fDetector; // Detector
Int_t fQuadrant; // Quadrant
Float_t fADCValue; // ADC channel value
// --- Standard libraries
#include <stdlib.h>
-#include <search.h>
// --- ROOT system
#include <TRandom.h>
// --- AliRoot classes
#include "AliZDCFragment.h"
-
-ClassImp(AliZDCFragment)
+ClassImp(AliZDCFragment)
+
+int comp(const void *i,const void *j) {return *(int *)i - *(int *)j;}
+
+
//_____________________________________________________________________________
AliZDCFragment::AliZDCFragment()
{
fNimf = 0;
fZmax = 0;
fTau = 0;
- fNalpha = 0;
- fZtot = 0;
- fNtot = 0;
for(Int_t i=0; i<=99; i++){
fZZ[i] = 0;
fNN[i] = 0;
}
+ fNalpha = 0;
+ fZtot = 0;
+ fNtot = 0;
}
void AliZDCFragment::GenerateIMF(Int_t* fZZ, Int_t &fNalpha)
{
// Coefficients of polynomial for average number of IMF
- Float_t ParamNimf[5]={0.011236,1.8364,56.572,-116.24,58.289};
+ const Float_t ParamNimf[5]={0.011236,1.8364,56.572,-116.24,58.289};
// Coefficients of polynomial for fluctuations on average number of IMF
- Float_t ParamFluctNimf[4]={-0.13176,2.9392,-5.2147,2.3092};
+ const Float_t ParamFluctNimf[4]={-0.13176,2.9392,-5.2147,2.3092};
// Coefficients of polynomial for average maximum Z of fragments
- Float_t ParamZmax[4]={0.16899,14.203,-2.8284,65.036};
+ const Float_t ParamZmax[4]={0.16899,14.203,-2.8284,65.036};
// Coefficients of polynomial for fluctuations on maximum Z of fragments
- Float_t ParamFluctZmax[5]={0.013782,-0.17282,1.5065,1.0654,-2.4317};
+ const Float_t ParamFluctZmax[5]={0.013782,-0.17282,1.5065,1.0654,-2.4317};
// Coefficients of polynomial for exponent tau of fragments Z distribution
- Float_t ParamTau[3]={6.7233,-15.85,13.047};
+ const Float_t ParamTau[3]={6.7233,-15.85,13.047};
//Coefficients of polynomial for average number of alphas
- Float_t ParamNalpha[4]={-0.68554,39.605,-68.311,30.165};
+ const Float_t ParamNalpha[4]={-0.68554,39.605,-68.311,30.165};
// Coefficients of polynomial for fluctuations on average number of alphas
- Float_t ParamFluctNalpha[5]={0.283,6.2141,-17.113,17.394,-6.6084};
+ const Float_t ParamFluctNalpha[5]={0.283,6.2141,-17.113,17.394,-6.6084};
// Coefficients of function for Pb nucleus skin
- Float_t ParamSkinPb[2]={0.93,11.05};
+ const Float_t ParamSkinPb[2]={0.93,11.05};
// Thickness of nuclear surface
- Float_t NuclearThick = 0.52;
+ const Float_t NuclearThick = 0.52;
// Maximum impact parameter for U [r0*A**(1/3)]
- Float_t bMaxU = 14.87;
+ const Float_t bMaxU = 14.87;
// Maximum impact parameter for Pb [r0*A**(1/3)]
- Float_t bMaxPb = 14.22;
+ const Float_t bMaxPb = 14.22;
// Z of the projectile
- Float_t ZProj = 82.;
+ const Float_t ZProj = 82.;
// From b(Pb) to b(U)
Float_t bU = fB*bMaxU/bMaxPb;
Float_t y = gRandom->Rndm();
if(y < ((AverageNimf+FluctNimf)-fNimf)) fNimf += 1;
if(fNimf ==0 && ZbNorm>0.75) fNimf = 1;
-// printf("\n fNimf = %d\n", fNimf);
Float_t AverageZmax = ParamZmax[0]+ParamZmax[1]*ZbNorm+ParamZmax[2]*
TMath::Power(ZbNorm,2)+ParamZmax[3]*TMath::Power(ZbNorm,3);
FluctZmax = FluctZmax*xg;
fZmax = AverageZmax+FluctZmax;
if(fZmax>ZProj) fZmax = ZProj;
-// printf("\n fZmax = %f\n", fZmax);
+
+// printf("\n\n ------------------------------------------------------------");
+// printf("\n Generation of nuclear fragments\n");
+// printf("\n fNimf = %d\n", fNimf);
+// printf("\n fZmax = %f\n", fZmax);
// Find the number of alpha particles
// from Singh et al. : Pb+emulsion
Float_t yy = gRandom->Rndm();
if(yy < ((AverageAlpha+FluctAlpha)-fNalpha)) fNalpha += 1;
- Int_t Choice = 0;
- Float_t ZbFrag = 0, SumZ = 0.;
// 2 possibilities:
// 1) for bNorm < 0.9 ==> first remove alphas, then fragments
// 2) for bNorm > 0.9 ==> first remove fragments, then alphas
+ Int_t Choice = 0;
+ Float_t ZbFrag = 0, SumZ = 0.;
+
if(bNorm<=0.9) {
// remove alpha from zbound to find zbound for fragments (Z>=3)
ZbFrag = fZbAverage-fNalpha*2;
}
delete funTau;
- // Sorting vector in ascending order
-
-// int comp(const void*, const void*);
+ // Sorting vector in ascending order with C function QSORT
qsort((void*)zz,fNimf,sizeof(float),comp);
- for(Int_t i=0; i<fNimf; i++){
+
+// for(Int_t i=0; i<fNimf; i++){
// printf("\n After sorting -> zz[%d] = %f \n",i,zz[i]);
- }
+// }
// Rescale the maximum charge to fZmax
for(Int_t j=0; j<fNimf; j++){
fZZ[j] = Int_t (zz[j]*fZmax/zz[fNimf-1]);
if(fZZ[j]<3) fZZ[j] = 3;
-// printf("\n fZZ[%d] = %d \n",j,fZZ[j]);
+// printf("\n fZZ[%d] = %d \n",j,fZZ[j]);
}
- // Check that the sum of the bound charges is not greater than Zbound-Zalfa
+ // Check that the sum of the bound charges is not > than Zbound-Zalfa
for(Int_t ii=0; ii<fNimf; ii++){
SumZ += fZZ[ii];
for(Int_t i=0; i<fNimf; i++){
SumZ += fZZ[i];
}
-// printf("\n The END -> fNimf = %d, SumZ = %f, fZmax = %f\n",
-// fNimf, SumZ, fZmax);
-// printf("\n fNalpha = %d\n", fNalpha);
-// for(Int_t j=0; j<fNimf; j++){
-// printf("\n fZZ[%d] = %d \n",j,fZZ[j]);
-// }
}
-int comp(const void *i,const void *j){return (int*)j-(int*)i;}
-
//_____________________________________________________________________________
void AliZDCFragment::AttachNeutrons(Int_t *fZZ, Int_t *fNN, Int_t &fZtot,Int_t &fNtot)
{
- Float_t AIon[68]={1.87612,2.80943,3.7284,5.60305,6.53536,
+ const Float_t AIon[68]={1.87612,2.80943,3.7284,5.60305,6.53536,
6.53622,8.39479,9.32699,10.2551,11.17793,
13.04378,14.89917,17.6969,18.62284,21.41483,
22.34193,25.13314,26.06034,28.85188,29.7818,
168.554,171.349,173.4536,177.198,179.0518,
180.675,183.473,188.1345,190.77,193.729,
221.74295};
- Int_t ZIon[68]={1,1,2,3,3,
+ const Int_t ZIon[68]={1,1,2,3,3,
4,4,5,5,6,
7,8,9,10,11,
12,13,14,15,16,
92};
Int_t iZ, iA;
-// printf("\n jfNimf=%d\n",fNimf);
+// printf("\n fNimf=%d\n",fNimf);
+
for(Int_t i=0; i<fNimf; i++) {
for(Int_t j=0; j<68; j++) {
iZ = ZIon[j];
if((fZZ[i]-iZ) == 0){
iA = Int_t(AIon[j]/0.93149432+0.5);
fNN[i] = iA - iZ;
-// printf("\n j=%d,iA=%d,fZZ[%d]=%d,fNN[%d]=%d\n",j,iA,i,fZZ[i],i,fNN[i]);
break;
}
else if((fZZ[i]-iZ) < 0){
fZZ[i] = ZIon[j-1];
iA = Int_t (AIon[j-1]/0.93149432+0.5);
fNN[i] = iA - ZIon[j-1];
-// printf("\n j=%d,iA=%d,fZZ[%d]=%d,fNN[%d]=%d\n",j,iA,i,fZZ[i],i,fNN[i]);
break;
}
}
#include <TMath.h>
-int comp(const void*, const void*);
+extern int comp(const void *, const void *);
class AliZDCFragment : public TNamed {
Float_t fZmax; // Mean value of maximum Z of fragment
Float_t fTau; // Exponent of charge distribution: dN/dZ = Z*exp(-fTau)
Int_t fZZ[100]; // Array of atomic numbers of fragments
- Int_t fNalpha; // Number of alpha particles
-
Int_t fNN[100]; // Array of number of neutrons of fragments
+ Int_t fNalpha; // Number of alpha particles
Int_t fZtot; // Total number of bound protons
Int_t fNtot; // Total number of bound neutrons
// Add a ZDC hit
//
Int_t i;
- for(i=0; i<2; i++) fVolume[i] = vol[i];
- fX = hits[0];
- fY = hits[1];
- fZ = hits[2];
- fPrimKinEn = hits[3];
- fXImpact = hits[4];
- fYImpact = hits[5];
- fSFlag = hits[6];
- fLightPMQ = hits[7];
- fLightPMC = hits[8];
- fEnergy = hits[9];
+ for(i=0; i<2; i++) {
+ fVolume[i] = vol[i];
+ }
+ fX = hits[0];
+ fY = hits[1];
+ fZ = hits[2];
+ fPrimKinEn = hits[3];
+ fXImpact = hits[4];
+ fYImpact = hits[5];
+ fSFlag = hits[6];
+ fLightPMQ = hits[7];
+ fLightPMC = hits[8];
+ fEnergy = hits[9];
}
Int_t operator == (AliZDCHit &quad) {
Int_t i;
// Superfluo finche' c'e' shunt = 1 !?!?
-// if(fTrack!=quad.GetTrack()) return 0;
+ if(fTrack!=quad.GetTrack()) return 0;
for(i=0; i<2; i++) if(fVolume[i]!=quad.GetVolume(i)) return 0;
return 1;
}
/*
$Log$
+Revision 1.15 2001/03/12 17:47:56 hristov
+Changes needed on Sun with CC 5.0
+
Revision 1.14 2001/02/23 16:48:28 coppedis
Correct bug in ZEM hit definition
Some syntax corrections for non standard HP aCC
Revision 1.1 2000/07/10 13:58:01 fca
-New version of ZDC from E.Scomparin & C.Oppedisano
+New version of ZDC from E.Scomparin & C.Oppedisano
Revision 1.7 2000/01/19 17:17:40 fca
fMedSensF2 = 0;
fMedSensZN = 0;
fMedSensZP = 0;
- fMedSensGR = 0;
fMedSensZEM = 0;
+ fMedSensGR = 0;
fMedSensPI = 0;
- fNoShower = 0;
}
//_____________________________________________________________________________
// Standard constructor for Zero Degree Calorimeter
//
- fDigits = new TClonesArray("AliZDCDigit",1000);
-
fMedSensF1 = 0;
fMedSensF2 = 0;
fMedSensZN = 0;
fMedSensZP = 0;
- fMedSensGR = 0;
fMedSensZEM = 0;
+ fMedSensGR = 0;
fMedSensPI = 0;
- fNoShower = 0;
+
+
+ // Parameters for light tables
+ fNalfan = 90; // Number of Alfa (neutrons)
+ fNalfap = 90; // Number of Alfa (protons)
+ fNben = 18; // Number of beta (neutrons)
+ fNbep = 28; // Number of beta (protons)
+ Int_t ip,jp,kp;
+ for(ip=0; ip<4; ip++){
+ for(kp=0; kp<fNalfap; kp++){
+ for(jp=0; jp<fNbep; jp++){
+ fTablep[ip][kp][jp] = 0;
+ }
+ }
+ }
+ Int_t in,jn,kn;
+ for(in=0; in<4; in++){
+ for(kn=0; kn<fNalfan; kn++){
+ for(jn=0; jn<fNben; jn++){
+ fTablen[in][kn][jn] = 0;
+ }
+ }
+ }
+
+ // Parameters for hadronic calorimeters geometry
+ fDimZP[0] = 11.2;
+ fDimZP[1] = 6.;
+ fDimZP[2] = 75.;
+ fPosZN[0] = 0.;
+ fPosZN[1] = 0.;
+ fPosZN[2] = 11650.;
+ fPosZP[0] = -23.;
+ fPosZP[1] = 0.;
+ fPosZP[2] = 11600.;
+ fFibZN[0] = 0.;
+ fFibZN[1] = 0.01825;
+ fFibZN[2] = 50.;
+ fFibZP[0] = 0.;
+ fFibZP[1] = 0.0275;
+ fFibZP[2] = 75.;
+
+ // Parameters for EM calorimeter geometry
+ fPosZEM[0] = 0.;
+ fPosZEM[1] = 5.8;
+ fPosZEM[2] = 11600.;
+
+
+ fDigits = new TClonesArray("AliZDCDigit",1000);
}
//_____________________________________________________________________________
void AliZDCv1::CreateBeamLine()
{
- Float_t angle;
- Float_t zq, conpar[9], elpar[3], tubpar[3];
+ Float_t angle, zq, conpar[9], elpar[3], tubpar[3], zd1, zd2;
Int_t im1, im2;
- Float_t zd1, zd2;
-
Int_t *idtmed = fIdtmed->GetArray();
- // -- Mother of the ZDC
+ // -- Mother of the ZDCs
conpar[0] = 0.;
conpar[1] = 360.;
gMC->Gsvolu("P002", "TUBE", idtmed[5], tubpar, 3);
gMC->Gspos("P002", 1, "ZDC ", 0., 0., tubpar[2] + zd1, 0, "ONLY");
- zd1 += tubpar[2] * 2.;
+ zd1 += tubpar[2]*2.;
tubpar[0] = 10./2.;
tubpar[1] = 10.4/2.;
void AliZDCv1::CreateZDC()
{
- Int_t *idtmed = fIdtmed->GetArray();
Int_t irot1, irot2;
Float_t DimPb[6], DimVoid[6];
+
+ Int_t *idtmed = fIdtmed->GetArray();
+
+ // Parameters for hadronic calorimeters geometry
+ // NB -> parameters used ONLY in CreateZDC()
+ Float_t fDimZN[3] = {3.52, 3.52, 50.}; // Dimensions of neutron detector
+ Float_t fGrvZN[3] = {0.03, 0.03, 50.}; // Grooves for neutron detector
+ Float_t fGrvZP[3] = {0.04, 0.04, 75.}; // Grooves for proton detector
+ Int_t fDivZN[3] = {11, 11, 0}; // Division for neutron detector
+ Int_t fDivZP[3] = {7, 15, 0}; // Division for proton detector
+ Int_t fTowZN[2] = {2, 2}; // Tower for neutron detector
+ Int_t fTowZP[2] = {4, 1}; // Tower for proton detector
+
+ // Parameters for EM calorimeter geometry
+ // NB -> parameters used ONLY in CreateZDC()
+ Float_t fDimZEMPb = 0.15*(TMath::Sqrt(2.)); // z-dimension of the Pb slice
+ Float_t fDimZEMAir = 0.001; // scotch
+ Float_t fFibRadZEM = 0.0315; // External fiber radius (including cladding)
+ Int_t fDivZEM[3] = {92, 0, 20}; // Divisions for EM detector
+ Float_t fDimZEM0 = 2*fDivZEM[2]*(fDimZEMPb+fDimZEMAir+fFibRadZEM*(TMath::Sqrt(2.)));
+ Float_t fDimZEM[6] = {fDimZEM0, 3.5, 3.5, 45., 0., 0.}; // Dimensions of EM detector
+ Float_t fFibZEM2 = fDimZEM[2]/TMath::Sin(fDimZEM[3]*kDegrad)-fFibRadZEM;
+ Float_t fFibZEM[3] = {0., 0.0275, fFibZEM2}; // Fibers for EM calorimeter
//-- Create calorimeters geometry
Int_t *idtmed = fIdtmed->GetArray();
Float_t dens, ubuf[1], wmat[2], a[2], z[2], epsil=0.001, stmin=0.01;
- Int_t i, isvolActive, isvol, inofld;
+ Float_t deemax = -1, stemax;
Float_t fieldm = gAlice->Field()->Max();
Float_t tmaxfd=gAlice->Field()->Max();
+ Int_t i, isvolActive, isvol, inofld;
Int_t isxfld = gAlice->Field()->Integ();
- Float_t deemax=-1;
- Float_t stemax;
// --- Store in UBUF r0 for nuclear radius calculation R=r0*A**1/3
void AliZDCv1::InitTables()
{
Int_t k, j;
- //Initialize parameters for light tables and read them
- fNalfan = 90;
- fNalfap = 90;
- fNben = 18;
- fNbep = 28;
-
+
char *lightfName1,*lightfName2,*lightfName3,*lightfName4,
*lightfName5,*lightfName6,*lightfName7,*lightfName8;
FILE *fp1, *fp2, *fp3, *fp4, *fp5, *fp6, *fp7, *fp8;
printf("\n Digitize -> Det = %d, Quad = %d, Light = %d\n", Det, Quad, Light);
}
+ // Parameters for conversion of light yield in ADC channels
+ Float_t fPMGain[3][5]; // PM gain
+ Float_t fADCRes; // ADC conversion factor
+
Int_t j,i;
for(i=0; i<3; i++){
for(j=0; j<5; j++){
-// fPedMean[i][j] = 50.;
-// fPedSigma[i][j] = 10.;
fPMGain[i][j] = 100000.;
}
}
fADCRes = 0.00000064; // ADC Resolution: 250 fC/ADCch
-// Float_t Ped = gRandom->Gaus(fPedMean[Det-1][Quad],fPedSigma[Det-1][Quad]);
-// Int_t ADCch = Int_t(Light*fPMGain[Det-1][Quad]*fADCRes+Ped);
Int_t ADCch = Int_t(Light*fPMGain[Det-1][Quad]*fADCRes);
-
-// if(fDebug == 1){
-// printf(" Ped = %f, ADCch = %d\n", Ped, ADCch);
-// }
-
+
return ADCch;
}
-//____________________________________________________________________________
-//void AliZDCv1::FinishEvent()
-//{
-// Code moved to Hits2SDigits();
-//}
//_____________________________________________________________________________
void AliZDCv1::SDigits2Digits()
//
Int_t j;
-
Int_t vol[2], ibeta=0, ialfa, ibe, nphe;
Float_t x[3], xdet[3], destep, hits[10], m, ekin, um[3], ud[3], be, radius, out;
TLorentzVector s, p;
if(vol[0]==1){
xdet[0] = x[0]-fPosZN[0];
xdet[1] = x[1]-fPosZN[1];
- if((xdet[0]<=0.) && (xdet[1]>=0.)) vol[1]=1;
- if((xdet[0]>0.) && (xdet[1]>0.)) vol[1]=2;
- if((xdet[0]<0.) && (xdet[1]<0.)) vol[1]=3;
- if((xdet[0]>0.) && (xdet[1]<0.)) vol[1]=4;
+ if((xdet[0]<=0.) && (xdet[1]>=0.)) vol[1]=1;
+ if((xdet[0]>0.) && (xdet[1]>0.)) vol[1]=2;
+ if((xdet[0]<0.) && (xdet[1]<0.)) vol[1]=3;
+ if((xdet[0]>0.) && (xdet[1]<0.)) vol[1]=4;
}
//Quadrant in ZP
if(xqZP>=(i-3) && xqZP<(i-2)){
vol[1] = i;
break;
- }
- }
+ }
+ }
}
//ZEM has only 1 quadrant
vol[1] = 1;
xdet[0] = x[0]-fPosZEM[0];
xdet[1] = x[1]-fPosZEM[1];
-// printf("x %f %f xdet %f %f\n",x[0],x[1],xdet[0],xdet[1]);
}
-// if(vol[1]>4){
-// printf("\n-> Det. %d Quad. %d \n", vol[0], vol[1]);
-// printf("x %f %f xdet %f %f\n",x[0],x[1],xdet[0],xdet[1]);}
// Store impact point and kinetic energy of the ENTERING particle
-// Int_t Curtrack = gAlice->CurrentTrack();
-// Int_t Prim = gAlice->GetPrimary(Curtrack);
-// printf ("Primary: %d, Current Track: %d \n", Prim, Curtrack);
-
// if(Curtrack==Prim){
if(gMC->IsTrackEntering()){
//Particle energy
gMC->TrackMomentum(p);
-// printf("p[0] = %f, p[1] = %f, p[2] = %f, p[3] = %f \n",
-// p[0], p[1], p[2], p[3]);
hits[3] = p[3];
// Impact point on ZDC
gMC->TrackMomentum(p);
Float_t ptot=TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
Float_t beta = ptot/p[3];
-// Int_t pcID = gMC->TrackPid();
-// printf(" Pc %d in quadrant %d -> beta = %f \n", pcID, vol[1], beta);
if(beta<0.67) return;
if((beta>=0.67) && (beta<=0.75)) ibeta = 0;
if((beta>0.75) && (beta<=0.85)) ibeta = 1;
}
// Connect the Root Galice file containing Geometry, Kine, Hits and Digits
- TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("Hijing_b2.root");
+ TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
if (!file) {
printf("\n Creating galice.root \n");
- file = new TFile("Hijing_b2.root");
+ file = new TFile("galice.root");
} else {
printf("\n galice.root found in file list");
}
hPMQ4zp -> SetXTitle("phe");
- TH1F *hEzem = new TH1F("hEzem","Energy deposited in ZEM",100,0,1.e3);
+ TH1F *hEzem = new TH1F("hEzem","Energy deposited in ZEM",100,0,5.e3);
hEzem -> SetXTitle("E (GeV)");
- TH1F *hPMzem = new TH1F("hPMzem","Light produced in ZEM PM",100,0,4.e3);
+ TH1F *hPMzem = new TH1F("hPMzem","Light produced in ZEM PM",100,0,5.e3);
hPMzem -> SetXTitle("phe");
//
//
// Loop over events
//
+// NB -> Il TClonesArray delle particelle va inizializzato prima del loop
+ TParticle *particle;
+ TClonesArray *Particles = gAlice->Particles();
for (Int_t evNumber=0; evNumber<evTot; evNumber++){
Float_t energy, EtotZN=0, EtotZP=0, LightCzn=0, LightCzp=0, LtotZN=0, LtotZP=0;
Int_t nbytes=0, nbytesd=0, ipart, nhits, ndigits, pdgcode, ADCzn, ADCzp;
- TParticle *particle;
AliZDCHit *ZDChit;
AliZDCDigit *ZDCdigit;
// Get pointers to Alice detectors and Hits containers
AliDetector *ZDC = gAlice->GetDetector("ZDC");
- TClonesArray *Particles = gAlice->Particles();
if (ZDC) {
TClonesArray *ZDChits = ZDC->Hits();
- TClonesArray *ZDCdigits = ZDC->Digits();
+// TClonesArray *ZDCdigits = ZDC->Digits();
}
// # of entries in Hits tree
Int_t ntracks = TH->GetEntries();
// # of entries in Digits tree
- TTree *TD = gAlice->TreeD();
- Int_t ndigen = TD->GetEntries();
+// TTree *TD = gAlice->TreeD();
+// Int_t ndigen = TD->GetEntries();
- gAlice->ResetDigits();
- nbytesd += TD->GetEvent(ndigen-1);
+// gAlice->ResetDigits();
+// nbytesd += TD->GetEvent(ndigen-1);
- if (ZDC) {
- ndigits = ZDCdigits->GetEntries();
- printf("\n Digits Tree --- # of entries: %d; # of digits: %d\n",ndigen, ndigits);
- }
- for(Int_t digit=0; digit<ndigits; digit++) {
- ZDCdigit = (AliZDCDigit*)ZDCdigits->UncheckedAt(digit);
- printf("\n Digit# %d, fDetector = %d, fVolume = %d, fADCValue = %f\n",
- digit,ZDCdigit->fDetector,ZDCdigit->fQuadrant,ZDCdigit->fADCValue);
- }
+// if (ZDC) {
+// ndigits = ZDCdigits->GetEntries();
+// printf("\n Digits Tree --- # of entries: %d; # of digits: %d\n",ndigen, ndigits);
+// }
+// for(Int_t digit=0; digit<ndigits; digit++) {
+// ZDCdigit = (AliZDCDigit*)ZDCdigits->UncheckedAt(digit);
+// printf("\n Digit# %d, fDetector = %d, fVolume = %d, fADCValue = %f\n",
+// digit,ZDCdigit->fDetector,ZDCdigit->fQuadrant,ZDCdigit->fADCValue);
+// }
// Start loop on tracks in the hits containers
for (Int_t track=0; track<ntracks; track++) {
}
if(ZDChit->fVolume[0]==3) { //ZEM
hEzem->Fill(ZDChit->fEnergy);
- hPMzem->Fill(ZDChit->fLightPMQ);
+ hPMzem->Fill(ZDChit->fLightPMC);
hspotzem->Fill(ZDChit->fXImpact,ZDChit->fYImpact);
}
hEzp->Fill(EtotZP);
hLzp->Fill(LtotZP);
hPMCzp->Fill(LightCzp);
- printf("\n Histos var -> Ezn = %f, Lzn = %f, Ezp = %f, Lzp = %f \n\n",
- EtotZN, LtotZN, EtotZP, LtotZP);
+// printf("\n Histos var -> Ezn = %f, Lzn = %f, Ezp = %f, Lzp = %f \n\n",
+// EtotZN, LtotZN, EtotZP, LtotZP);
}
}//ZDC
}//Track loop
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff