--- /dev/null
+#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()
+//______________________________________________________________________________
--- /dev/null
+#ifndef AliRICHParam_h
+#define AliRICHParam_h
+
+#include <TObject.h>
+#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
--- /dev/null
+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)
--- /dev/null
+#ifdef __CINT__
+void Opticals()
+{
+ gROOT->Reset();
+#endif
+ int i;
+ const Int_t kNbins=26;
+ Float_t aPckov[kNbins];
+ for(i=0;i<kNbins;i++){ //Photons energy intervals
+ aPckov[i]=(Float_t(i)*0.1+5.5)*1e-9;
+ }
+
+ Float_t aIndexFreon[kNbins];
+ Float_t aIndexQuartz[kNbins];
+ Float_t aIndexOpaqueQuartz[kNbins];
+ Float_t aIndexCH4[kNbins];
+ Float_t aIndexGrid[kNbins];
+
+ Float_t e1= 10.666;Float_t e2= 18.125; Float_t f1= 46.411; Float_t f2= 228.71;
+ for (i=0;i<kNbins;i++){
+ aIndexFreon[i] = aPckov[i] * .0172 * 1e9 + 1.177;
+ Float_t ene=aPckov[i]*1e9;
+ Float_t a=f1/(e1*e1 - ene*ene);
+ Float_t b=f2/(e2*e2 - ene*ene);
+ aIndexQuartz[i] = TMath::Sqrt(1. + a + b );
+ aIndexOpaqueQuartz[i] =1;
+ aIndexCH4[i] =1.000444;
+ aIndexGrid[i] =1;
+ }
+
+ Float_t aAbsFreon[kNbins]={179.0987, 179.0987, 179.0987, 179.0987, 179.0987,
+ 179.0987, 179.0987, 179.0987, 179.0987, 142.9206,
+ 56.64957, 25.58622, 13.95293, 12.03905, 10.42953,
+ 8.804196, 7.069031, 4.461292, 2.028366, 1.293013,
+ 0.577267, 0.40746, 0.334964, 0.0, 0.0,
+ 0};
+
+
+ Float_t aAbsQuartz[kNbins]={105.8, 65.52, 48.58, 42.85, 35.79,
+ 31.262, 28.598, 27.527, 25.007, 22.815,
+ 21.004, 19.266, 17.525, 15.878, 14.177,
+ 11.719, 9.282, 6.62, 4.0925, 2.601,
+ 1.149, 0.667, 0.3627, 0.192, 0.1497,
+ 0.10857};
+
+ Float_t aAbsOpaqueQuartz[kNbins];
+ Float_t aAbsCH4[kNbins];
+ Float_t aAbsGrid[kNbins];
+ Float_t aAbsCsI[kNbins];
+ for (i=0;i<kNbins;i++){
+ aAbsCH4[i] =AbsoCH4(aPckov[i]*1e9);
+ aAbsOpaqueQuartz[i]=1e-5;
+ aAbsCsI[i] =1e-4;
+ aAbsGrid[i] =1e-4;
+ }
+
+ Float_t aQeCsI[kNbins] = {0.000199999995, 0.000600000028, 0.000699999975, 0.00499999989, 0.00749999983,
+ 0.010125, 0.0242999997, 0.0405000001, 0.0688500032, 0.105299994,
+ 0.121500008, 0.141749993, 0.157949999, 0.162, 0.166050002,
+ 0.167669997, 0.174299985, 0.176789999, 0.179279998, 0.182599992,
+ 0.18592, 0.187579989, 0.189239994, 0.190899998, 0.207499996,
+ 0.215799987};
+ Float_t aQeAll[kNbins];
+ for(i=0;i<kNbins;i++){
+ aQeCsI[i]/= (1.0-Fresnel(aPckov[i]*1e9,1.0,0)); //FRESNEL LOSS CORRECTION
+ aQeAll[i]=1; //QE for all other materials except for PC must be 1.
+ }
+
+
+#ifdef __CINT__
+
+
+//Freon, Quartz, Opaque Quartz,Methane,CsI,Grid
+ const Int_t kFreonMarker= 24; const Int_t kFreonColor=kRed;
+ const Int_t kCH4Marker= 25; const Int_t kCH4Color=kGreen;
+ const Int_t kSiO2Marker= 26; const Int_t kSiO2Color=kBlue;
+ const Int_t kCsIMarker = 2; const Int_t kCsIColor =kMagenta;
+
+ TCanvas *pC=new TCanvas("c1","RICH optics to check",800,900);
+
+ pC->Divide(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
--- /dev/null
+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;iEventN<iNevents;iEventN++){// loop on events
+ Int_t iNparticles=gAlice->GetEvent(iEventN);
+ Int_t iNtracks=gAlice->TreeH()->GetEntries();
+ cout<<"Event "<<iEventN<<" contains "<<iNparticles<<" particles and "<<iNtracks<<" tracks\n";
+
+ for(Int_t iTrackN=0;iTrackN<iNtracks;iTrackN++){ // loop on tracks
+ gAlice->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<<endl;
+ Double_t phi=90*deg+kRot,theta=90*deg-kT;
+ Info(" menu for 0","r=%8.3f theta=%5.1f phi=%5.1f x=%8.3f y=%8.3f z=%8.3f",
+ r, theta*r2g, phi*r2g,
+ r*sin(theta)*cos(phi),
+ r*sin(theta)*sin(phi),
+ r*cos(theta));
+
+ phi=90*deg+kRot+kP,theta=90*deg;
+ Info(" menu for 1","r=%8.3f theta=%5.1f phi=%5.1f x=%8.3f y=%8.3f z=%8.3f",
+ r, theta*r2g, phi*r2g,
+ r*sin(theta)*cos(phi),
+ r*sin(theta)*sin(phi),
+ r*cos(theta));
+
+ phi=90*deg+kRot,theta=90*deg;
+ Info(" menu for 2","r=%8.3f theta=%5.1f phi=%5.1f x=%8.3f y=%8.3f z=%8.3f",
+ r, theta*r2g, phi*r2g,
+ r*sin(theta)*cos(phi),
+ r*sin(theta)*sin(phi),
+ r*cos(theta));
+
+
+ phi=90*deg+kRot-kP,theta=90*deg;
+ Info(" menu for 3","r=%8.3f theta=%5.1f phi=%5.1f x=%8.3f y=%8.3f z=%8.3f",
+ r, theta*r2g, phi*r2g,
+ r*sin(theta)*cos(phi),
+ r*sin(theta)*sin(phi),
+ r*cos(theta));
+
+
+ phi=90*deg+kRot+kP,theta=90*deg+kT;
+ Info(" menu for 4","r=%8.3f theta=%5.1f phi=%5.1f x=%8.3f y=%8.3f z=%8.3f",
+ r, theta*r2g, phi*r2g,
+ r*sin(theta)*cos(phi),
+ r*sin(theta)*sin(phi),
+ r*cos(theta));
+
+
+ phi=90*deg+kRot,theta=90*deg+kT;
+ Info(" menu for 5","r=%8.3f theta=%5.1f phi=%5.1f x=%8.3f y=%8.3f z=%8.3f",
+ r, theta*r2g, phi*r2g,
+ r*sin(theta)*cos(phi),
+ r*sin(theta)*sin(phi),
+ r*cos(theta));
+
+ phi=90*deg+kRot-kP,theta=90*deg+kT;
+ Info(" menu for 6","r=%8.3f theta=%5.1f phi=%5.1f x=%8.3f y=%8.3f z=%8.3f",
+ r, theta*r2g, phi*r2g,
+ r*sin(theta)*cos(phi),
+ r*sin(theta)*sin(phi),
+ r*cos(theta));
+
+ delete p;
+}//void PrintGeo()
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff