RICH/Opticals.h

   1 #ifdef __CINT__
   2 void Opticals()
   3 {
   4   gROOT->Reset();
   5 #endif
   6   int i;
   7   const Int_t kNbins=30;
   8
   9   Float_t aPckov[kNbins];  for(i=0;i<kNbins;i++) aPckov[i]=(0.1*i+5.5)*1e-9;//Photons energy bins 5.5 eV - 8.5 eV step 0.1 eV
  10
  11
  12   Float_t aIdxC6F14[kNbins];
  13   Float_t aIdxCH4[kNbins];
  14   Float_t aIdxGrid[kNbins];
  15   Float_t aIdxSiO2[kNbins];
  16   Float_t aIdxOpSiO2[kNbins];
  17
  18   Float_t  e1=10.666,e2=18.125,f1=46.411,f2= 228.71; //RICH TDR page 35
  19   for (i=0;i<kNbins;i++){
  20     aIdxC6F14[i] = aPckov[i] *0.0172*1e9+1.177;
  21     aIdxSiO2[i]  = TMath::Sqrt(1.+f1/(e1*e1-TMath::Power(aPckov[i]*1e9,2))+f2/(e2*e2-TMath::Power(aPckov[i]*1e9,2)));//TDR p.35
  22     aIdxCH4[i]   =1.000444;
  23     aIdxGrid[i]  =1;
  24     aIdxOpSiO2[i]  =1;
  25   }
  26
  27   Float_t aAbsC6F14[kNbins]={//New values from A.DiMauro 28.10.03 total 31
  28     32701.4219, 17996.1141, 10039.7281, 1799.1230, 1799.1231, 1799.1231, 1241.4091, 179.0987, 179.0986, 179.0987,
  29       179.0987,   118.9800,    39.5058,   23.7244,   11.1283,    7.1573,    3.6249,   2.1236,   0.7362,   0.5348,
  30         0.3387,     0.3074,     0.3050,    0.0001,    0.0001,    0.0001,    0.0001,   0.0001,   0.0001,   0.0001};
  31   Float_t aAbsSiO2[kNbins]={//New values from A.DiMauro 28.10.03 total 31
  32        34.4338, 30.5424, 30.2584, 31.4928, 31.7868, 17.8397, 9.3410, 6.4492, 6.1128, 5.8128,
  33         5.5589,  5.2877,  5.0162,  4.7999,  4.5734,  4.2135, 3.7471, 2.6033, 1.5223, 0.9658,
  34         0.4242,  0.2500,  0.1426,  0.0863,  0.0793,  0.0724, 0.0655, 0.0587, 0.0001, 0.0001};
  35
  36   Float_t aAbsOpSiO2[kNbins];
  37   Float_t aAbsCH4[kNbins];
  38   Float_t aAbsGrid[kNbins];
  39   Float_t aAbsCsI[kNbins];
  40     for (i=0;i<kNbins;i++){
  41       aAbsCH4[i]         =AbsoCH4(aPckov[i]*1e9);
  42       aAbsOpSiO2[i]      =1e-5;
  43       aAbsCsI[i]         =1e-4;
  44       aAbsGrid[i]        =1e-4;
  45     }
  46
  47   Float_t aQeCsI[kNbins] = {//New values from A.DiMauro 28.10.03 total 31
  48                             0.0002, 0.0006, 0.0007, 0.0010, 0.0049, 0.0073, 0.0104, 0.0519, 0.0936, 0.1299,
  49                             0.1560, 0.1768, 0.1872, 0.1976, 0.2142, 0.2288, 0.2434, 0.2599, 0.2673, 0.2808,
  50                             0.2859, 0.2954, 0.3016, 0.3120, 0.3172, 0.3224, 0.3266, 0.3328, 0.3359, 0.3390};
  51 //                             0.3431};
  52
  53
  54   Float_t aQeAll[kNbins];
  55   for(i=0;i<kNbins;i++){
  56     aQeCsI[i]/= (1.0-Fresnel(aPckov[i]*1e9,1.0,0)); //FRESNEL LOSS CORRECTION
  57     aQeAll[i]=1; //QE for all other materials except for PC must be 1.
  58   }
  59
  60
  61 #ifdef __CINT__
  62
  63   Float_t aAbsC6F14old[kNbins]={//previous values 26 in total added 0.0001 to the 30
  64       179.0987, 179.0987, 179.0987, 179.0987, 179.0987, 179.0987, 179.0987, 179.0987, 179.0987, 142.9206,
  65        56.6496,  25.5862,  13.9529,  12.0391,  10.4295,   8.8042,   7.0690,   4.4611,   2.0284,   1.2930,
  66         0.5772,   0.4074,   0.3350,   0.0001,   0.0001,   0.0001,   0.0001,   0.0001,   0.0001,   0.0001};
  67   Float_t aAbsSiO2old[kNbins]={//previous values 26 in total added 0.0001 to the 30
  68       105.8000,  65.5200,  48.5800,  42.8500,  35.7900,  31.2620,  28.5980,  27.5270,  25.0070,  22.8150,
  69        21.0040,  19.2660,  17.5250,  15.8780,  14.1770,  11.7190,   9.2820,   6.6200,   4.0930,   2.6010,
  70         1.1490,   0.6670,   0.3630,   0.1920,   0.1500,   0.1090,   0.0001,   0.0001,   0.0001,   0.0001};
  71   Float_t aQeCsIold[kNbins]={
  72     0.0002, 0.0006, 0.0007, 0.0050, 0.0075, 0.0101, 0.0243, 0.0405, 0.0689, 0.1053,
  73     0.1215, 0.1417, 0.1579, 0.1620, 0.1661, 0.1677, 0.1743, 0.1768, 0.1793, 0.1826,
  74     0.1859, 0.1876, 0.1892, 0.1909, 0.2075, 0.2158, 0.0001, 0.0001, 0.0001, 0.0001 };
  75
  76 //Now plot all the thigs
  77 //Freon, Quartz, Opaque ,Methane,CsI,Grid
  78   const Int_t kC6F14M= 24; const Int_t kC6F14C=kRed;
  79   const Int_t kCH4M= 25;   const Int_t kCH4Color=kGreen;
  80   const Int_t kSiO2M= 26;  const Int_t kSiO2C=kBlue;
  81   const Int_t kCsIMarker = 2;   const Int_t kCsIColor =kMagenta;
  82
  83   TCanvas *pC=new TCanvas("c1","RICH optics to check",1100,900);
  84
  85   pC->Divide(3,2);
  86
  87   pC->cd(1);
  88   TGraph *pIdxC6F14G=new TGraph(kNbins,aPckov,aIdxC6F14);pIdxC6F14G->SetMarkerStyle(kC6F14M); pIdxC6F14G->SetMarkerColor(kC6F14C);
  89   TGraph *pIdxSiO2G =new TGraph(kNbins,aPckov,aIdxSiO2); pIdxSiO2G->SetMarkerStyle(kSiO2M);   pIdxSiO2G->SetMarkerColor(kSiO2C);
  90   TGraph *pIdxCH4G  =new TGraph(kNbins,aPckov,aIdxCH4);  pIdxCH4G->SetMarkerStyle(kCH4M);     pIdxCH4G->SetMarkerColor(kCH4Color);
  91   TMultiGraph *pIdxMG=new TMultiGraph();  TLegend *pIdxLe=new TLegend(0.6,0.2,0.85,0.4);
  92   pIdxMG->Add(pIdxC6F14G);     pIdxLe->AddEntry(pIdxC6F14G,   "C6F14","p");
  93   pIdxMG->Add(pIdxSiO2G);      pIdxLe->AddEntry(pIdxSiO2G,    "SiO2","p");
  94   pIdxMG->Add(pIdxCH4G);       pIdxLe->AddEntry(pIdxCH4G,     "CH4","p");
  95   pIdxMG->Draw("APL");pIdxMG->GetXaxis()->SetTitle("ref index versus energy, GeV");pIdxMG->Draw("APL");
  96   pIdxLe->Draw();
  97
  98
  99   pC->cd(2);
 100   TGraph *pAbsC6F14G=new TGraph(kNbins,aPckov,aAbsC6F14);pAbsC6F14G->SetMarkerStyle(kC6F14M); pAbsC6F14G->SetMarkerColor(kC6F14C);
 101   TGraph *pAbsSiO2G =new TGraph(kNbins,aPckov,aAbsSiO2); pAbsSiO2G->SetMarkerStyle(kSiO2M);   pAbsSiO2G->SetMarkerColor(kSiO2C);
 102   TMultiGraph *pAbsMG=new TMultiGraph();  TLegend *pAbsLegend=new TLegend(0.6,0.3,0.85,0.5);
 103   pAbsMG->Add(pAbsC6F14G);      pAbsLegend->AddEntry(pAbsC6F14G,  "C6F14","p");
 104   pAbsMG->Add(pAbsSiO2G);       pAbsLegend->AddEntry(pAbsSiO2G,  "SiO2", "p");
 105   pAbsMG->Draw("APL"); pAbsMG->GetXaxis()->SetTitle("absorbtion length, cm versus energy, GeV");  pAbsMG->Draw("APL");
 106   pAbsLegend->Draw();
 107
 108   pC->cd(3);
 109   TGraph *pAbsCH4Gr=new TGraph(kNbins,aPckov,aAbsCH4);
 110   pAbsCH4Gr->SetMarkerStyle(kCH4M); pAbsCH4Gr->SetMarkerColor(kCH4Color);
 111   pAbsCH4Gr->Draw("APL");
 112   pAbsCH4Gr->GetXaxis()->SetTitle("energy, GeV");
 113   pAbsCH4Gr->SetTitle("CH4 absorption length, cm");
 114   pAbsCH4Gr->Draw("APL");
 115
 116   pC->cd(4);
 117   TGraph *pQeCsIG=new TGraph(kNbins,aPckov,aQeCsI);
 118   pQeCsIG->SetMarkerStyle(kCsIMarker); pQeCsIG->SetMarkerColor(kCsIColor);
 119   pQeCsIG->Draw("APL");
 120   pQeCsIG->GetXaxis()->SetTitle("energy, GeV");
 121   pQeCsIG->SetTitle("QE");
 122   pQeCsIG->Draw("APL");
 123
 124 //transmission
 125   pC->cd(5);
 126   Float_t aTrC6F14[kNbins],aTrSiO2[kNbins],aTrCH4[kNbins];
 127   Float_t aTotTr[kNbins];
 128   for(Int_t i=0;i<kNbins;i++){
 129     aTrC6F14[i]=TMath::Exp(-AliRICHParam::FreonThickness() /(aAbsC6F14[i]+0.0001));
 130     aTrSiO2[i] =TMath::Exp(-AliRICHParam::QuartzThickness()/(aAbsSiO2[i] +0.0001));
 131     aTrCH4[i]  =TMath::Exp(-AliRICHParam::GapThickness()   /(aAbsCH4[i]  +0.0001));
 132     aTotTr[i]    =aTrC6F14[i]*aTrSiO2[i]*aTrCH4[i]*aQeCsI[i];
 133   }
 134   TGraph *pTrC6F14G=new TGraph(kNbins,aPckov,aTrC6F14);pTrC6F14G->SetMarkerStyle(kC6F14M);pTrC6F14G->SetMarkerColor(kC6F14C);
 135   TGraph *pTrSiO2G=new TGraph(kNbins,aPckov,aTrSiO2);pTrSiO2G->SetMarkerStyle(kSiO2M); pTrSiO2G->SetMarkerColor(kSiO2C);
 136   TGraph *pTrCH4G=new TGraph(kNbins,aPckov,aTrCH4);pTrCH4G->SetMarkerStyle(kCH4M); pTrCH4G->SetMarkerColor(kCH4Color);
 137   TGraph *pTotTrG=new TGraph(kNbins,aPckov,aTotTr);pTotTrG->SetMarkerStyle(30);pTotTrG->SetMarkerColor(kYellow);
 138   TMultiGraph *pTrMG=new TMultiGraph();
 139   TLegend *pTrLegend=new TLegend(0.2,0.4,0.35,0.6);
 140   pTrMG->Add(pQeCsIG);       pTrLegend->AddEntry(pQeCsIG,   "CsI QE", "p");
 141   pTrMG->Add(pTrC6F14G);     pTrLegend->AddEntry(pTrC6F14G, "C6F14", "p");
 142   pTrMG->Add(pTrSiO2G);      pTrLegend->AddEntry(pTrSiO2G,  "SiO2",  "p");
 143   pTrMG->Add(pTrCH4G);       pTrLegend->AddEntry(pTrCH4G,   "CH4",   "p");
 144   pTrMG->Add(pTotTrG);       pTrLegend->AddEntry(pTotTrG,   "total", "p");
 145   pTrMG->Draw("APL");pTrMG->GetXaxis()->SetTitle("transmission versus energy, GeV");pTrMG->Draw("APL");
 146   pTrLegend->Draw();
 147
 148   pC->cd(6);
 149
 150   TMultiGraph *pCompMG=new TMultiGraph;
 151   pCompMG->Add(pQeCsIG);
 152   pCompMG->Add(new TGraph(kNbins,aPckov,aQeCsIold));
 153   pCompMG->Draw("APL");pCompMG->GetXaxis()->SetTitle("comparison of QE versus energy, GeV");pCompMG->Draw("APL");
 154
 155   return;
 156   TCanvas *pQeC=new TCanvas("pQeC","CsI QE currently all the same",800,900);
 157   pQeC->Divide(2,4);
 158   for(int i=1;i<=7;i++){
 159     pQeC->cd(i);
 160     switch(i){
 161       case 1: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI);pQeCsIGr->SetTitle("Module 1");break;
 162       case 2: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI);pQeCsIGr->SetTitle("Module 2");break;
 163       case 3: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI);pQeCsIGr->SetTitle("Module 3");break;
 164       case 4: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI);pQeCsIGr->SetTitle("Module 4");break;
 165       case 5: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI);pQeCsIGr->SetTitle("Module 5");break;
 166       case 6: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI);pQeCsIGr->SetTitle("Module 6");break;
 167       case 7: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI);pQeCsIGr->SetTitle("Module 7");break;
 168     }
 169     pQeCsIGr->SetMarkerStyle(kCsIMarker); pQeCsIGr->SetMarkerColor(kCsIColor);
 170     pQeCsIGr->Draw("APL");
 171     pQeCsIGr->GetXaxis()->SetTitle("energy, GeV");
 172     pQeCsIGr->Draw("APL");
 173   }
 174 }//main
 175
 176 //__________________________________________________________________________________________________
 177 Float_t AbsoCH4(Float_t x)
 178 {//Evaluate the absorbtion lenght of CH4
 179   Float_t sch4[9] = {.12,.16,.23,.38,.86,2.8,7.9,28.,80.};              //MB X 10^22
 180   Float_t em[9] = {8.1,8.158,8.212,8.267,8.322,8.378,8.435,8.493,8.55};
 181   const Float_t kLoschmidt=2.686763e19;                                      // LOSCHMIDT NUMBER IN CM-3
 182   const Float_t kPressure=750.,kTemperature=283.;
 183   const Float_t pn=kPressure/760.;
 184   const Float_t tn=kTemperature/273.16;
 185   const Float_t c0=-1.655279e-1;
 186   const Float_t c1=6.307392e-2;
 187   const Float_t c2=-8.011441e-3;
 188   const Float_t c3=3.392126e-4;
 189
 190   Float_t crossSection=0;
 191   if (x<7.75)
 192     crossSection=.06e-22;
 193   else if(x>=7.75 && x<=8.1){                 //------ METHANE CROSS SECTION cm-2 ASTROPH. J. 214, L47 (1978)
 194         crossSection=(c0+c1*x+c2*x*x+c3*x*x*x)*1.e-18;
 195   }else if (x> 8.1){
 196     Int_t j=0;
 197     while (x<=em[j] || x>=em[j+1]){
 198       j++;
 199       Float_t a=(sch4[j+1]-sch4[j])/(em[j+1]-em[j]);
 200       crossSection=(sch4[j]+a*(x-em[j]))*1e-22;
 201     }
 202   }//if
 203
 204     Float_t density=kLoschmidt*pn/tn; //CH4 molecular density 1/cm-3
 205     return 1./(density*crossSection);
 206 }//AbsoCH4()
 207 //__________________________________________________________________________________________________
 208 Float_t Fresnel(Float_t ene,Float_t pdoti, Bool_t pola)
 209 {//ENE(EV), PDOTI=COS(INC.ANG.), PDOTR=COS(POL.PLANE ROT.ANG.)
 210
 211     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,
 212                       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,
 213                       7.8,7.9,8.0,8.1,8.2,8.3,8.4,8.5};
 214
 215
 216     Float_t csin[36] = {2.14,2.21,2.33,2.48,2.76,2.97,2.99,2.59,2.81,3.05,
 217                         2.86,2.53,2.55,2.66,2.79,2.96,3.18,3.05,2.84,2.81,2.38,2.11,
 218                         2.01,2.13,2.39,2.73,3.08,3.15,2.95,2.73,2.56,2.41,2.12,1.95,
 219                         1.72,1.53};
 220
 221     Float_t csik[36] = {0.,0.,0.,0.,0.,0.196,0.408,0.208,0.118,0.49,0.784,0.543,
 222                         0.424,0.404,0.371,0.514,0.922,1.102,1.139,1.376,1.461,1.253,0.878,
 223                         0.69,0.612,0.649,0.824,1.347,1.571,1.678,1.763,1.857,1.824,1.824,
 224                         1.714,1.498};
 225     Float_t xe=ene;
 226     Int_t  j=Int_t(xe*10)-49;
 227     Float_t cn=csin[j]+((csin[j+1]-csin[j])/0.1)*(xe-en[j]);
 228     Float_t ck=csik[j]+((csik[j+1]-csik[j])/0.1)*(xe-en[j]);
 229
 230     //FORMULAE FROM HANDBOOK OF OPTICS, 33.23 OR
 231     //W.R. HUNTER, J.O.S.A. 54 (1964),15 , J.O.S.A. 55(1965),1197
 232
 233     Float_t sinin=TMath::Sqrt(1-pdoti*pdoti);
 234     Float_t tanin=sinin/pdoti;
 235
 236     Float_t c1=cn*cn-ck*ck-sinin*sinin;
 237     Float_t c2=4*cn*cn*ck*ck;
 238     Float_t aO=TMath::Sqrt(0.5*(TMath::Sqrt(c1*c1+c2)+c1));
 239     Float_t b2=0.5*(TMath::Sqrt(c1*c1+c2)-c1);
 240
 241     Float_t rs=((aO-pdoti)*(aO-pdoti)+b2)/((aO+pdoti)*(aO+pdoti)+b2);
 242     Float_t rp=rs*((aO-sinin*tanin)*(aO-sinin*tanin)+b2)/((aO+sinin*tanin)*(aO+sinin*tanin)+b2);
 243
 244
 245     //CORRECTION FACTOR FOR SURFACE ROUGHNESS
 246     //B.J. STAGG  APPLIED OPTICS, 30(1991),4113
 247
 248     Float_t sigraf=18.;
 249     Float_t lamb=1240/ene;
 250     Float_t fresn;
 251
 252     Float_t  rO=TMath::Exp(-(4*TMath::Pi()*pdoti*sigraf/lamb)*(4*TMath::Pi()*pdoti*sigraf/lamb));
 253
 254     if(pola)
 255     {
 256         Float_t pdotr=0.8;                                 //DEGREE OF POLARIZATION : 1->P , -1->S
 257         fresn=0.5*(rp*(1+pdotr)+rs*(1-pdotr));
 258     }
 259     else
 260         fresn=0.5*(rp+rs);
 261
 262     fresn = fresn*rO;
 263     return(fresn);
 264 }//Fresnel()
 265 #endif