Initial import
authorkir <kir@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 16 Jul 2003 13:13:03 +0000 (13:13 +0000)
committerkir <kir@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 16 Jul 2003 13:13:03 +0000 (13:13 +0000)
RICH/AliRICHParam.cxx [new file with mode: 0644]
RICH/AliRICHParam.h [new file with mode: 0644]
RICH/Makefile [new file with mode: 0644]
RICH/Opticals.h [new file with mode: 0644]
RICH/menu.C [new file with mode: 0644]

diff --git a/RICH/AliRICHParam.cxx b/RICH/AliRICHParam.cxx
new file mode 100644 (file)
index 0000000..0b609a4
--- /dev/null
@@ -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 (file)
index 0000000..109eaea
--- /dev/null
@@ -0,0 +1,129 @@
+#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
diff --git a/RICH/Makefile b/RICH/Makefile
new file mode 100644 (file)
index 0000000..690b02f
--- /dev/null
@@ -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 (file)
index 0000000..ca77453
--- /dev/null
@@ -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;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
diff --git a/RICH/menu.C b/RICH/menu.C
new file mode 100644 (file)
index 0000000..865ac75
--- /dev/null
@@ -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;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()