RICH/Opticals.h

   1 #ifdef __CINT__
   2 void Opticals()
   3 {
   4   gROOT->Reset();
   5 #endif
   6   int i;
   7   const Int_t kNbins=26;
   8   Float_t aPckov[kNbins];
   9   for(i=0;i<kNbins;i++){    //Photons energy intervals
  10     aPckov[i]=(Float_t(i)*0.1+5.5)*1e-9;
  11   }
  12
  13   Float_t aIndexFreon[kNbins];
  14   Float_t aIndexQuartz[kNbins];
  15   Float_t aIndexOpaqueQuartz[kNbins];
  16   Float_t aIndexCH4[kNbins];
  17   Float_t aIndexGrid[kNbins];
  18
  19   Float_t  e1= 10.666;Float_t  e2= 18.125;  Float_t  f1= 46.411; Float_t  f2= 228.71;
  20   for (i=0;i<kNbins;i++){
  21     aIndexFreon[i]        = aPckov[i] * .0172 * 1e9 + 1.177;
  22     Float_t ene=aPckov[i]*1e9;
  23     Float_t a=f1/(e1*e1 - ene*ene);
  24     Float_t b=f2/(e2*e2 - ene*ene);
  25     aIndexQuartz[i]        = TMath::Sqrt(1. + a + b );
  26     aIndexOpaqueQuartz[i]  =1;
  27     aIndexCH4[i]      =1.000444;
  28     aIndexGrid[i]         =1;
  29   }
  30
  31   Float_t aAbsFreon[kNbins]={179.0987,  179.0987,    179.0987,  179.0987,   179.0987,
  32                              179.0987,  179.0987,    179.0987,  179.0987,   142.9206,
  33                               56.64957,  25.58622,    13.95293,  12.03905,   10.42953,
  34                                8.804196,  7.069031,    4.461292,  2.028366,   1.293013,
  35                                0.577267,  0.40746,     0.334964,  0.0,        0.0,
  36                                0};
  37
  38
  39   Float_t aAbsQuartz[kNbins]={105.8,   65.52,  48.58,  42.85,   35.79,
  40                                31.262, 28.598, 27.527, 25.007,  22.815,
  41                                21.004, 19.266, 17.525, 15.878,  14.177,
  42                                11.719, 9.282,   6.62,   4.0925,  2.601,
  43                                1.149,  0.667,   0.3627, 0.192,   0.1497,
  44                                0.10857};
  45
  46   Float_t aAbsOpaqueQuartz[kNbins];
  47   Float_t aAbsCH4[kNbins];
  48   Float_t aAbsGrid[kNbins];
  49   Float_t aAbsCsI[kNbins];
  50     for (i=0;i<kNbins;i++){
  51       aAbsCH4[i]    =AbsoCH4(aPckov[i]*1e9);
  52       aAbsOpaqueQuartz[i]=1e-5;
  53       aAbsCsI[i]        =1e-4;
  54       aAbsGrid[i]       =1e-4;
  55     }
  56
  57   Float_t aQeCsI1[kNbins] = {0.000199999995, 0.000600000028, 0.000699999975, 0.00499999989, 0.00749999983,
  58                             0.010125,       0.0242999997,   0.0405000001,   0.0688500032,  0.105299994,
  59                             0.121500008,    0.141749993,    0.157949999,    0.162,         0.166050002,
  60                             0.167669997,    0.174299985,    0.176789999,    0.179279998,   0.182599992,
  61                             0.18592,        0.187579989,    0.189239994,    0.190899998,   0.207499996,
  62                             0.215799987};
  63   Float_t aQeCsI2[kNbins] = {0.000199999995, 0.000600000028, 0.000699999975, 0.00499999989, 0.00749999983,
  64                             0.010125,       0.0242999997,   0.0405000001,   0.0688500032,  0.105299994,
  65                             0.121500008,    0.141749993,    0.157949999,    0.162,         0.166050002,
  66                             0.167669997,    0.174299985,    0.176789999,    0.179279998,   0.182599992,
  67                             0.18592,        0.187579989,    0.189239994,    0.190899998,   0.207499996,
  68                             0.215799987};
  69   Float_t aQeCsI3[kNbins] = {0.000199999995, 0.000600000028, 0.000699999975, 0.00499999989, 0.00749999983,
  70                             0.010125,       0.0242999997,   0.0405000001,   0.0688500032,  0.105299994,
  71                             0.121500008,    0.141749993,    0.157949999,    0.162,         0.166050002,
  72                             0.167669997,    0.174299985,    0.176789999,    0.179279998,   0.182599992,
  73                             0.18592,        0.187579989,    0.189239994,    0.190899998,   0.207499996,
  74                             0.215799987};
  75   Float_t aQeCsI4[kNbins] = {0.000199999995, 0.000600000028, 0.000699999975, 0.00499999989, 0.00749999983,
  76                             0.010125,       0.0242999997,   0.0405000001,   0.0688500032,  0.105299994,
  77                             0.121500008,    0.141749993,    0.157949999,    0.162,         0.166050002,
  78                             0.167669997,    0.174299985,    0.176789999,    0.179279998,   0.182599992,
  79                             0.18592,        0.187579989,    0.189239994,    0.190899998,   0.207499996,
  80                             0.215799987};
  81   Float_t aQeCsI5[kNbins] = {0.000199999995, 0.000600000028, 0.000699999975, 0.00499999989, 0.00749999983,
  82                             0.010125,       0.0242999997,   0.0405000001,   0.0688500032,  0.105299994,
  83                             0.121500008,    0.141749993,    0.157949999,    0.162,         0.166050002,
  84                             0.167669997,    0.174299985,    0.176789999,    0.179279998,   0.182599992,
  85                             0.18592,        0.187579989,    0.189239994,    0.190899998,   0.207499996,
  86                             0.215799987};
  87   Float_t aQeCsI6[kNbins] = {0.000199999995, 0.000600000028, 0.000699999975, 0.00499999989, 0.00749999983,
  88                             0.010125,       0.0242999997,   0.0405000001,   0.0688500032,  0.105299994,
  89                             0.121500008,    0.141749993,    0.157949999,    0.162,         0.166050002,
  90                             0.167669997,    0.174299985,    0.176789999,    0.179279998,   0.182599992,
  91                             0.18592,        0.187579989,    0.189239994,    0.190899998,   0.207499996,
  92                             0.215799987};
  93   Float_t aQeCsI7[kNbins] = {0.000199999995, 0.000600000028, 0.000699999975, 0.00499999989, 0.00749999983,
  94                             0.010125,       0.0242999997,   0.0405000001,   0.0688500032,  0.105299994,
  95                             0.121500008,    0.141749993,    0.157949999,    0.162,         0.166050002,
  96                             0.167669997,    0.174299985,    0.176789999,    0.179279998,   0.182599992,
  97                             0.18592,        0.187579989,    0.189239994,    0.190899998,   0.207499996,
  98                             0.215799987};
  99   Float_t aQeAll[kNbins];
 100   for(i=0;i<kNbins;i++){
 101     aQeCsI1[i]/= (1.0-Fresnel(aPckov[i]*1e9,1.0,0)); //FRESNEL LOSS CORRECTION
 102     aQeCsI2[i]/= (1.0-Fresnel(aPckov[i]*1e9,1.0,0)); //FRESNEL LOSS CORRECTION
 103     aQeCsI3[i]/= (1.0-Fresnel(aPckov[i]*1e9,1.0,0)); //FRESNEL LOSS CORRECTION
 104     aQeCsI4[i]/= (1.0-Fresnel(aPckov[i]*1e9,1.0,0)); //FRESNEL LOSS CORRECTION
 105     aQeCsI5[i]/= (1.0-Fresnel(aPckov[i]*1e9,1.0,0)); //FRESNEL LOSS CORRECTION
 106     aQeCsI6[i]/= (1.0-Fresnel(aPckov[i]*1e9,1.0,0)); //FRESNEL LOSS CORRECTION
 107     aQeCsI7[i]/= (1.0-Fresnel(aPckov[i]*1e9,1.0,0)); //FRESNEL LOSS CORRECTION
 108     aQeAll[i]=1; //QE for all other materials except for PC must be 1.
 109   }
 110
 111
 112 #ifdef __CINT__
 113
 114
 115 //Freon, Quartz, Opaque Quartz,Methane,CsI,Grid
 116   const Int_t kFreonMarker= 24; const Int_t kFreonColor=kRed;
 117   const Int_t kCH4Marker= 25;   const Int_t kCH4Color=kGreen;
 118   const Int_t kSiO2Marker= 26;  const Int_t kSiO2Color=kBlue;
 119   const Int_t kCsIMarker = 2;   const Int_t kCsIColor =kMagenta;
 120
 121   TCanvas *pC=new TCanvas("c1","RICH optics to check",800,900);
 122
 123   pC->Divide(2,2);
 124
 125   pC->cd(1);
 126   TGraph *pAbsFreonGr=new TGraph(kNbins,aPckov,aAbsFreon);
 127   pAbsFreonGr->SetMarkerStyle(kFreonMarker); pAbsFreonGr->SetMarkerColor(kFreonColor);
 128   TGraph *pAbsSiO2Gr=new TGraph(kNbins,aPckov,aAbsQuartz);
 129   pAbsSiO2Gr->SetMarkerStyle(kSiO2Marker); pAbsSiO2Gr->SetMarkerColor(kSiO2Color);
 130   TMultiGraph *pAbsMG=new TMultiGraph();
 131   TLegend *pAbsLegend=new TLegend(0.6,0.3,0.85,0.5);
 132   pAbsMG->Add(pAbsFreonGr);     pAbsLegend->AddEntry(pAbsFreonGr,   "freon","p");            //1
 133   pAbsMG->Add(pAbsSiO2Gr);      pAbsLegend->AddEntry(pAbsSiO2Gr,  "quartz","p");          //2
 134   pAbsMG->Draw("APL");
 135   pAbsMG->GetXaxis()->SetTitle("energy, GeV");
 136   pAbsMG->GetYaxis()->SetTitle("absorption length, cm");
 137   pAbsMG->Draw("APL");
 138   pAbsLegend->Draw();
 139
 140
 141   pC->cd(2);
 142   TGraph *pIndexFreonGr=new TGraph(kNbins,aPckov,aIndexFreon);
 143   pIndexFreonGr->SetMarkerStyle(kFreonMarker); pIndexFreonGr->SetMarkerColor(kFreonColor);
 144   TGraph *pIndexSiO2Gr=new TGraph(kNbins,aPckov,aIndexQuartz);
 145   pIndexSiO2Gr->SetMarkerStyle(kSiO2Marker); pIndexSiO2Gr->SetMarkerColor(kSiO2Color);
 146   TGraph *pIndexCH4Gr=new TGraph(kNbins,aPckov,aIndexCH4);
 147   pIndexCH4Gr->SetMarkerStyle(kCH4Marker); pIndexCH4Gr->SetMarkerColor(kCH4Color);
 148   TMultiGraph *pIndexMG=new TMultiGraph();
 149   TLegend *pIndexLegend=new TLegend(0.6,0.2,0.85,0.4);
 150   pIndexMG->Add(pIndexFreonGr);     pIndexLegend->AddEntry(pIndexFreonGr,   "freon","p");            //1
 151   pIndexMG->Add(pIndexSiO2Gr);      pIndexLegend->AddEntry(pIndexSiO2Gr,  "quartz","p");
 152   pIndexMG->Add(pIndexCH4Gr);       pIndexLegend->AddEntry(pIndexCH4Gr,   "CH4","p");
 153   pIndexMG->Draw("APL");
 154   pIndexMG->GetXaxis()->SetTitle("energy, GeV");
 155   pIndexMG->GetYaxis()->SetTitle("refraction index");
 156   pIndexMG->Draw("APL");
 157   pIndexLegend->Draw();
 158
 159   pC->cd(3);
 160   TGraph *pAbsCH4Gr=new TGraph(kNbins,aPckov,aAbsCH4);
 161   pAbsCH4Gr->SetMarkerStyle(kCH4Marker); pAbsCH4Gr->SetMarkerColor(kCH4Color);
 162   pAbsCH4Gr->Draw("APL");
 163   pAbsCH4Gr->GetXaxis()->SetTitle("energy, GeV");
 164   pAbsCH4Gr->SetTitle("CH4 absorption length, cm");
 165   pAbsCH4Gr->Draw("APL");
 166
 167   TCanvas *pQeC=new TCanvas("pQeC","CsI QE currently all the same",800,900);
 168   pQeC->Divide(2,4);
 169   for(int i=1;i<=7;i++){
 170     pQeC->cd(i);
 171     switch(i){
 172       case 1: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI1);pQeCsIGr->SetTitle("Module 1");break;
 173       case 2: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI2);pQeCsIGr->SetTitle("Module 2");break;
 174       case 3: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI3);pQeCsIGr->SetTitle("Module 3");break;
 175       case 4: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI4);pQeCsIGr->SetTitle("Module 4");break;
 176       case 5: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI5);pQeCsIGr->SetTitle("Module 5");break;
 177       case 6: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI6);pQeCsIGr->SetTitle("Module 6");break;
 178       case 7: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI7);pQeCsIGr->SetTitle("Module 7");break;
 179     }
 180     pQeCsIGr->SetMarkerStyle(kCsIMarker); pQeCsIGr->SetMarkerColor(kCsIColor);
 181     pQeCsIGr->Draw("APL");
 182     pQeCsIGr->GetXaxis()->SetTitle("energy, GeV");
 183     pQeCsIGr->Draw("APL");
 184   }
 185 }//main
 186
 187 Float_t AbsoCH4(Float_t x)
 188 {
 189
 190     //KLOSCH,SCH4(9),WL(9),EM(9),ALENGTH(31)
 191     Float_t sch4[9] = {.12,.16,.23,.38,.86,2.8,7.9,28.,80.};              //MB X 10^22
 192     //Float_t wl[9] = {153.,152.,151.,150.,149.,148.,147.,146.,145};
 193     Float_t em[9] = {8.1,8.158,8.212,8.267,8.322,8.378,8.435,8.493,8.55};
 194     const Float_t kLosch=2.686763E19;                                      // LOSCHMIDT NUMBER IN CM-3
 195     const Float_t kIgas1=100, kIgas2=0, kOxy=10., kWater=5., kPressure=750.,kTemperature=283.;
 196     Float_t pn=kPressure/760.;
 197     Float_t tn=kTemperature/273.16;
 198
 199
 200 // ------- METHANE CROSS SECTION -----------------
 201 // ASTROPH. J. 214, L47 (1978)
 202
 203     Float_t sm=0;
 204     if (x<7.75)
 205         sm=.06e-22;
 206
 207     if(x>=7.75 && x<=8.1)
 208     {
 209         Float_t c0=-1.655279e-1;
 210         Float_t c1=6.307392e-2;
 211         Float_t c2=-8.011441e-3;
 212         Float_t c3=3.392126e-4;
 213         sm=(c0+c1*x+c2*x*x+c3*x*x*x)*1.e-18;
 214     }
 215
 216     if (x> 8.1)
 217     {
 218         Int_t j=0;
 219         while (x<=em[j] && x>=em[j+1])
 220         {
 221             j++;
 222             Float_t a=(sch4[j+1]-sch4[j])/(em[j+1]-em[j]);
 223             sm=(sch4[j]+a*(x-em[j]))*1e-22;
 224         }
 225     }
 226
 227     Float_t dm=(kIgas1/100.)*(1.-((kOxy+kWater)/1.e6))*kLosch*pn/tn;
 228     Float_t abslm=1./sm/dm;
 229
 230 //    ------- ISOBUTHANE CROSS SECTION --------------
 231 //     i-C4H10 (ai) abs. length from curves in
 232 //     Lu-McDonald paper for BARI RICH workshop .
 233 //     -----------------------------------------------------------
 234
 235     Float_t ai;
 236     Float_t absli;
 237     if (kIgas2 != 0)
 238     {
 239         if (x<7.25)
 240             ai=100000000.;
 241
 242         if(x>=7.25 && x<7.375)
 243             ai=24.3;
 244
 245         if(x>=7.375)
 246             ai=.0000000001;
 247
 248         Float_t si = 1./(ai*kLosch*273.16/293.);                    // ISOB. CRO.SEC.IN CM2
 249         Float_t di=(kIgas2/100.)*(1.-((kOxy+kWater)/1.e6))*kLosch*pn/tn;
 250         absli =1./si/di;
 251     }
 252     else
 253         absli=1.e18;
 254 //    ---------------------------------------------------------
 255 //
 256 //       transmission of O2
 257 //
 258 //       y= path in cm, x=energy in eV
 259 //       so= cross section for UV absorption in cm2
 260 //       do= O2 molecular density in cm-3
 261 //    ---------------------------------------------------------
 262
 263     Float_t abslo;
 264     Float_t so=0;
 265     if(x>=6.0)
 266     {
 267         if(x>=6.0 && x<6.5)
 268         {
 269             so=3.392709e-13 * TMath::Exp(2.864104 *x);
 270             so=so*1e-18;
 271         }
 272
 273         if(x>=6.5 && x<7.0)
 274         {
 275             so=2.910039e-34 * TMath::Exp(10.3337*x);
 276             so=so*1e-18;
 277         }
 278
 279
 280         if (x>=7.0)
 281         {
 282             Float_t a0=-73770.76;
 283             Float_t a1=46190.69;
 284             Float_t a2=-11475.44;
 285             Float_t a3=1412.611;
 286             Float_t a4=-86.07027;
 287             Float_t a5=2.074234;
 288             so= a0+(a1*x)+(a2*x*x)+(a3*x*x*x)+(a4*x*x*x*x)+(a5*x*x*x*x*x);
 289             so=so*1e-18;
 290         }
 291
 292         Float_t dox=(kOxy/1e6)*kLosch*pn/tn;
 293         abslo=1./so/dox;
 294     }
 295     else
 296         abslo=1.e18;
 297 //     ---------------------------------------------------------
 298 //
 299 //       transmission of H2O
 300 //
 301 //       y= path in cm, x=energy in eV
 302 //       sw= cross section for UV absorption in cm2
 303 //       dw= H2O molecular density in cm-3
 304 //     ---------------------------------------------------------
 305
 306     Float_t abslw;
 307
 308     Float_t b0=29231.65;
 309     Float_t b1=-15807.74;
 310     Float_t b2=3192.926;
 311     Float_t b3=-285.4809;
 312     Float_t b4=9.533944;
 313
 314     if(x>6.75)
 315     {
 316         Float_t sw= b0+(b1*x)+(b2*x*x)+(b3*x*x*x)+(b4*x*x*x*x);
 317         sw=sw*1e-18;
 318         Float_t dw=(kWater/1e6)*kLosch*pn/tn;
 319         abslw=1./sw/dw;
 320     }
 321     else
 322         abslw=1.e18;
 323
 324 //    ---------------------------------------------------------
 325
 326     Float_t alength=1./(1./abslm+1./absli+1./abslo+1./abslw);
 327     return (alength);
 328 }//AbsoCH4
 329
 330
 331 Float_t Fresnel(Float_t ene,Float_t pdoti, Bool_t pola)
 332 {
 333
 334     //ENE(EV), PDOTI=COS(INC.ANG.), PDOTR=COS(POL.PLANE ROT.ANG.)
 335
 336     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,
 337                       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,
 338                       7.8,7.9,8.0,8.1,8.2,8.3,8.4,8.5};
 339
 340
 341     Float_t csin[36] = {2.14,2.21,2.33,2.48,2.76,2.97,2.99,2.59,2.81,3.05,
 342                         2.86,2.53,2.55,2.66,2.79,2.96,3.18,3.05,2.84,2.81,2.38,2.11,
 343                         2.01,2.13,2.39,2.73,3.08,3.15,2.95,2.73,2.56,2.41,2.12,1.95,
 344                         1.72,1.53};
 345
 346     Float_t csik[36] = {0.,0.,0.,0.,0.,0.196,0.408,0.208,0.118,0.49,0.784,0.543,
 347                         0.424,0.404,0.371,0.514,0.922,1.102,1.139,1.376,1.461,1.253,0.878,
 348                         0.69,0.612,0.649,0.824,1.347,1.571,1.678,1.763,1.857,1.824,1.824,
 349                         1.714,1.498};
 350     Float_t xe=ene;
 351     Int_t  j=Int_t(xe*10)-49;
 352     Float_t cn=csin[j]+((csin[j+1]-csin[j])/0.1)*(xe-en[j]);
 353     Float_t ck=csik[j]+((csik[j+1]-csik[j])/0.1)*(xe-en[j]);
 354
 355     //FORMULAE FROM HANDBOOK OF OPTICS, 33.23 OR
 356     //W.R. HUNTER, J.O.S.A. 54 (1964),15 , J.O.S.A. 55(1965),1197
 357
 358     Float_t sinin=TMath::Sqrt(1-pdoti*pdoti);
 359     Float_t tanin=sinin/pdoti;
 360
 361     Float_t c1=cn*cn-ck*ck-sinin*sinin;
 362     Float_t c2=4*cn*cn*ck*ck;
 363     Float_t aO=TMath::Sqrt(0.5*(TMath::Sqrt(c1*c1+c2)+c1));
 364     Float_t b2=0.5*(TMath::Sqrt(c1*c1+c2)-c1);
 365
 366     Float_t rs=((aO-pdoti)*(aO-pdoti)+b2)/((aO+pdoti)*(aO+pdoti)+b2);
 367     Float_t rp=rs*((aO-sinin*tanin)*(aO-sinin*tanin)+b2)/((aO+sinin*tanin)*(aO+sinin*tanin)+b2);
 368
 369
 370     //CORRECTION FACTOR FOR SURFACE ROUGHNESS
 371     //B.J. STAGG  APPLIED OPTICS, 30(1991),4113
 372
 373     Float_t sigraf=18.;
 374     Float_t lamb=1240/ene;
 375     Float_t fresn;
 376
 377     Float_t  rO=TMath::Exp(-(4*TMath::Pi()*pdoti*sigraf/lamb)*(4*TMath::Pi()*pdoti*sigraf/lamb));
 378
 379     if(pola)
 380     {
 381         Float_t pdotr=0.8;                                 //DEGREE OF POLARIZATION : 1->P , -1->S
 382         fresn=0.5*(rp*(1+pdotr)+rs*(1-pdotr));
 383     }
 384     else
 385         fresn=0.5*(rp+rs);
 386
 387     fresn = fresn*rO;
 388     return(fresn);
 389 }//Fresnel(...)
 390 #endif