No return in void function (Alpha)
[u/mrichter/AliRoot.git] / RICH / Opticals.h
1 #ifdef __CINT__
2 void Opticals()
3 {
4   gROOT->Reset();
5 #endif   
6   int i;
7   const Int_t kNbins=30;//number of photon energy points
8   
9   Float_t aPckov[kNbins];  for(i=0;i<kNbins;i++) aPckov[i]=1e-9*AliRICHParam::PhotonEnergy(i);
10     
11   
12   Float_t aIdxC6F14[kNbins];  //Freon ref index
13   Float_t aIdxCH4[kNbins];    //Methane ref index
14   Float_t aIdxSiO2[kNbins];   //Quartz ref index
15   Float_t aIdxMetal[kNbins];  //Metal ref index always 0
16   Float_t aIdxGel[kNbins];    //Aerogel ref index
17         
18   for (i=0;i<kNbins;i++){
19     aIdxC6F14[i] = AliRICHParam::RefIdxC6F14(aPckov[i]*1e9);
20     aIdxSiO2[i]  = AliRICHParam::RefIdxSiO2(aPckov[i]*1e9);
21     aIdxCH4[i]   = AliRICHParam::RefIdxCH4(5.5);
22     aIdxGel[i]   = AliRICHParam::RefIdxGel(5.5);
23     aIdxMetal[i]  =0;//metal
24   } 
25     
26   Float_t aAbsC6F14[kNbins]={//New values from A.DiMauro 28.10.03 total 30
27     32701.4219, 17996.1141, 10039.7281, 1799.1230, 1799.1231, 1799.1231, 1241.4091, 179.0987, 179.0986, 179.0987,
28       179.0987,   118.9800,    39.5058,   23.7244,   11.1283,    7.1573,    3.6249,   2.1236,   0.7362,   0.5348,
29         0.3387,     0.3074,     0.3050,    0.0001,    0.0001,    0.0001,    0.0001,   0.0001,   0.0001,   0.0001};    
30   Float_t aAbsSiO2[kNbins]={//New values from A.DiMauro 28.10.03 total 30
31        34.4338, 30.5424, 30.2584, 31.4928, 31.7868, 17.8397, 9.3410, 6.4492, 6.1128, 5.8128,
32         5.5589,  5.2877,  5.0162,  4.7999,  4.5734,  4.2135, 3.7471, 2.6033, 1.5223, 0.9658,
33         0.4242,  0.2500,  0.1426,  0.0863,  0.0793,  0.0724, 0.0655, 0.0587, 0.0001, 0.0001};
34     
35   Float_t aAbsOpSiO2[kNbins];
36   Float_t aAbsCH4[kNbins];
37   Float_t aAbsGrid[kNbins];
38   Float_t aAbsCsI[kNbins];
39   Float_t aAbsGel[kNbins];
40   Float_t aAbsRef[kNbins];//abs prob on reflector
41   
42     for (i=0;i<kNbins;i++){
43       aAbsCH4[i]         =AliRICHParam::AbsCH4(AliRICHParam::PhotonEnergy(i)); 
44       aAbsOpSiO2[i]      =1e-5; 
45       aAbsCsI[i]         =1e-4; 
46       aAbsGrid[i]        =1e-4;       
47       aAbsGel[i]         =AliRICHParam::AbsGel(AliRICHParam::PhotonEnergy(i)); 
48       aAbsRef[i]         =0.01;//1% reflector abs prob 
49     }
50     
51   Float_t aQeCsI[kNbins] = {//New values from A.DiMauro 28.10.03 total 31
52                             0.0002, 0.0006, 0.0007, 0.0010, 0.0049, 0.0073, 0.0104, 0.0519, 0.0936, 0.1299,
53                             0.1560, 0.1768, 0.1872, 0.1976, 0.2142, 0.2288, 0.2434, 0.2599, 0.2673, 0.2808,
54                             0.2859, 0.2954, 0.3016, 0.3120, 0.3172, 0.3224, 0.3266, 0.3328, 0.3359, 0.3390};
55 //                             0.3431};
56
57   
58   Float_t aQeAll[kNbins];
59   for(i=0;i<kNbins;i++){
60     aQeCsI[i]/= (1.0-Fresnel(AliRICHParam::PhotonEnergy(i),1.0,0)); //FRESNEL LOSS CORRECTION
61     aQeAll[i]=1; //QE for all other materials except for PC must be 1.
62   }
63        
64                     
65 #ifdef __CINT__
66
67   Float_t aAbsC6F14old[kNbins]={//previous values 26 in total added 0.0001 to the 30
68       179.0987, 179.0987, 179.0987, 179.0987, 179.0987, 179.0987, 179.0987, 179.0987, 179.0987, 142.9206, 
69        56.6496,  25.5862,  13.9529,  12.0391,  10.4295,   8.8042,   7.0690,   4.4611,   2.0284,   1.2930, 
70         0.5772,   0.4074,   0.3350,   0.0001,   0.0001,   0.0001,   0.0001,   0.0001,   0.0001,   0.0001};
71   Float_t aAbsSiO2old[kNbins]={//previous values 26 in total added 0.0001 to the 30
72       105.8000,  65.5200,  48.5800,  42.8500,  35.7900,  31.2620,  28.5980,  27.5270,  25.0070,  22.8150, 
73        21.0040,  19.2660,  17.5250,  15.8780,  14.1770,  11.7190,   9.2820,   6.6200,   4.0930,   2.6010, 
74         1.1490,   0.6670,   0.3630,   0.1920,   0.1500,   0.1090,   0.0001,   0.0001,   0.0001,   0.0001};
75   Float_t aQeCsIold[kNbins]={
76     0.0002, 0.0006, 0.0007, 0.0050, 0.0075, 0.0101, 0.0243, 0.0405, 0.0689, 0.1053, 
77     0.1215, 0.1417, 0.1579, 0.1620, 0.1661, 0.1677, 0.1743, 0.1768, 0.1793, 0.1826,
78     0.1859, 0.1876, 0.1892, 0.1909, 0.2075, 0.2158, 0.0001, 0.0001, 0.0001, 0.0001 };      
79
80 //Now plot all the thigs  
81 //Markers and colors for different curves:
82     
83   const Double_t kWidth=0.25,kHeight=0.2;  
84   const Int_t kC6F14Marker=24 , kC6F14Color=kRed;  
85   const Int_t kCH4Marker  =25 , kCH4Color  =kGreen;  
86   const Int_t kSiO2M      =26 , kSiO2Color =kBlue;  
87   const Int_t kCsIMarker  = 2 , kCsIColor  =kMagenta;  
88   const Int_t kGelMarker  =27 , kGelColor  =46;  
89   const Int_t kRefMarker  =28 , kRefColor  =47;  
90   
91   TCanvas *pC=new TCanvas("c1","RICH optics to check",1100,900);
92
93   pC->Divide(3,2);           
94 //Ref index           
95   pC->cd(1);                 
96   TGraph *pIdxC6F14=new TGraph(kNbins,aPckov,aIdxC6F14);pIdxC6F14->SetMarkerStyle(kC6F14Marker);pIdxC6F14->SetMarkerColor(kC6F14Color);
97   TGraph *pIdxSiO2 =new TGraph(kNbins,aPckov,aIdxSiO2); pIdxSiO2 ->SetMarkerStyle(kSiO2M)      ;pIdxSiO2 ->SetMarkerColor(kSiO2Color);  
98   TGraph *pIdxCH4  =new TGraph(kNbins,aPckov,aIdxCH4);  pIdxCH4  ->SetMarkerStyle(kCH4Marker)       ;pIdxCH4  ->SetMarkerColor(kCH4Color);
99   TGraph *pIdxGel  =new TGraph(kNbins,aPckov,aIdxGel);  pIdxGel  ->SetMarkerStyle(kGelMarker)  ;pIdxGel  ->SetMarkerColor(kGelColor);
100   
101   TMultiGraph *pIdxMG=new TMultiGraph();  TLegend *pIdxLe=new TLegend(0.5,0.21,0.5+kWidth,0.21+kHeight);
102   pIdxMG->Add(pIdxC6F14); pIdxLe->AddEntry(pIdxC6F14,"C6F14"  ,"p");            
103   pIdxMG->Add(pIdxSiO2) ; pIdxLe->AddEntry(pIdxSiO2 ,"SiO2"   ,"p");          
104   pIdxMG->Add(pIdxCH4)  ; pIdxLe->AddEntry(pIdxCH4  ,"CH4"    ,"p");          
105   pIdxMG->Add(pIdxGel)  ; pIdxLe->AddEntry(pIdxGel  ,"Aerogel","p");          
106   pIdxMG->Draw("APL");
107   pIdxMG->SetTitle("Refractive index");  pIdxMG->GetXaxis()->SetTitle("energy, GeV");
108   pIdxMG->Draw("APL");
109   pIdxLe->Draw();      
110 //Absorbtion
111   pC->cd(2);
112   gPad->SetLogy();
113   TGraph *pAbsC6F14=new TGraph(kNbins,aPckov,aAbsC6F14);pAbsC6F14->SetMarkerStyle(kC6F14Marker); pAbsC6F14->SetMarkerColor(kC6F14Color);
114   TGraph *pAbsSiO2 =new TGraph(kNbins,aPckov,aAbsSiO2) ;pAbsSiO2 ->SetMarkerStyle(kSiO2M)      ; pAbsSiO2 ->SetMarkerColor(kSiO2Color);
115   TGraph *pAbsCH4  =new TGraph(kNbins,aPckov,aAbsCH4)  ;pAbsCH4  ->SetMarkerStyle(kCH4Marker)  ; pAbsCH4  ->SetMarkerColor(kCH4Color);
116   TGraph *pAbsGel  =new TGraph(kNbins,aPckov,aAbsGel)  ;pAbsGel  ->SetMarkerStyle(kGelMarker)  ; pAbsGel  ->SetMarkerColor(kGelColor);
117   TGraph *pAbsRef  =new TGraph(kNbins,aPckov,aAbsRef)  ;pAbsRef  ->SetMarkerStyle(kRefMarker)  ; pAbsRef  ->SetMarkerColor(kRefColor);
118   
119   TMultiGraph *pAbsMG=new TMultiGraph();  TLegend *pAbsLe=new TLegend(0.2,0.15,0.2+kWidth,0.15+kHeight);
120   pAbsMG->Add(pAbsC6F14);      pAbsLe->AddEntry(pAbsC6F14,  "C6F14"    ,"p"); 
121   pAbsMG->Add(pAbsSiO2);       pAbsLe->AddEntry(pAbsSiO2 ,  "SiO2"     ,"p"); 
122   pAbsMG->Add(pAbsCH4);        pAbsLe->AddEntry(pAbsCH4  ,  "CH4"      ,"p"); 
123   pAbsMG->Add(pAbsGel);        pAbsLe->AddEntry(pAbsGel  ,  "Aerogel"  ,"p"); 
124   pAbsMG->Add(pAbsRef);        pAbsLe->AddEntry(pAbsRef  ,  "Reflector","p"); 
125   pAbsMG->Draw("APL"); 
126   pAbsMG->SetTitle("Absorbtion length,cm"); pAbsMG->GetXaxis()->SetTitle("energy, GeV");  
127   pAbsMG->Draw("APL");
128   pAbsLe->Draw();   
129 //QE  
130   pC->cd(4);
131   TGraph *pQeCsI=new TGraph(kNbins,aPckov,aQeCsI); pQeCsI->SetMarkerStyle(kCsIMarker); pQeCsI->SetMarkerColor(kCsIColor);
132   pQeCsI->Draw("APL");
133   pQeCsI->SetTitle("QE");  pQeCsI->GetXaxis()->SetTitle("energy, GeV");
134   pQeCsI->Draw("APL");  
135 //transmission  
136   pC->cd(5);
137   Float_t mm =0.1;Float_t cm=1.; 
138   Float_t aTrC6F14[kNbins],aTrSiO2[kNbins],aTrCH4[kNbins];
139   Float_t aTotTr[kNbins];
140   for(Int_t i=0;i<kNbins;i++){//calculate probability for photon to survive during transversing a volume of material with absorption length  
141     aTrC6F14[i]=TMath::Exp(-15*mm/(aAbsC6F14[i]+0.0001));
142     aTrSiO2[i] =TMath::Exp(-5*mm/(aAbsSiO2[i] +0.0001));
143     aTrCH4[i]  =TMath::Exp(-8*cm/(aAbsCH4[i]  +0.0001));    
144     aTotTr[i]    =aTrC6F14[i]*aTrSiO2[i]*aTrCH4[i]*aQeCsI[i];
145   }
146   TGraph *pTrC6F14=new TGraph(kNbins,aPckov,aTrC6F14);pTrC6F14->SetMarkerStyle(kC6F14Marker);pTrC6F14->SetMarkerColor(kC6F14Color);  
147   TGraph *pTrSiO2 =new TGraph(kNbins,aPckov,aTrSiO2) ;pTrSiO2 ->SetMarkerStyle(kSiO2M)      ;pTrSiO2 ->SetMarkerColor(kSiO2Color);  
148   TGraph *pTrCH4  =new TGraph(kNbins,aPckov,aTrCH4)  ;pTrCH4  ->SetMarkerStyle(kCH4Marker)  ;pTrCH4  ->SetMarkerColor(kCH4Color);   
149   TGraph *pTotTr  =new TGraph(kNbins,aPckov,aTotTr)  ;pTotTr  ->SetMarkerStyle(30)          ;pTotTr  ->SetMarkerColor(kYellow);
150   
151   TMultiGraph *pTrMG=new TMultiGraph();  TLegend *pTrLe=new TLegend(0.2,0.4,0.2+kWidth,0.4+kHeight);
152   pTrMG->Add(pQeCsI);       pTrLe->AddEntry(pQeCsI,   "CsI QE", "p");  
153   pTrMG->Add(pTrC6F14);     pTrLe->AddEntry(pTrC6F14, "C6F14", "p");  
154   pTrMG->Add(pTrSiO2);      pTrLe->AddEntry(pTrSiO2,  "SiO2",  "p");          
155   pTrMG->Add(pTrCH4);       pTrLe->AddEntry(pTrCH4,   "CH4",   "p");          
156   pTrMG->Add(pTotTr);       pTrLe->AddEntry(pTotTr,   "total", "p");          
157   pTrMG->Draw("APL");
158   pTrMG->SetTitle("Transmission");  pTrMG->GetXaxis()->SetTitle("energy, GeV");
159   pTrMG->Draw("APL");
160   pTrLe->Draw();      
161 //comparison of new and old staff
162   pC->cd(6);
163   TGraph *pQeCsIold=new TGraph(kNbins,aPckov,aQeCsIold); pQeCsIold->SetMarkerStyle(kC6F14Marker);pQeCsIold->SetMarkerColor(kC6F14Color);  
164   
165   TMultiGraph *pCompMG=new TMultiGraph;  TLegend *pCompLe=new TLegend(0.2,0.6,0.2+kWidth,0.6+kHeight);
166   pCompMG->Add(pQeCsI);       pCompLe->AddEntry(pQeCsI,    "QE new 30.10.03", "p");  
167   pCompMG->Add(pQeCsIold);    pCompLe->AddEntry(pQeCsIold, "QE old 01.01.02", "p");
168   pCompMG->Draw("APL");
169   pCompMG->SetTitle("Comparison of new and old staff");
170   pCompMG->GetXaxis()->SetTitle("energy, GeV");
171   pCompMG->Draw("APL");
172   pCompLe->Draw();      
173     
174   return;
175   TCanvas *pQeC=new TCanvas("pQeC","CsI QE currently all the same",800,900); 
176   pQeC->Divide(2,4);
177   for(int i=1;i<=7;i++){
178     pQeC->cd(i);  
179     switch(i){
180       case 1: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI);pQeCsIGr->SetTitle("Module 1");break;
181       case 2: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI);pQeCsIGr->SetTitle("Module 2");break;
182       case 3: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI);pQeCsIGr->SetTitle("Module 3");break;
183       case 4: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI);pQeCsIGr->SetTitle("Module 4");break;
184       case 5: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI);pQeCsIGr->SetTitle("Module 5");break;
185       case 6: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI);pQeCsIGr->SetTitle("Module 6");break;
186       case 7: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI);pQeCsIGr->SetTitle("Module 7");break;
187     }
188     pQeCsIGr->SetMarkerStyle(kCsIMarker); pQeCsIGr->SetMarkerColor(kCsIColor);  
189     pQeCsIGr->Draw("APL");
190     pQeCsIGr->GetXaxis()->SetTitle("energy, GeV");
191     pQeCsIGr->Draw("APL");
192   }
193 }//main
194 //__________________________________________________________________________________________________
195 Float_t Fresnel(Float_t ene,Float_t pdoti, Bool_t pola)
196 {//ENE(EV), PDOTI=COS(INC.ANG.), PDOTR=COS(POL.PLANE ROT.ANG.)
197     
198     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,
199                       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,
200                       7.8,7.9,8.0,8.1,8.2,8.3,8.4,8.5};
201      
202
203     Float_t csin[36] = {2.14,2.21,2.33,2.48,2.76,2.97,2.99,2.59,2.81,3.05,
204                         2.86,2.53,2.55,2.66,2.79,2.96,3.18,3.05,2.84,2.81,2.38,2.11,
205                         2.01,2.13,2.39,2.73,3.08,3.15,2.95,2.73,2.56,2.41,2.12,1.95,
206                         1.72,1.53};
207       
208     Float_t csik[36] = {0.,0.,0.,0.,0.,0.196,0.408,0.208,0.118,0.49,0.784,0.543,
209                         0.424,0.404,0.371,0.514,0.922,1.102,1.139,1.376,1.461,1.253,0.878,
210                         0.69,0.612,0.649,0.824,1.347,1.571,1.678,1.763,1.857,1.824,1.824,
211                         1.714,1.498};
212     Float_t xe=ene;
213     Int_t  j=Int_t(xe*10)-49;
214     Float_t cn=csin[j]+((csin[j+1]-csin[j])/0.1)*(xe-en[j]);
215     Float_t ck=csik[j]+((csik[j+1]-csik[j])/0.1)*(xe-en[j]);
216
217     //FORMULAE FROM HANDBOOK OF OPTICS, 33.23 OR
218     //W.R. HUNTER, J.O.S.A. 54 (1964),15 , J.O.S.A. 55(1965),1197
219
220     Float_t sinin=TMath::Sqrt(1-pdoti*pdoti);
221     Float_t tanin=sinin/pdoti;
222
223     Float_t c1=cn*cn-ck*ck-sinin*sinin;
224     Float_t c2=4*cn*cn*ck*ck;
225     Float_t aO=TMath::Sqrt(0.5*(TMath::Sqrt(c1*c1+c2)+c1));
226     Float_t b2=0.5*(TMath::Sqrt(c1*c1+c2)-c1);
227     
228     Float_t rs=((aO-pdoti)*(aO-pdoti)+b2)/((aO+pdoti)*(aO+pdoti)+b2);
229     Float_t rp=rs*((aO-sinin*tanin)*(aO-sinin*tanin)+b2)/((aO+sinin*tanin)*(aO+sinin*tanin)+b2);
230     
231
232     //CORRECTION FACTOR FOR SURFACE ROUGHNESS
233     //B.J. STAGG  APPLIED OPTICS, 30(1991),4113
234
235     Float_t sigraf=18.;
236     Float_t lamb=1240/ene;
237     Float_t fresn;
238  
239     Float_t  rO=TMath::Exp(-(4*TMath::Pi()*pdoti*sigraf/lamb)*(4*TMath::Pi()*pdoti*sigraf/lamb));
240
241     if(pola)
242     {
243         Float_t pdotr=0.8;                                 //DEGREE OF POLARIZATION : 1->P , -1->S
244         fresn=0.5*(rp*(1+pdotr)+rs*(1-pdotr));
245     }
246     else
247         fresn=0.5*(rp+rs);
248       
249     fresn = fresn*rO;
250     return(fresn);
251 }//Fresnel()
252 #endif