From d48cca747d0bfce22e6fe6b58d1e577e7467bc88 Mon Sep 17 00:00:00 2001 From: kir Date: Wed, 16 Jul 2003 13:13:03 +0000 Subject: [PATCH] Initial import --- RICH/AliRICHParam.cxx | 55 +++++++ RICH/AliRICHParam.h | 129 ++++++++++++++++ RICH/Makefile | 71 +++++++++ RICH/Opticals.h | 338 ++++++++++++++++++++++++++++++++++++++++++ RICH/menu.C | 230 ++++++++++++++++++++++++++++ 5 files changed, 823 insertions(+) create mode 100644 RICH/AliRICHParam.cxx create mode 100644 RICH/AliRICHParam.h create mode 100644 RICH/Makefile create mode 100644 RICH/Opticals.h create mode 100644 RICH/menu.C diff --git a/RICH/AliRICHParam.cxx b/RICH/AliRICHParam.cxx new file mode 100644 index 00000000000..0b609a45c65 --- /dev/null +++ b/RICH/AliRICHParam.cxx @@ -0,0 +1,55 @@ +#include "AliRICHParam.h" +#include "AliRICHConst.h" //for units ???? + +ClassImp(AliRICHParam) + +//______________________________________________________________________________ +// RICH main parameters manipulator +AliRICHParam::AliRICHParam() +{//defines the default parameters + Segmentation (144,160); //nx,ny + DeadZone (3*cm); //spacer between PC planes + PadSize (8.4*mm,8.0*mm); + Size (80*cm,7*cm,60*cm); //full length, not GEANT half notation + AngleRot (0*deg); //rotation of the whole RICH around Z + AnglesDeg (20,19.5); //XY angle, YZ angle deg + Offset (490*cm+1.267*cm); //1.267???????cm distance from IP to the center of module + GapThickness (8*cm); //Gap Thickness + ProximityGapThickness(0.4*cm); //Proximity Gap Thickness + QuartzLength (133*cm); //Quartz Length + QuartzWidth (127.9*cm); //Quartz Width + QuartzThickness (0.5*cm); //Quartz Thickness + OuterFreonLength (133*cm); //Outer Freon Length + OuterFreonWidth (41.3*cm); //Outer Freon Width + InnerFreonLength (133*cm); //Inner Freon Length + InnerFreonWidth (41.3*cm); //Inner Freon Width + FreonThickness (1.5*cm); //Freon Thickness + RadiatorToPads (0); //Distance from radiator to pads + + SigmaIntegration(5.); + ChargeSlope(27.); + ChargeSpreadX(0.18);ChargeSpreadY(0.18); + MaxAdc(4096); + AlphaFeedback(0.036); + EIonisation(26.e-9); + SqrtKx3(0.77459667); + Kx2(0.962); + Kx4(0.379); + SqrtKy3(0.77459667); + Ky2(0.962); + Ky4(0.379); + Pitch(0.25); + WireSag(1); // 1->On, 0->Off + Voltage(2150); // Should only be 2000, 2050, 2100 or 2150 + + Recalc(); +}//AliRICHParam::named ctor +//______________________________________________________________________________ +void AliRICHParam::Recalc() +{//recalculate + Float_t csi_length=fNy*fPadY+fDeadZone; + Float_t csi_width =fNx*fPadX+2*fDeadZone; + fPadPlaneWidth = (csi_width - 2*fDeadZone)/3; + fPadPlaneLength = (csi_length - fDeadZone)/2; +}//void AliRICHParam::Recalc() +//______________________________________________________________________________ diff --git a/RICH/AliRICHParam.h b/RICH/AliRICHParam.h new file mode 100644 index 00000000000..109eaea7d63 --- /dev/null +++ b/RICH/AliRICHParam.h @@ -0,0 +1,129 @@ +#ifndef AliRICHParam_h +#define AliRICHParam_h + +#include +#include "AliRICHConst.h" + +class AliRICHParam :public TObject +{ +public: + AliRICHParam(); + void Recalc();//Recalculate dependent parameters after changes + inline void Segmentation(Int_t Nx, Int_t Ny) {fNx=Nx;fNy=Ny;Recalc();} + inline Int_t Nx() const{return fNx;} + inline Int_t Ny() const{return fNy;} + inline void DeadZone(Float_t a) { fDeadZone=a;Recalc();} + inline Float_t DeadZone() const{return fDeadZone;} + inline void PadSize(Float_t x,Float_t y) { fPadX=x;fPadY=y;Recalc();} + inline Float_t PadX() const{return fPadX;} + inline Float_t PadY() const{return fPadY;} + inline Float_t PadPlaneWidth() const{return fPadPlaneWidth;} + inline Float_t PadPlaneLength() const{return fPadPlaneLength;} + inline void Size(Float_t x,Float_t y,Float_t z){fSizeX=x;fSizeY=y;fSizeZ=z;} + inline void GeantSize(Float_t *pParam) const{pParam[0]=fSizeX/2;pParam[1]=fSizeY/2;pParam[2]=fSizeZ/2;} + inline Float_t SizeX() const{return fSizeX;} + inline Float_t SizeY() const{return fSizeY;} + inline Float_t SizeZ() const{return fSizeZ;} + inline void Offset(Float_t offset) { fOffset=offset;} + inline Float_t Offset() const{return fOffset;} + inline void AnglesDeg(Float_t xy,Float_t yz) { fAngleXY=xy;fAngleYZ=yz;} + inline Float_t AngleYZ() const{return fAngleYZ*d2r;} + inline Float_t AngleXY() const{return fAngleXY*d2r;} + inline void AngleRot(Float_t angle) { fAngleRot=angle;} + inline Float_t AngleRot() const{return fAngleRot*d2r;} + inline void GapThickness(Float_t a) { fGapThickness=a;} + inline Float_t GapThickness() const{return fGapThickness;} + inline void ProximityGapThickness(Float_t a) { fProximityGapThickness=a;} + inline Float_t ProximityGapThickness() const{return fProximityGapThickness;} + inline void QuartzLength(Float_t a) { fQuartzLength=a;} + inline Float_t QuartzLength() const{return fQuartzLength;} + inline void QuartzWidth(Float_t a) { fQuartzWidth=a;} + inline Float_t QuartzWidth() const{return fQuartzWidth;} + inline void QuartzThickness(Float_t a) { fQuartzThickness=a;} + inline Float_t QuartzThickness() const{return fQuartzThickness;} + inline void OuterFreonLength(Float_t a) { fOuterFreonLength=a;} + inline Float_t OuterFreonLength() const{return fOuterFreonLength;} + inline void OuterFreonWidth(Float_t a) { fOuterFreonWidth=a;} + inline Float_t OuterFreonWidth() const{return fOuterFreonWidth;} + inline void InnerFreonLength(Float_t a) { fInnerFreonLength=a;} + inline Float_t InnerFreonLength() const{return fInnerFreonLength;} + inline void InnerFreonWidth(Float_t a) { fInnerFreonWidth=a;} + inline Float_t InnerFreonWidth() const{return fInnerFreonWidth;} + inline void FreonThickness(Float_t a) { fFreonThickness=a;} + inline Float_t FreonThickness() const{return fFreonThickness;} + inline void RadiatorToPads(Float_t a) { fRadiatorToPads=a;} + inline Float_t RadiatorToPads() const{return fRadiatorToPads;} + + inline void SigmaIntegration(Float_t a) { fSigmaIntegration=a;} + inline Float_t SigmaIntegration() const{return fSigmaIntegration;} + inline void ChargeSpreadX(Float_t a) { fChargeSpreadX=a;} + inline Float_t ChargeSpreadX() const{return fChargeSpreadX;} + inline void ChargeSpreadY(Float_t a) { fChargeSpreadY=a;} + inline Float_t ChargeSpreadY() const{return fChargeSpreadY;} + inline void ChargeSlope(Float_t a) { fChargeSlope=a;} + inline Float_t ChargeSlope() {return fChargeSlope;} + inline void MaxAdc(Float_t a) { fMaxAdc=a;} + inline Float_t MaxAdc() const{return fMaxAdc;} + inline void Pitch(Float_t a) { fPitch=a;}; + inline Float_t Pitch() const{return fPitch;} + inline void AlphaFeedback(Float_t a) { fAlphaFeedback=a;} + inline Float_t AlphaFeedback() const{return fAlphaFeedback;} + inline void EIonisation(Float_t a) { fEIonisation=a;} + inline Float_t EIonisation() const{return fEIonisation;} + inline void SqrtKx3(Float_t a) { fSqrtKx3=a;}; + inline void Kx2(Float_t a) { fKx2=a;}; + inline void Kx4(Float_t a) { fKx4=a;}; + inline void SqrtKy3(Float_t a) { fSqrtKy3=a;}; + inline void Ky2(Float_t a) { fKy2=a;}; + inline void Ky4(Float_t a) { fKy4=a;}; + inline void WireSag(Int_t a) { fWireSag=a;}; + inline void Voltage(Int_t a) { fVoltage=a;}; +protected: + Int_t fNx; //number of pads along X + Int_t fNy; //number of pads along Y + Float_t fDeadZone; //spacer between PC planes, cm + Float_t fPadX; //pad width, cm + Float_t fPadY; //pad lenght, cm + Float_t fPadPlaneWidth; //pad plane width, cm + Float_t fPadPlaneLength; //pad plane length, cm + + Float_t fSizeX; //chamber length, cm + Float_t fSizeY; //chamber thickness, cm + Float_t fSizeZ; //chamber width, cm + Float_t fAngleRot; //azimuthal rotation angle in X-Y plane, grad + Float_t fAngleYZ; //angle between RICH chambers in YZ plane, grad + Float_t fAngleXY; //angle between RICH chambers in XY plane, grad + Float_t fOffset; //chambers offset from IP, cm + Float_t fGapThickness; //gap thickness, cm + Float_t fProximityGapThickness; //proximity gap thickness, cm + Float_t fQuartzLength; //quartz length + Float_t fQuartzWidth; //quartz width + Float_t fQuartzThickness; //quartz thickness + Float_t fOuterFreonLength; //outer freon length + Float_t fOuterFreonWidth; //outer freon width + Float_t fInnerFreonLength; //inner freon length + Float_t fInnerFreonWidth; //inner freon width + Float_t fFreonThickness; //freon thickness + Float_t fRadiatorToPads; //distance from radiator to pads + + Float_t fChargeSlope; //Slope of the charge distribution + Float_t fChargeSpreadX; //Width of the charge distribution in x + Float_t fChargeSpreadY; //Width of the charge distribution in y + Float_t fSigmaIntegration; //Number of sigma's used for charge distribution + Float_t fAlphaFeedback; //Feedback photons coefficient + Float_t fEIonisation; //Mean ionisation energy + Float_t fMaxAdc; //Maximum ADC channel + Float_t fSqrtKx3; //Mathieson parameters for x + Float_t fKx2; //Mathieson parameters for x + Float_t fKx4; //Mathieson parameters for x + Float_t fSqrtKy3; //Mathieson parameters for y + Float_t fKy2; //Mathieson parameters for y + Float_t fKy4; //Mathieson parameters for y + Float_t fPitch; //Anode-cathode pitch + Int_t fWireSag; //Flag to turn on/off (0/1) wire sag + Int_t fVoltage; //Working voltage (2000, 2050, 2100, 2150) + + ClassDef(AliRICHParam,1) //RICH main parameters +}; + +#endif //AliRICHParam_h diff --git a/RICH/Makefile b/RICH/Makefile new file mode 100644 index 00000000000..690b02f9026 --- /dev/null +++ b/RICH/Makefile @@ -0,0 +1,71 @@ +Module :=RICH + +include lib$(Module).pkg + +Sources :=$(SRCS) + +OutputDir :=$(KIR_DIR)/$(Module) +SoLib :=$(KIR_DIR)/lib$(Module).so + +CxxFlags := -g -Wall -fPIC -pipe -I$(ROOTSYS)/include -I$(ALICE_ROOT)/include + +ifdef ALIVERBOSE + Mute := +else + Mute :=@ +endif + +Headers := $(Sources:.cxx=.h) $(Module)LinkDef.h + +DictSource := $(OutputDir)/G__$(Module).cxx +DictHeader := $(DictSource:.cxx=.h) +DictObject := $(DictSource:.cxx=.o) + + +Objects := $(patsubst %.cxx,$(OutputDir)/%.o,$(Sources)) $(DictObject) + + +DepFile := $(OutputDir)/$(Module).d + +##### TARGETS ##### + +$(SoLib): $(Objects) + @echo "Creating $@" + $(Mute)g++ -Wl,-soname,lib$(Module).so -shared -O -g $^ -o $@ + +$(DepFile): $(Headers) + @[ -d $(OutputDir) ] || mkdir -p $(OutputDir) + @touch $(DepFile) + @echo "Generating dependency $@" + $(Mute)rmkdepend -f$(DepFile) -p$(OutputDir)/ -- $(CxxFlags) -- $(Sources) 2>/dev/null + +$(DictSource): $(Headers) + @echo "Generating $@" + $(Mute)rootcint -f $@ -c $(filter -I%,$(CxxFlags)) $^ + +$(DictObject) : $(DictSource) + @echo "Compiling $^" + $(Mute)g++ $(CxxFlags) -I. -c $^ -o $@ + +show: + @echo "Headers: $(Headers)" + @echo "Sources: $(Sources)" + @echo "Depend : $(DepFile)" + @echo "Objects: $(Objects)" + @echo "Library: $(SoLib)" + +spec: $(Sources) + @echo $^ + @echo $@ + +clean: + @echo "Cleaning..." + $(Mute)rm -rf $(Objects) $(DictSource) $(DictHeader) $(SoLib) $(DepFile) + +############################ cxx rule ######################################### +$(OutputDir)/%.o : %.cxx + @echo $*.cxx + $(Mute)g++ $(CxxFlags) -c $*.cxx -o $(OutputDir)/$*.o +############################ Dependencies ##################################### + +-include $(DepFile) diff --git a/RICH/Opticals.h b/RICH/Opticals.h new file mode 100644 index 00000000000..ca774539f39 --- /dev/null +++ b/RICH/Opticals.h @@ -0,0 +1,338 @@ +#ifdef __CINT__ +void Opticals() +{ + gROOT->Reset(); +#endif + int i; + const Int_t kNbins=26; + Float_t aPckov[kNbins]; + for(i=0;iDivide(2,2); + + pC->cd(1); + TGraph *pAbsFreonGr=new TGraph(kNbins,aPckov,aAbsFreon); + pAbsFreonGr->SetMarkerStyle(kFreonMarker); pAbsFreonGr->SetMarkerColor(kFreonColor); + TGraph *pAbsSiO2Gr=new TGraph(kNbins,aPckov,aAbsQuartz); + pAbsSiO2Gr->SetMarkerStyle(kSiO2Marker); pAbsSiO2Gr->SetMarkerColor(kSiO2Color); + TMultiGraph *pAbsMG=new TMultiGraph(); + TLegend *pAbsLegend=new TLegend(0.6,0.3,0.85,0.5); + pAbsMG->Add(pAbsFreonGr); pAbsLegend->AddEntry(pAbsFreonGr, "freon","p"); //1 + pAbsMG->Add(pAbsSiO2Gr); pAbsLegend->AddEntry(pAbsSiO2Gr, "quartz","p"); //2 + pAbsMG->Draw("APL"); + pAbsMG->GetXaxis()->SetTitle("energy, GeV"); + pAbsMG->GetYaxis()->SetTitle("absorption length, cm"); + pAbsMG->Draw("APL"); + pAbsLegend->Draw(); + + + pC->cd(2); + TGraph *pIndexFreonGr=new TGraph(kNbins,aPckov,aIndexFreon); + pIndexFreonGr->SetMarkerStyle(kFreonMarker); pIndexFreonGr->SetMarkerColor(kFreonColor); + TGraph *pIndexSiO2Gr=new TGraph(kNbins,aPckov,aIndexQuartz); + pIndexSiO2Gr->SetMarkerStyle(kSiO2Marker); pIndexSiO2Gr->SetMarkerColor(kSiO2Color); + TGraph *pIndexCH4Gr=new TGraph(kNbins,aPckov,aIndexCH4); + pIndexCH4Gr->SetMarkerStyle(kCH4Marker); pIndexCH4Gr->SetMarkerColor(kCH4Color); + TMultiGraph *pIndexMG=new TMultiGraph(); + TLegend *pIndexLegend=new TLegend(0.6,0.2,0.85,0.4); + pIndexMG->Add(pIndexFreonGr); pIndexLegend->AddEntry(pIndexFreonGr, "freon","p"); //1 + pIndexMG->Add(pIndexSiO2Gr); pIndexLegend->AddEntry(pIndexSiO2Gr, "quartz","p"); + pIndexMG->Add(pIndexCH4Gr); pIndexLegend->AddEntry(pIndexCH4Gr, "CH4","p"); + pIndexMG->Draw("APL"); + pIndexMG->GetXaxis()->SetTitle("energy, GeV"); + pIndexMG->GetYaxis()->SetTitle("refraction index"); + pIndexMG->Draw("APL"); + pIndexLegend->Draw(); + + pC->cd(3); + TGraph *pAbsCH4Gr=new TGraph(kNbins,aPckov,aAbsCH4); + pAbsCH4Gr->SetMarkerStyle(kCH4Marker); pAbsCH4Gr->SetMarkerColor(kCH4Color); + pAbsCH4Gr->Draw("APL"); + pAbsCH4Gr->GetXaxis()->SetTitle("energy, GeV"); + pAbsCH4Gr->SetTitle("CH4 absorption length, cm"); + pAbsCH4Gr->Draw("APL"); + + + pC->cd(4); + TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI); + pQeCsIGr->SetMarkerStyle(kCsIMarker); pQeCsIGr->SetMarkerColor(kCsIColor); + pQeCsIGr->Draw("APL"); + pQeCsIGr->GetXaxis()->SetTitle("energy, GeV"); + pQeCsIGr->SetTitle("CsI QE, all others are 1"); + pQeCsIGr->Draw("APL"); +}//main + +Float_t AbsoCH4(Float_t x) +{ + + //KLOSCH,SCH4(9),WL(9),EM(9),ALENGTH(31) + Float_t sch4[9] = {.12,.16,.23,.38,.86,2.8,7.9,28.,80.}; //MB X 10^22 + //Float_t wl[9] = {153.,152.,151.,150.,149.,148.,147.,146.,145}; + Float_t em[9] = {8.1,8.158,8.212,8.267,8.322,8.378,8.435,8.493,8.55}; + const Float_t kLosch=2.686763E19; // LOSCHMIDT NUMBER IN CM-3 + const Float_t kIgas1=100, kIgas2=0, kOxy=10., kWater=5., kPressure=750.,kTemperature=283.; + Float_t pn=kPressure/760.; + Float_t tn=kTemperature/273.16; + + +// ------- METHANE CROSS SECTION ----------------- +// ASTROPH. J. 214, L47 (1978) + + Float_t sm=0; + if (x<7.75) + sm=.06e-22; + + if(x>=7.75 && x<=8.1) + { + Float_t c0=-1.655279e-1; + Float_t c1=6.307392e-2; + Float_t c2=-8.011441e-3; + Float_t c3=3.392126e-4; + sm=(c0+c1*x+c2*x*x+c3*x*x*x)*1.e-18; + } + + if (x> 8.1) + { + Int_t j=0; + while (x<=em[j] && x>=em[j+1]) + { + j++; + Float_t a=(sch4[j+1]-sch4[j])/(em[j+1]-em[j]); + sm=(sch4[j]+a*(x-em[j]))*1e-22; + } + } + + Float_t dm=(kIgas1/100.)*(1.-((kOxy+kWater)/1.e6))*kLosch*pn/tn; + Float_t abslm=1./sm/dm; + +// ------- ISOBUTHANE CROSS SECTION -------------- +// i-C4H10 (ai) abs. length from curves in +// Lu-McDonald paper for BARI RICH workshop . +// ----------------------------------------------------------- + + Float_t ai; + Float_t absli; + if (kIgas2 != 0) + { + if (x<7.25) + ai=100000000.; + + if(x>=7.25 && x<7.375) + ai=24.3; + + if(x>=7.375) + ai=.0000000001; + + Float_t si = 1./(ai*kLosch*273.16/293.); // ISOB. CRO.SEC.IN CM2 + Float_t di=(kIgas2/100.)*(1.-((kOxy+kWater)/1.e6))*kLosch*pn/tn; + absli =1./si/di; + } + else + absli=1.e18; +// --------------------------------------------------------- +// +// transmission of O2 +// +// y= path in cm, x=energy in eV +// so= cross section for UV absorption in cm2 +// do= O2 molecular density in cm-3 +// --------------------------------------------------------- + + Float_t abslo; + Float_t so=0; + if(x>=6.0) + { + if(x>=6.0 && x<6.5) + { + so=3.392709e-13 * TMath::Exp(2.864104 *x); + so=so*1e-18; + } + + if(x>=6.5 && x<7.0) + { + so=2.910039e-34 * TMath::Exp(10.3337*x); + so=so*1e-18; + } + + + if (x>=7.0) + { + Float_t a0=-73770.76; + Float_t a1=46190.69; + Float_t a2=-11475.44; + Float_t a3=1412.611; + Float_t a4=-86.07027; + Float_t a5=2.074234; + so= a0+(a1*x)+(a2*x*x)+(a3*x*x*x)+(a4*x*x*x*x)+(a5*x*x*x*x*x); + so=so*1e-18; + } + + Float_t dox=(kOxy/1e6)*kLosch*pn/tn; + abslo=1./so/dox; + } + else + abslo=1.e18; +// --------------------------------------------------------- +// +// transmission of H2O +// +// y= path in cm, x=energy in eV +// sw= cross section for UV absorption in cm2 +// dw= H2O molecular density in cm-3 +// --------------------------------------------------------- + + Float_t abslw; + + Float_t b0=29231.65; + Float_t b1=-15807.74; + Float_t b2=3192.926; + Float_t b3=-285.4809; + Float_t b4=9.533944; + + if(x>6.75) + { + Float_t sw= b0+(b1*x)+(b2*x*x)+(b3*x*x*x)+(b4*x*x*x*x); + sw=sw*1e-18; + Float_t dw=(kWater/1e6)*kLosch*pn/tn; + abslw=1./sw/dw; + } + else + abslw=1.e18; + +// --------------------------------------------------------- + + Float_t alength=1./(1./abslm+1./absli+1./abslo+1./abslw); + return (alength); +}//AbsoCH4 + + +Float_t Fresnel(Float_t ene,Float_t pdoti, Bool_t pola) +{ + + //ENE(EV), PDOTI=COS(INC.ANG.), PDOTR=COS(POL.PLANE ROT.ANG.) + + Float_t en[36] = {5.0,5.1,5.2,5.3,5.4,5.5,5.6,5.7,5.8,5.9,6.0,6.1,6.2, + 6.3,6.4,6.5,6.6,6.7,6.8,6.9,7.0,7.1,7.2,7.3,7.4,7.5,7.6,7.7, + 7.8,7.9,8.0,8.1,8.2,8.3,8.4,8.5}; + + + Float_t csin[36] = {2.14,2.21,2.33,2.48,2.76,2.97,2.99,2.59,2.81,3.05, + 2.86,2.53,2.55,2.66,2.79,2.96,3.18,3.05,2.84,2.81,2.38,2.11, + 2.01,2.13,2.39,2.73,3.08,3.15,2.95,2.73,2.56,2.41,2.12,1.95, + 1.72,1.53}; + + Float_t csik[36] = {0.,0.,0.,0.,0.,0.196,0.408,0.208,0.118,0.49,0.784,0.543, + 0.424,0.404,0.371,0.514,0.922,1.102,1.139,1.376,1.461,1.253,0.878, + 0.69,0.612,0.649,0.824,1.347,1.571,1.678,1.763,1.857,1.824,1.824, + 1.714,1.498}; + Float_t xe=ene; + Int_t j=Int_t(xe*10)-49; + Float_t cn=csin[j]+((csin[j+1]-csin[j])/0.1)*(xe-en[j]); + Float_t ck=csik[j]+((csik[j+1]-csik[j])/0.1)*(xe-en[j]); + + //FORMULAE FROM HANDBOOK OF OPTICS, 33.23 OR + //W.R. HUNTER, J.O.S.A. 54 (1964),15 , J.O.S.A. 55(1965),1197 + + Float_t sinin=TMath::Sqrt(1-pdoti*pdoti); + Float_t tanin=sinin/pdoti; + + Float_t c1=cn*cn-ck*ck-sinin*sinin; + Float_t c2=4*cn*cn*ck*ck; + Float_t aO=TMath::Sqrt(0.5*(TMath::Sqrt(c1*c1+c2)+c1)); + Float_t b2=0.5*(TMath::Sqrt(c1*c1+c2)-c1); + + Float_t rs=((aO-pdoti)*(aO-pdoti)+b2)/((aO+pdoti)*(aO+pdoti)+b2); + Float_t rp=rs*((aO-sinin*tanin)*(aO-sinin*tanin)+b2)/((aO+sinin*tanin)*(aO+sinin*tanin)+b2); + + + //CORRECTION FACTOR FOR SURFACE ROUGHNESS + //B.J. STAGG APPLIED OPTICS, 30(1991),4113 + + Float_t sigraf=18.; + Float_t lamb=1240/ene; + Float_t fresn; + + Float_t rO=TMath::Exp(-(4*TMath::Pi()*pdoti*sigraf/lamb)*(4*TMath::Pi()*pdoti*sigraf/lamb)); + + if(pola) + { + Float_t pdotr=0.8; //DEGREE OF POLARIZATION : 1->P , -1->S + fresn=0.5*(rp*(1+pdotr)+rs*(1-pdotr)); + } + else + fresn=0.5*(rp+rs); + + fresn = fresn*rO; + return(fresn); +}//Fresnel(...) +#endif diff --git a/RICH/menu.C b/RICH/menu.C new file mode 100644 index 00000000000..865ac75168f --- /dev/null +++ b/RICH/menu.C @@ -0,0 +1,230 @@ +void Show() +{ +//How to get number of events: + gpRL->LoadHeader();//without this the following will be zero: + Info("RICH/menu.C::Show","3 ways to get number of events 1=%i 2=%i 3=%i", + gAlice->GetEventsPerRun(), + gpRL->TreeE()->GetEntries(), + gAlice->TreeE()->GetEntries()); +// + return; + for(Int_t iEventN=0;iEventNGetEvent(iEventN); + Int_t iNtracks=gAlice->TreeH()->GetEntries(); + cout<<"Event "<TreeH()->GetEntry(iTrackN); + // gAlice->Particle(iTrackN)->Print(); + for(AliRICHhit *pRICHhit=(AliRICHhit*)pRICH->FirstHit(-1);pRICHhit;pRICHhit->(AliRICHhit*)pRICH->NextHit()){ + TVector3 mrsV3(pRICHhit->X(),pRICHhit->Y(),pRICHhit->Z()); + cout<<"Before\n";mrsV3.Dump(); + TVector3 armV3=pRICH->ToArm(mrsV3); + cout<<"After\n";armV3.Dump(); + }//loop on hits of given track + }// loop on tracks + }//loop on events + +}//void Show() + +void menu(Int_t iNevents=5)// How many events to generate. +{ + Info("RICH/menu.C","%i event(s) are requested",iNevents); + + TString runString="gAlice->Run("; runString=runString+iNevents+");"; + + TControlBar *pMenu = new TControlBar("vertical","RICH main"); + + if(GetAlice()){//it's from file, reconstruct + pMenu->AddButton("Show","Show()","Shows the structure of events in files"); + pMenu->AddButton("RingViewer","RingViewer()","Show rings with reconstructed info"); + }else{//it's aliroot, simulate + pMenu->AddButton("Init", "Init()", "Normal init"); + pMenu->AddButton("Run", runString.Data(), "Process!"); + } + pMenu->AddButton("Debug ON", "gAlice->SetDebug(1)", "Switch debug on"); + pMenu->AddButton("Debug OFF", "gAlice->SetDebug()", "Switch debug off"); + pMenu->AddButton("Geo", "Geo()", "Geomentry submenu"); + pMenu->AddButton("Browser", "new TBrowser;", "Start ROOT TBrowser"); + pMenu->AddButton("Quit", ".q", "Close session"); + pMenu->Show(); +}//void menu(Int_t iNevents) + +void Init() +{ + gAlice->Init("~/my/AliceConfig.C"); +} + +AliRunLoader *rl,*gpRL; +AliLoader *rrl,*gpRichLoader; +AliRICH *r,*gpRich; +AliRun *a; +Bool_t GetAlice() +{ + if(gAlice){//it's aliroot + Info("RICH/menu.C::GetAlice","gAlice!=NULL, it's aliroot, execute Init()"); + Init(); + return kFALSE; + }else{//it's root with ALICE libs loaded + Info("RICH/menu.C::GetAlice","gAlice=0, getting it from simulated file."); + + gpRL=AliRunLoader::Open("galice.root","AlicE","update"); + if(!gpRL) Fatal("GetAlice","Can't get AliRunLoader"); + + gpRL->LoadgAlice(); + if(!gAlice) Fatal("RICH/menu.C::GetAlice","No gAlice in file"); +// gAlice=gpRL->GetAliRun(); + + gpRich=(AliRICH*)gAlice->GetDetector("RICH"); + if(!gpRich) Fatal("RICH/menu.C::GetAlice","No RICH in file"); + + gpRichLoader=gpRL->GetLoader("RICHLoader"); + if(!gpRichLoader) Fatal("RICH/menu.C::GetAlice","No RICH loader in file"); + + Info("RICH/menu.C::GetAlice","contains %i event(s)",gAlice->GetEventsPerRun()); + rl=gpRL;rrl=gpRichLoader;r=gpRich;a=gAlice;//for manual manipulation convinience + return kTRUE; + } +}//void GetAlice() + + + +void Geo() +{ + TControlBar *pMenu = new TControlBar("vertical","RICH draw"); + pMenu->AddButton("Draw Front", "gMC->Gdraw(\"ALIC\", 0,0,0, 10,10, 0.01,0.01)","Draws ALIC volume in XY view"); + pMenu->AddButton("Draw Side", "gMC->Gdraw(\"ALIC\",90,180, 0, 10,10, 0.01,0.01)","Draws ALIC volume in YZ view"); + pMenu->AddButton("Draw Top", "gMC->Gdraw(\"ALIC\",90, 90, 0, 10,10, 0.03,0.03)","Draws ALIC volume in XZ view"); + pMenu->AddButton("ALICE Tree", "((TGeant3*)gMC)->Gdtree(\"ALIC\")","Draws ALICE tree"); + pMenu->AddButton("RICH Tree", "((TGeant3*)gMC)->Gdtree(\"RICH\")","Draws RICH tree"); + pMenu->AddButton("Geo test", "GeoTest()", "Draw RICH geo as a macro"); + pMenu->AddButton("Print ref", "PrintGeo()", "Print RICH chambers default position"); + pMenu->AddButton("Print act", "r->Print();", "Print RICH chambers default position"); + pMenu->Show(); +}//void Draw() + +void GeoTest() +{ + + TBRIK *pAliceBRIK=new TBRIK("aliceBRIK","ALICE mother volume","void",500,500,500); + AliRICH *pRICH=(AliRICH*)gAlice->GetDetector("RICH"); + TBRIK *pArmBRIK=new TBRIK("armBRIK","RICH arm1","void",pRICH->GetSizeX(),pRICH->GetSizeY(),pRICH->GetSizeZ()); + + TNode *pAliceNode=new TNode("aliceNode","Mother volume","aliceBRIK"); + pAliceNode->cd(); + +// ARM 1 LEFT + TRotation rot1; + TVector3 v1(0,pRICH->GetOffset(),0); + + + rot1.RotateX(pRICH->GetYZAngle()*kDegrad); v1.RotateX(pRICH->GetYZAngle()*kDegrad); + rot1.RotateZ(pRICH->GetXYAngle()*kDegrad); v1.RotateZ(pRICH->GetXYAngle()*kDegrad); + rot1.RotateZ(pRICH->GetRotAngle()*kDegrad); v1.RotateZ(pRICH->GetRotAngle()*kDegrad); + + TRotMatrix *pArm1RotMatrix=new TRotMatrix("rotArm1","rotArm1",rot1.ThetaX()*kRaddeg, rot1.PhiX()*kRaddeg, + rot1.ThetaY()*kRaddeg, rot1.PhiY()*kRaddeg, + rot1.ThetaZ()*kRaddeg, rot1.PhiZ()*kRaddeg); + + TNode *pArm1Node=new TNode("arm1Node","Left arm","armBRIK",v1.X(),v1.Y(),v1.Z(),"rotArm1"); + arm1Node->SetLineColor(kRed); + +// ARM 2 LEFT + TRotation rot2; + TVector3 v2(0,pRICH->GetOffset(),0); + + + rot2.RotateX(pRICH->GetYZAngle()*kDegrad); v2.RotateX(pRICH->GetYZAngle()*kDegrad); + rot2.RotateZ(-pRICH->GetXYAngle()*kDegrad); v2.RotateZ(-pRICH->GetXYAngle()*kDegrad); + rot2.RotateZ(pRICH->GetRotAngle()*kDegrad); v2.RotateZ(pRICH->GetRotAngle()*kDegrad); + + TRotMatrix *pArm2RotMatrix=new TRotMatrix("rotArm2","rotArm2",rot2.ThetaX()*kRaddeg, rot2.PhiX()*kRaddeg, + rot2.ThetaY()*kRaddeg, rot2.PhiY()*kRaddeg, + rot2.ThetaZ()*kRaddeg, rot2.PhiZ()*kRaddeg); + + TNode *pArm2Node=new TNode("arm2Node","Left arm","armBRIK",v2.X(),v2.Y(),v2.Z(),"rotArm2"); + arm2Node->SetLineColor(kBlue); + + aliceNode->Draw(); +}//void GeoTest() + +void RingViewer() +{ + gStyle->SetPalette(1); + TCanvas *view = new TCanvas("Display","ALICE RICH Display",0,0,1200,750); + + TH2F *p2F = new TH2F("p2F","RICH DISPLAY",160,0,160,144,0,144); + p2F->SetStats(0); + p2F->SetMaximum(100); + + Int_t Nevents = gAlice->GetEventsPerRun(); +} + +void PrintGeo(Float_t rotDeg=0) +{ + AliRICHParam *p=new AliRICHParam; + Double_t r=p->Offset(); + Double_t kP=p->AngleXY(); + Double_t kT=p->AngleYZ(); + Double_t kRot; + + if(rotDeg==0) + kRot=p->AngleRot(); + else + kRot=rotDeg*deg; + + cout<