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