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   Float_t aPckov[kNbins];
   9   for(i=0;i<kNbins;i++){    //Photons energy bins 5.5 eV - 8.5 eV step 0.1 eV
  10     aPckov[i]=(0.1*i+5.5)*1e-9;
  11   }
  12
  13   Float_t aIndexFreon[kNbins];
  14   Float_t aIndexSiO2[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;//RICH TDR page 35
  20   for (i=0;i<kNbins;i++){
  21     aIndexFreon[i]        = aPckov[i] * .0172 * 1e9 + 1.177;
  22     Float_t ene=aPckov[i]*1e9;
  23     aIndexSiO2[i]        = TMath::Sqrt(1. + f1/(e1*e1 - ene*ene) + f2/(e2*e2 - ene*ene) );
  24     aIndexOpaqueQuartz[i]  =1;
  25     aIndexCH4[i]      =1.000444;
  26     aIndexGrid[i]         =1;
  27   }
  28
  29   Float_t aAbsFreon[kNbins]={179.0987,  179.0987,    179.0987,  179.0987,   179.0987,   //Previous values from code
  30                              179.0987,  179.0987,    179.0987,  179.0987,   142.9206,
  31                               56.64957,  25.58622,    13.95293,  12.03905,   10.42953,
  32                                8.804196,  7.069031,    4.461292,  2.028366,   1.293013,
  33                                0.577267,  0.40746,     0.334964,  0.00001,        0.00001,
  34                                0.00001};
  35
  36 //   Float_t aAbsFreon[kNbins]={32701.4219, 17996.1141, 10039.7281, 1799.1231, 1799.1231,//New values from A.DiMauro 28.10.03
  37 //                               1799.1231,  1241.4091,   179.0987,  179.0987,  179.0987,  179.0987,   179.0987,
  38 //                              179.0987,  179.0987,    179.0987,  179.0987,   142.9206,
  39 //                               56.64957,  25.58622,    13.95293,  12.03905,   10.42953,
  40 //                                8.804196,  7.069031,    4.461292,  2.028366,   1.293013,
  41 //                                0.577267,  0.40746,     0.334964,  0.0,        0.0,
  42 //                                0.0000, 0.0000, 0.0000, 0.0000, 0.0000};
  43
  44   Float_t aAbsSiO2[kNbins]={105.8,   65.52,  48.58,  42.85,   35.79,
  45                                31.262, 28.598, 27.527, 25.007,  22.815,
  46                                21.004, 19.266, 17.525, 15.878,  14.177,
  47                                11.719, 9.282,   6.62,   4.0925,  2.601,
  48                                1.149,  0.667,   0.3627, 0.192,   0.1497,
  49                                0.10857};
  50
  51   Float_t aAbsOpaqueQuartz[kNbins];
  52   Float_t aAbsCH4[kNbins];
  53   Float_t aAbsGrid[kNbins];
  54   Float_t aAbsCsI[kNbins];
  55     for (i=0;i<kNbins;i++){
  56       aAbsCH4[i]    =AbsoCH4(aPckov[i]*1e9);
  57       aAbsOpaqueQuartz[i]=1e-5;
  58       aAbsCsI[i]        =1e-4;
  59       aAbsGrid[i]       =1e-4;
  60     }
  61
  62   Float_t aQeCsI1[kNbins] = {0.0002, 0.0003, 0.0005, 0.0010, 0.0020, 0.0020, 0.0099, 0.0500, 0.1800, 0.1850,
  63                              0.1899, 0.1920, 0.1979, 0.2099, 0.2179, 0.2249, 0.2349, 0.2599, 0.2800, 0.3000,
  64                              0.3100, 0.3249, 0.3400, 0.3549, 0.3600, 0.3700, 0.3799, 0.3899, 0.4000, 0.4099};
  65   Float_t aQeCsI2[kNbins] = {0.0002, 0.0003, 0.0005, 0.0010, 0.0020, 0.0020, 0.0099, 0.0500, 0.1800, 0.1850,
  66                              0.1899, 0.1920, 0.1979, 0.2099, 0.2179, 0.2249, 0.2349, 0.2599, 0.2800, 0.3000,
  67                              0.3100, 0.3249, 0.3400, 0.3549, 0.3600, 0.3700, 0.3799, 0.3899, 0.4000, 0.4099};
  68   Float_t aQeCsI3[kNbins] = {0.0002, 0.0003, 0.0005, 0.0010, 0.0020, 0.0020, 0.0099, 0.0500, 0.1800, 0.1850,
  69                              0.1899, 0.1920, 0.1979, 0.2099, 0.2179, 0.2249, 0.2349, 0.2599, 0.2800, 0.3000,
  70                              0.3100, 0.3249, 0.3400, 0.3549, 0.3600, 0.3700, 0.3799, 0.3899, 0.4000, 0.4099};
  71   Float_t aQeCsI4[kNbins] = {0.0002, 0.0003, 0.0005, 0.0010, 0.0020, 0.0020, 0.0099, 0.0500, 0.1800, 0.1850,
  72                              0.1899, 0.1920, 0.1979, 0.2099, 0.2179, 0.2249, 0.2349, 0.2599, 0.2800, 0.3000,
  73                              0.3100, 0.3249, 0.3400, 0.3549, 0.3600, 0.3700, 0.3799, 0.3899, 0.4000, 0.4099};
  74   Float_t aQeCsI5[kNbins] = {0.0002, 0.0003, 0.0005, 0.0010, 0.0020, 0.0020, 0.0099, 0.0500, 0.1800, 0.1850,
  75                              0.1899, 0.1920, 0.1979, 0.2099, 0.2179, 0.2249, 0.2349, 0.2599, 0.2800, 0.3000,
  76                              0.3100, 0.3249, 0.3400, 0.3549, 0.3600, 0.3700, 0.3799, 0.3899, 0.4000, 0.4099};
  77   Float_t aQeCsI6[kNbins] = {0.0002, 0.0003, 0.0005, 0.0010, 0.0020, 0.0020, 0.0099, 0.0500, 0.1800, 0.1850,
  78                              0.1899, 0.1920, 0.1979, 0.2099, 0.2179, 0.2249, 0.2349, 0.2599, 0.2800, 0.3000,
  79                              0.3100, 0.3249, 0.3400, 0.3549, 0.3600, 0.3700, 0.3799, 0.3899, 0.4000, 0.4099};
  80   Float_t aQeCsI7[kNbins] = {0.0002, 0.0003, 0.0005, 0.0010, 0.0020, 0.0020, 0.0099, 0.0500, 0.1800, 0.1850,
  81                              0.1899, 0.1920, 0.1979, 0.2099, 0.2179, 0.2249, 0.2349, 0.2599, 0.2800, 0.3000,
  82                              0.3100, 0.3249, 0.3400, 0.3549, 0.3600, 0.3700, 0.3799, 0.3899, 0.4000, 0.4099};
  83
  84   Float_t aQeAll[kNbins];
  85   for(i=0;i<kNbins;i++){
  86     aQeCsI1[i]/= (1.0-Fresnel(aPckov[i]*1e9,1.0,0)); //FRESNEL LOSS CORRECTION
  87     aQeCsI2[i]/= (1.0-Fresnel(aPckov[i]*1e9,1.0,0)); //FRESNEL LOSS CORRECTION
  88     aQeCsI3[i]/= (1.0-Fresnel(aPckov[i]*1e9,1.0,0)); //FRESNEL LOSS CORRECTION
  89     aQeCsI4[i]/= (1.0-Fresnel(aPckov[i]*1e9,1.0,0)); //FRESNEL LOSS CORRECTION
  90     aQeCsI5[i]/= (1.0-Fresnel(aPckov[i]*1e9,1.0,0)); //FRESNEL LOSS CORRECTION
  91     aQeCsI6[i]/= (1.0-Fresnel(aPckov[i]*1e9,1.0,0)); //FRESNEL LOSS CORRECTION
  92     aQeCsI7[i]/= (1.0-Fresnel(aPckov[i]*1e9,1.0,0)); //FRESNEL LOSS CORRECTION
  93     aQeAll[i]=1; //QE for all other materials except for PC must be 1.
  94   }
  95
  96
  97 #ifdef __CINT__
  98
  99 //Now plot all the thigs
 100 //Freon, Quartz, Opaque ,Methane,CsI,Grid
 101   const Int_t kFreonMarker= 24; const Int_t kFreonColor=kRed;
 102   const Int_t kCH4Marker= 25;   const Int_t kCH4Color=kGreen;
 103   const Int_t kSiO2Marker= 26;  const Int_t kSiO2Color=kBlue;
 104   const Int_t kCsIMarker = 2;   const Int_t kCsIColor =kMagenta;
 105
 106   TCanvas *pC=new TCanvas("c1","RICH optics to check",1100,900);
 107
 108   pC->Divide(3,2);
 109
 110   pC->cd(1);
 111   TGraph *pAbsFreonGr=new TGraph(kNbins,aPckov,aAbsFreon);
 112   pAbsFreonGr->SetMarkerStyle(kFreonMarker); pAbsFreonGr->SetMarkerColor(kFreonColor);
 113   TGraph *pAbsSiO2Gr=new TGraph(kNbins,aPckov,aAbsSiO2);
 114   pAbsSiO2Gr->SetMarkerStyle(kSiO2Marker); pAbsSiO2Gr->SetMarkerColor(kSiO2Color);
 115   TMultiGraph *pAbsMG=new TMultiGraph();
 116   TLegend *pAbsLegend=new TLegend(0.6,0.3,0.85,0.5);
 117   pAbsMG->Add(pAbsFreonGr);     pAbsLegend->AddEntry(pAbsFreonGr,   "freon","p");            //1
 118   pAbsMG->Add(pAbsSiO2Gr);      pAbsLegend->AddEntry(pAbsSiO2Gr,  "SiO2","p");          //2
 119   pAbsMG->Draw("APL");
 120   pAbsMG->GetXaxis()->SetTitle("energy, GeV");
 121   pAbsMG->GetYaxis()->SetTitle("absorption length, cm");
 122   pAbsMG->Draw("APL");
 123   pAbsLegend->Draw();
 124
 125
 126   pC->cd(2);
 127   TGraph *pIndexFreonGr=new TGraph(kNbins,aPckov,aIndexFreon);
 128   pIndexFreonGr->SetMarkerStyle(kFreonMarker); pIndexFreonGr->SetMarkerColor(kFreonColor);
 129   TGraph *pIndexSiO2Gr=new TGraph(kNbins,aPckov,aIndexSiO2);
 130   pIndexSiO2Gr->SetMarkerStyle(kSiO2Marker); pIndexSiO2Gr->SetMarkerColor(kSiO2Color);
 131   TGraph *pIndexCH4Gr=new TGraph(kNbins,aPckov,aIndexCH4);
 132   pIndexCH4Gr->SetMarkerStyle(kCH4Marker); pIndexCH4Gr->SetMarkerColor(kCH4Color);
 133   TMultiGraph *pIndexMG=new TMultiGraph();
 134   TLegend *pIndexLegend=new TLegend(0.6,0.2,0.85,0.4);
 135   pIndexMG->Add(pIndexFreonGr);     pIndexLegend->AddEntry(pIndexFreonGr,   "C6F14","p");            //1
 136   pIndexMG->Add(pIndexSiO2Gr);      pIndexLegend->AddEntry(pIndexSiO2Gr,  "SiO2","p");
 137   pIndexMG->Add(pIndexCH4Gr);       pIndexLegend->AddEntry(pIndexCH4Gr,   "CH4","p");
 138   pIndexMG->Draw("APL");
 139   pIndexMG->GetXaxis()->SetTitle("energy, GeV");
 140   pIndexMG->GetYaxis()->SetTitle("refraction index");
 141   pIndexMG->Draw("APL");
 142   pIndexLegend->Draw();
 143
 144   pC->cd(3);
 145   TGraph *pAbsCH4Gr=new TGraph(kNbins,aPckov,aAbsCH4);
 146   pAbsCH4Gr->SetMarkerStyle(kCH4Marker); pAbsCH4Gr->SetMarkerColor(kCH4Color);
 147   pAbsCH4Gr->Draw("APL");
 148   pAbsCH4Gr->GetXaxis()->SetTitle("energy, GeV");
 149   pAbsCH4Gr->SetTitle("CH4 absorption length, cm");
 150   pAbsCH4Gr->Draw("APL");
 151
 152   pC->cd(4);
 153   TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI1);
 154   pQeCsIGr->SetMarkerStyle(kCsIMarker); pQeCsIGr->SetMarkerColor(kCsIColor);
 155   pQeCsIGr->Draw("APL");
 156   pQeCsIGr->GetXaxis()->SetTitle("energy, GeV");
 157   pQeCsIGr->SetTitle("QE");
 158   pQeCsIGr->Draw("APL");
 159
 160
 161   pC->cd(5);
 162   Float_t aTrFreon[kNbins],aTrSiO2[kNbins],aTrCH4[kNbins];
 163   Float_t aTotTr[kNbins];
 164   for(Int_t i=0;i<kNbins;i++){
 165     aTrFreon[i]=TMath::Exp(-AliRICHParam::FreonThickness() /(aAbsFreon[i]+0.0001));
 166     aTrSiO2[i] =TMath::Exp(-AliRICHParam::QuartzThickness()/(aAbsSiO2[i] +0.0001));
 167     aTrCH4[i]  =TMath::Exp(-AliRICHParam::GapThickness()   /(aAbsCH4[i]  +0.0001));
 168     aTotTr[i]    =aTrFreon[i]*aTrSiO2[i]*aTrCH4[i]*aQeCsI1[i];
 169   }
 170   TGraph *pTrFreonG=new TGraph(kNbins,aPckov,aTrFreon);pTrFreonG->SetMarkerStyle(kFreonMarker);pTrFreonG->SetMarkerColor(kFreonColor);
 171   TGraph *pTrSiO2G=new TGraph(kNbins,aPckov,aTrSiO2);pTrSiO2G->SetMarkerStyle(kSiO2Marker); pTrSiO2G->SetMarkerColor(kSiO2Color);
 172   TGraph *pTrCH4G=new TGraph(kNbins,aPckov,aTrCH4);pTrCH4G->SetMarkerStyle(kCH4Marker); pTrCH4G->SetMarkerColor(kCH4Color);
 173   TGraph *pTotTrG=new TGraph(kNbins,aPckov,aTotTr);pTotTrG->SetMarkerStyle(30);
 174   TMultiGraph *pTrMG=new TMultiGraph();
 175   TLegend *pTrLegend=new TLegend(0.6,0.2,0.85,0.4);
 176   pTrMG->Add(pTrFreonG);     pTrLegend->AddEntry(pTrFreonG, "freon","p");            //1
 177   pTrMG->Add(pTrSiO2G);      pTrLegend->AddEntry(pTrSiO2G,  "SiO2","p");
 178   pTrMG->Add(pTrCH4G);       pTrLegend->AddEntry(pTrCH4G,   "CH4","p");
 179   pTrMG->Add(pTotTrG);       pTrLegend->AddEntry(pTotTrG,   "total","p");
 180   pTrMG->Draw("APL");
 181   pTrMG->GetXaxis()->SetTitle("energy, GeV");
 182   pTrMG->GetYaxis()->SetTitle("transmission");
 183   pTrMG->Draw("APL");
 184   pTrLegend->Draw();
 185
 186
 187   return;
 188   TCanvas *pQeC=new TCanvas("pQeC","CsI QE currently all the same",800,900);
 189   pQeC->Divide(2,4);
 190   for(int i=1;i<=7;i++){
 191     pQeC->cd(i);
 192     switch(i){
 193       case 1: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI1);pQeCsIGr->SetTitle("Module 1");break;
 194       case 2: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI2);pQeCsIGr->SetTitle("Module 2");break;
 195       case 3: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI3);pQeCsIGr->SetTitle("Module 3");break;
 196       case 4: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI4);pQeCsIGr->SetTitle("Module 4");break;
 197       case 5: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI5);pQeCsIGr->SetTitle("Module 5");break;
 198       case 6: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI6);pQeCsIGr->SetTitle("Module 6");break;
 199       case 7: TGraph *pQeCsIGr=new TGraph(kNbins,aPckov,aQeCsI7);pQeCsIGr->SetTitle("Module 7");break;
 200     }
 201     pQeCsIGr->SetMarkerStyle(kCsIMarker); pQeCsIGr->SetMarkerColor(kCsIColor);
 202     pQeCsIGr->Draw("APL");
 203     pQeCsIGr->GetXaxis()->SetTitle("energy, GeV");
 204     pQeCsIGr->Draw("APL");
 205   }
 206 }//main
 207
 208 //__________________________________________________________________________________________________
 209 Float_t AbsoCH4(Float_t x)
 210 {//Evaluate the absorbtion lenght of CH4
 211   Float_t sch4[9] = {.12,.16,.23,.38,.86,2.8,7.9,28.,80.};              //MB X 10^22
 212   Float_t em[9] = {8.1,8.158,8.212,8.267,8.322,8.378,8.435,8.493,8.55};
 213   const Float_t kLoschmidt=2.686763e19;                                      // LOSCHMIDT NUMBER IN CM-3
 214   const Float_t kPressure=750.,kTemperature=283.;
 215   const Float_t pn=kPressure/760.;
 216   const Float_t tn=kTemperature/273.16;
 217   const Float_t c0=-1.655279e-1;
 218   const Float_t c1=6.307392e-2;
 219   const Float_t c2=-8.011441e-3;
 220   const Float_t c3=3.392126e-4;
 221
 222   Float_t crossSection=0;
 223   if (x<7.75)
 224     crossSection=.06e-22;
 225   else if(x>=7.75 && x<=8.1){                 //------ METHANE CROSS SECTION cm-2 ASTROPH. J. 214, L47 (1978)
 226         crossSection=(c0+c1*x+c2*x*x+c3*x*x*x)*1.e-18;
 227   }else if (x> 8.1){
 228     Int_t j=0;
 229     while (x<=em[j] || x>=em[j+1]){
 230       j++;
 231       Float_t a=(sch4[j+1]-sch4[j])/(em[j+1]-em[j]);
 232       crossSection=(sch4[j]+a*(x-em[j]))*1e-22;
 233     }
 234   }//if
 235
 236     Float_t density=kLoschmidt*pn/tn; //CH4 molecular density 1/cm-3
 237     return 1./(density*crossSection);
 238 }//AbsoCH4()
 239 //__________________________________________________________________________________________________
 240 Float_t Fresnel(Float_t ene,Float_t pdoti, Bool_t pola)
 241 {//ENE(EV), PDOTI=COS(INC.ANG.), PDOTR=COS(POL.PLANE ROT.ANG.)
 242
 243     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,
 244                       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,
 245                       7.8,7.9,8.0,8.1,8.2,8.3,8.4,8.5};
 246
 247
 248     Float_t csin[36] = {2.14,2.21,2.33,2.48,2.76,2.97,2.99,2.59,2.81,3.05,
 249                         2.86,2.53,2.55,2.66,2.79,2.96,3.18,3.05,2.84,2.81,2.38,2.11,
 250                         2.01,2.13,2.39,2.73,3.08,3.15,2.95,2.73,2.56,2.41,2.12,1.95,
 251                         1.72,1.53};
 252
 253     Float_t csik[36] = {0.,0.,0.,0.,0.,0.196,0.408,0.208,0.118,0.49,0.784,0.543,
 254                         0.424,0.404,0.371,0.514,0.922,1.102,1.139,1.376,1.461,1.253,0.878,
 255                         0.69,0.612,0.649,0.824,1.347,1.571,1.678,1.763,1.857,1.824,1.824,
 256                         1.714,1.498};
 257     Float_t xe=ene;
 258     Int_t  j=Int_t(xe*10)-49;
 259     Float_t cn=csin[j]+((csin[j+1]-csin[j])/0.1)*(xe-en[j]);
 260     Float_t ck=csik[j]+((csik[j+1]-csik[j])/0.1)*(xe-en[j]);
 261
 262     //FORMULAE FROM HANDBOOK OF OPTICS, 33.23 OR
 263     //W.R. HUNTER, J.O.S.A. 54 (1964),15 , J.O.S.A. 55(1965),1197
 264
 265     Float_t sinin=TMath::Sqrt(1-pdoti*pdoti);
 266     Float_t tanin=sinin/pdoti;
 267
 268     Float_t c1=cn*cn-ck*ck-sinin*sinin;
 269     Float_t c2=4*cn*cn*ck*ck;
 270     Float_t aO=TMath::Sqrt(0.5*(TMath::Sqrt(c1*c1+c2)+c1));
 271     Float_t b2=0.5*(TMath::Sqrt(c1*c1+c2)-c1);
 272
 273     Float_t rs=((aO-pdoti)*(aO-pdoti)+b2)/((aO+pdoti)*(aO+pdoti)+b2);
 274     Float_t rp=rs*((aO-sinin*tanin)*(aO-sinin*tanin)+b2)/((aO+sinin*tanin)*(aO+sinin*tanin)+b2);
 275
 276
 277     //CORRECTION FACTOR FOR SURFACE ROUGHNESS
 278     //B.J. STAGG  APPLIED OPTICS, 30(1991),4113
 279
 280     Float_t sigraf=18.;
 281     Float_t lamb=1240/ene;
 282     Float_t fresn;
 283
 284     Float_t  rO=TMath::Exp(-(4*TMath::Pi()*pdoti*sigraf/lamb)*(4*TMath::Pi()*pdoti*sigraf/lamb));
 285
 286     if(pola)
 287     {
 288         Float_t pdotr=0.8;                                 //DEGREE OF POLARIZATION : 1->P , -1->S
 289         fresn=0.5*(rp*(1+pdotr)+rs*(1-pdotr));
 290     }
 291     else
 292         fresn=0.5*(rp+rs);
 293
 294     fresn = fresn*rO;
 295     return(fresn);
 296 }//Fresnel(...)
 297 #endif