RICH/AliRICHv1.cxx

   1 // **************************************************************************
   2 // * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 // *                                                                        *
   4 // * Author: The ALICE Off-line Project.                                    *
   5 // * Contributors are mentioned in the code where appropriate.              *
   6 // *                                                                        *
   7 // * Permission to use, copy, modify and distribute this software and its   *
   8 // * documentation strictly for non-commercial purposes is hereby granted   *
   9 // * without fee, provided that the above copyright notice appears in all   *
  10 // * copies and that both the copyright notice and this permission notice   *
  11 // * appear in the supporting documentation. The authors make no claims     *
  12 // * about the suitability of this software for any purpose. It is          *
  13 // * provided "as is" without express or implied warranty.                  *
  14 // **************************************************************************
  15
  16
  17 #include "AliRICHv1.h"     //class header
  18 #include "AliRICHParam.h"
  19 #include <TParticle.h>     //Hits2SDigits()
  20 #include <TRandom.h>
  21 #include <TVirtualMC.h>    //StepManager() for gMC
  22 #include <TPDGCode.h>      //StepHistory()
  23 #include <AliStack.h>      //StepManager(),Hits2SDigits()
  24 #include <AliLoader.h>        //Hits2SDigits()
  25 #include <AliRunLoader.h>     //Hits2SDigits()
  26 #include <AliConst.h>
  27 #include <AliPDG.h>
  28 #include <AliMC.h>            //StepManager()
  29 #include <AliRawDataHeader.h> //Digits2Raw()
  30 #include <AliRun.h>           //CreateMaterials()
  31 #include <AliMagF.h>          //CreateMaterials()
  32 #include <TGeoManager.h>      //CreateGeometry()
  33 #include <TMultiGraph.h>      //Optics()
  34 #include <TGraph.h>           //Optics()
  35 #include <TLegend.h>          //Optics()
  36 #include <TCanvas.h>          //Optics()
  37 #include <TF2.h>              //CreateMaterials()
  38 #include <AliCDBManager.h>    //CreateMaterials()
  39 #include <AliCDBEntry.h>      //CreateMaterials()
  40
  41 ClassImp(AliRICHv1)
  42
  43 void AliRICHv1::CreateMaterials()
  44 {
  45 // Definition of available RICH materials
  46 // Arguments: none
  47 //   Returns: none
  48   AliDebug(1,"Start v1 RICH.");
  49
  50   const Int_t kNbins=30;       //number of photon energy points
  51
  52   Float_t aAbsC6F14[kNbins]={//New values from A.DiMauro 28.10.03 total 30
  53     32701.4219, 17996.1141, 10039.7281, 1799.1230, 1799.1231, 1799.1231, 1241.4091, 179.0987, 179.0986, 179.0987,
  54       179.0987,   118.9800,    39.5058,   23.7244,   11.1283,    7.1573,    3.6249,   2.1236,   0.7362,   0.5348,
  55         0.3387,     0.3074,     0.3050,    0.0001,    0.0001,    0.0001,    0.0001,   0.0001,   0.0001,   0.0001};
  56
  57   Float_t aAbsSiO2[kNbins]={//New values from A.DiMauro 28.10.03 total 30
  58        34.4338, 30.5424, 30.2584, 31.4928, 31.7868, 17.8397, 9.3410, 6.4492, 6.1128, 5.8128,
  59         5.5589,  5.2877,  5.0162,  4.7999,  4.5734,  4.2135, 3.7471, 2.6033, 1.5223, 0.9658,
  60         0.4242,  0.2500,  0.1426,  0.0863,  0.0793,  0.0724, 0.0655, 0.0587, 0.0001, 0.0001};
  61
  62   Float_t aQeCsI[kNbins] = {//New values from A.DiMauro 28.10.03 total 31 the last one cut to provide 30
  63                             0.0002, 0.0006, 0.0007, 0.0010, 0.0049, 0.0073, 0.0104, 0.0519, 0.0936, 0.1299,
  64                             0.1560, 0.1768, 0.1872, 0.1976, 0.2142, 0.2288, 0.2434, 0.2599, 0.2673, 0.2808,
  65                             0.2859, 0.2954, 0.3016, 0.3120, 0.3172, 0.3224, 0.3266, 0.3328, 0.3359, 0.3390}; //0.3431};
  66
  67   Float_t aQeCsIold[kNbins]={//previous values 26 in total added 0.0001 to the 30
  68     0.0002, 0.0006, 0.0007, 0.0050, 0.0075, 0.0101, 0.0243, 0.0405, 0.0689, 0.1053,
  69     0.1215, 0.1417, 0.1579, 0.1620, 0.1661, 0.1677, 0.1743, 0.1768, 0.1793, 0.1826,
  70     0.1859, 0.1876, 0.1892, 0.1909, 0.2075, 0.2158, 0.0001, 0.0001, 0.0001, 0.0001 };
  71
  72 //         radiator            window             gas               metal
  73   Float_t aIdxC6F14[kNbins] , aIdxSiO2[kNbins] , aIdxCH4[kNbins] , aIdx0[kNbins]   , aIdx1[kNbins] ;
  74   Float_t                                        aAbsCH4[kNbins] , aAbsMet[kNbins]                 ;
  75   Float_t                                                                            aQe1[kNbins]  ;    //QE for all but PC
  76   Float_t aEckov[kNbins];  //Ckov energy in GeV
  77
  78 //Read all the staff from CDB
  79   TF2 *pIdxC6F14;
  80   AliCDBEntry *pEntry=AliCDBManager::Instance()->Get("RICH/RICHConfig/RefIdxC6F14",0,0,0); //0-0-0 is for simulation
  81   if(pEntry){
  82     pIdxC6F14=(TF2*)pEntry->GetObject(); delete pEntry;
  83   }else{
  84     AliWarning("No valid calibarion, the hardcoded will be used!");
  85     pIdxC6F14=new TF2("RidxC4F14","sqrt(1+0.554*(1239.84e-9/x)^2/((1239.84e-9/x)^2-5796)-0.0005*(y-20))",5.5e-9,8.5e-9,0,50); //DiMauro mail
  86   }
  87
  88   Double_t eCkovMin=5.5e-9,eCkovMax=8.5e-9; //in GeV
  89   TF1 idxSiO2("RidxSiO2" ,"sqrt(1+46.411/(10.666*10.666-x*x*1e18)+228.71/(18.125*18.125-x*x*1e18))",eCkovMin,eCkovMax);  //TDR p.35
  90
  91
  92   for(Int_t i=0;i<kNbins;i++){
  93     aEckov     [i] =AliRICHParam::EckovMin()+0.1e-9*i;//Ckov energy in GeV
  94
  95     aIdxC6F14  [i] =pIdxC6F14->Eval(aEckov[i],pIdxC6F14->GetUniqueID());
  96     aIdxSiO2   [i] =idxSiO2   .Eval(aEckov[i]);
  97     aIdxCH4    [i] =AliRICHParam::IdxCH4   (aEckov[i]); aAbsCH4 [i]  =AliRICHParam::AbsCH4     (aEckov[i]);
  98     aAbsMet    [i] =0.0001;                              //metal has absorption probability
  99     aIdx0      [i] =0;                                   //metal ref idx must be 0 in order to reflect photon
 100     aIdx1      [i] =1;                                   //metal ref idx must be 1 in order to apply photon to QE conversion
 101 //    aQeCsI[i]     /= (1.0-Fresnel(aEckov[i]*1e9,1.0,0)) ; //FRESNEL LOSS CORRECTION ?????????????
 102     aQe1       [i] =1                                   ; //QE for all other materials except for PC must be 1.
 103   }
 104
 105 //data from PDG booklet 2002     density [gr/cm^3] rad len [cm] abs len [cm]
 106   Float_t   aAir[4]={12,14,16,36}    ,   zAir[4]={6,7,8,18} ,   wAir[4]={0.000124,0.755267,0.231781,0.012827} , dAir=0.00120479; Int_t nAir=4;//mixture 0.9999999
 107   Float_t aC6F14[2]={ 12.01 , 18.99} , zC6F14[2]={ 6 , 9}   , wC6F14[2]={6 , 14} , dC6F14=1.68    ; Int_t nC6F14=-2;
 108   Float_t  aSiO2[2]={ 28.09 , 15.99} ,  zSiO2[2]={14 , 8}   ,  wSiO2[2]={1 ,  2} ,  dSiO2=2.64    ; Int_t  nSiO2=-2;
 109   Float_t   aCH4[2]={ 12.01 ,  1.01} ,   zCH4[2]={ 6 , 1}   ,   wCH4[2]={1 ,  4} ,   dCH4=7.17e-4 ; Int_t   nCH4=-2;
 110   Float_t   aCsI[2]={132.90 ,126.90} ,   zCsI[2]={55 ,53}   ,   wCsI[2]={1 ,  1} ,   dCsI=0.1     ; Int_t   nCsI=-2;
 111   Float_t     aRoha= 12.01 ,               zRoha=  6 ,                              dRoha=  0.10 , radRoha= 18.80 , absRoha=  86.3/dRoha; //special material- quazi carbon
 112   Float_t       aCu= 63.55 ,                 zCu= 29 ,                                dCu=  8.96 ,   radCu=  1.43 ,   absCu= 134.9/dCu  ;
 113   Float_t        aW=183.84 ,                  zW= 74 ,                                 dW= 19.30 ,    radW=  0.35 ,    absW= 185.0/dW   ;
 114   Float_t       aAl= 26.98 ,                 zAl= 13 ,                                dAl=  2.70 ,   radAl=  8.90 ,   absAl= 106.4/dAl  ;
 115
 116   Int_t   matId=0;                           //tmp material id number
 117   Int_t   unsens =  0, sens=1;               //sensitive or unsensitive medium
 118   Int_t   itgfld = gAlice->Field()->Integ(); //type of field intergration 0 no field -1 user in guswim 1 Runge Kutta 2 helix 3 const field along z
 119   Float_t maxfld = gAlice->Field()->Max();   //max field value
 120   Float_t tmaxfd = -10.0;                    //max deflection angle due to magnetic field in one step
 121   Float_t deemax = - 0.2;                    //max fractional energy loss in one step
 122   Float_t stemax = - 0.1;                    //mas step allowed [cm]
 123   Float_t epsil  =   0.001;                  //abs tracking precision [cm]
 124   Float_t stmin  = - 0.001;                  //min step size [cm] in continius process transport, negative value: choose it automatically
 125   AliMixture(++matId,"Air"  ,aAir  ,zAir  ,dAir  ,nAir  ,wAir  ); AliMedium(kAir  ,"Air"  ,matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
 126   AliMixture(++matId,"C6F14",aC6F14,zC6F14,dC6F14,nC6F14,wC6F14); AliMedium(kC6F14,"C6F14",matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
 127   AliMixture(++matId,"SiO2" ,aSiO2 ,zSiO2 ,dSiO2 ,nSiO2 ,wSiO2 ); AliMedium(kSiO2 ,"SiO2" ,matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
 128   AliMixture(++matId,"CH4"  ,aCH4  ,zCH4  ,dCH4  ,nCH4  ,wCH4  ); AliMedium(kCH4  ,"CH4"  ,matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
 129   AliMixture(++matId,"CsI"  ,aCsI  ,zCsI  ,dCsI  ,nCsI  ,wCsI  ); AliMedium(kCsI  ,"CsI"  ,matId,   sens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);//sensitive
 130
 131   AliMaterial(++matId,"Roha",aRoha,zRoha,dRoha,radRoha,absRoha);  AliMedium(kRoha,"Roha", matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
 132   AliMaterial(++matId,"Cu"  ,aCu  ,zCu  ,dCu  ,radCu  ,absCu  );  AliMedium(kCu  ,"Cu"  , matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
 133   AliMaterial(++matId,"W"   ,aW   ,zW   ,dW   ,radW   ,absW   );  AliMedium(kW   ,"W"   , matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
 134   AliMaterial(++matId,"Al"  ,aAl  ,zAl  ,dAl  ,radAl  ,absAl  );  AliMedium(kAl  ,"Al"  , matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
 135
 136   gMC->SetCerenkov((*fIdtmed)[kC6F14]    , kNbins, aEckov, aAbsC6F14  , aQe1  , aIdxC6F14 );
 137   gMC->SetCerenkov((*fIdtmed)[kSiO2]     , kNbins, aEckov, aAbsSiO2   , aQe1  , aIdxSiO2  );
 138   gMC->SetCerenkov((*fIdtmed)[kCH4]      , kNbins, aEckov, aAbsCH4    , aQe1  , aIdxCH4   );
 139   gMC->SetCerenkov((*fIdtmed)[kCu]       , kNbins, aEckov, aAbsMet    , aQe1  , aIdx0     );
 140   gMC->SetCerenkov((*fIdtmed)[kW]        , kNbins, aEckov, aAbsMet    , aQe1  , aIdx0     ); //n=0 means reflect photons
 141   gMC->SetCerenkov((*fIdtmed)[kCsI]      , kNbins, aEckov, aAbsMet    , aQeCsI, aIdx1     ); //n=1 means convert photons
 142   gMC->SetCerenkov((*fIdtmed)[kAl]       , kNbins, aEckov, aAbsMet    , aQe1  , aIdx0     );
 143
 144   AliDebug(1,"Stop v1 RICH.");
 145   TString ttl=GetTitle();
 146   if(!ttl.Contains("ShowOptics")) return; //do not plot optical curves
 147   {
 148   const Double_t kWidth=0.25,kHeight=0.2;
 149   const Int_t kC6F14Marker=24 , kC6F14Color=kRed;
 150   const Int_t kCH4Marker  =25 , kCH4Color  =kGreen;
 151   const Int_t kSiO2M      =26 , kSiO2Color =kBlue;
 152   const Int_t kCsIMarker  = 2 , kCsIColor  =kMagenta;
 153
 154 //Ref index
 155   TGraph *pIdxC6F14=new TGraph(kNbins,aEckov,aIdxC6F14);pIdxC6F14->SetMarkerStyle(kC6F14Marker);pIdxC6F14->SetMarkerColor(kC6F14Color);
 156   TGraph *pIdxSiO2 =new TGraph(kNbins,aEckov,aIdxSiO2); pIdxSiO2 ->SetMarkerStyle(kSiO2M)      ;pIdxSiO2 ->SetMarkerColor(kSiO2Color);
 157   TGraph *pIdxCH4  =new TGraph(kNbins,aEckov,aIdxCH4);  pIdxCH4  ->SetMarkerStyle(kCH4Marker)  ;pIdxCH4  ->SetMarkerColor(kCH4Color);
 158   TMultiGraph *pIdxMG=new TMultiGraph("refidx","Ref index;E_{#check{C}} [GeV]");  TLegend *pIdxLe=new TLegend(0.5,0.21,0.5+kWidth,0.21+kHeight);
 159   pIdxMG->Add(pIdxC6F14); pIdxLe->AddEntry(pIdxC6F14,"C6F14"  ,"p");
 160   pIdxMG->Add(pIdxSiO2) ; pIdxLe->AddEntry(pIdxSiO2 ,"SiO2"   ,"p");
 161   pIdxMG->Add(pIdxCH4)  ; pIdxLe->AddEntry(pIdxCH4  ,"CH4"    ,"p");
 162 //Absorbtion
 163   TGraph *pAbsC6F14=new TGraph(kNbins,aEckov,aAbsC6F14);pAbsC6F14->SetMarkerStyle(kC6F14Marker); pAbsC6F14->SetMarkerColor(kC6F14Color);
 164   TGraph *pAbsSiO2 =new TGraph(kNbins,aEckov,aAbsSiO2) ;pAbsSiO2 ->SetMarkerStyle(kSiO2M)      ; pAbsSiO2 ->SetMarkerColor(kSiO2Color);
 165   TGraph *pAbsCH4  =new TGraph(kNbins,aEckov,aAbsCH4)  ;pAbsCH4  ->SetMarkerStyle(kCH4Marker)  ; pAbsCH4  ->SetMarkerColor(kCH4Color);
 166
 167   TMultiGraph *pAbsMG=new TMultiGraph("abs","Absorption [cm];E_{#check{C}} [GeV]");  TLegend *pAbsLe=new TLegend(0.2,0.15,0.2+kWidth,0.15+kHeight);
 168   pAbsMG->Add(pAbsC6F14);      pAbsLe->AddEntry(pAbsC6F14,  "C6F14"    ,"p");
 169   pAbsMG->Add(pAbsSiO2);       pAbsLe->AddEntry(pAbsSiO2 ,  "SiO2"     ,"p");
 170   pAbsMG->Add(pAbsCH4);        pAbsLe->AddEntry(pAbsCH4  ,  "CH4"      ,"p");
 171 //QE new and old
 172   TGraph *pQeCsI   =new TGraph(kNbins,aEckov,aQeCsI);    pQeCsI   ->SetMarkerStyle(kCsIMarker);  pQeCsI   ->SetMarkerColor(kCsIColor);
 173   TGraph *pQeCsIold=new TGraph(kNbins,aEckov,aQeCsIold); pQeCsIold->SetMarkerStyle(kC6F14Marker);pQeCsIold->SetMarkerColor(kC6F14Color);
 174   TMultiGraph *pCompMG=new TMultiGraph("qe","QE;E_{#check{C}} [GeV]");  TLegend *pCompLe=new TLegend(0.2,0.6,0.2+kWidth,0.6+kHeight);
 175   pCompMG->Add(pQeCsI);       pCompLe->AddEntry(pQeCsI,    "QE new 30.10.03", "p");
 176   pCompMG->Add(pQeCsIold);    pCompLe->AddEntry(pQeCsIold, "QE old 01.01.02", "p");
 177 //transmission
 178   Float_t aTrC6F14[kNbins],aTrSiO2[kNbins],aTrCH4[kNbins],aTrTotal[kNbins];
 179   for(Int_t i=0;i<kNbins;i++){//calculate probability for photon to survive during transversing a volume of material with absorption length
 180     aTrC6F14[i] =TMath::Exp(-AliRICHParam::RadThick() / (aAbsC6F14[i]+0.0001)); //radiator
 181     aTrSiO2[i]  =TMath::Exp(-AliRICHParam::WinThick() / (aAbsSiO2[i] +0.0001)); //window
 182     aTrCH4[i]   =TMath::Exp(-AliRICHParam::Pc2Win()   / (aAbsCH4[i]  +0.0001)); //from window to PC
 183     aTrTotal[i] =aTrC6F14[i]*aTrSiO2[i]*aTrCH4[i]*aQeCsI[i];
 184   }
 185   TGraph *pTrC6F14=new TGraph(kNbins,aEckov,aTrC6F14)  ;pTrC6F14->SetMarkerStyle(kC6F14Marker);pTrC6F14->SetMarkerColor(kC6F14Color);
 186   TGraph *pTrSiO2 =new TGraph(kNbins,aEckov,aTrSiO2)   ;pTrSiO2 ->SetMarkerStyle(kSiO2M)      ;pTrSiO2 ->SetMarkerColor(kSiO2Color);
 187   TGraph *pTrCH4  =new TGraph(kNbins,aEckov,aTrCH4)    ;pTrCH4  ->SetMarkerStyle(kCH4Marker)  ;pTrCH4  ->SetMarkerColor(kCH4Color);
 188   TGraph *pTrTotal=new TGraph(kNbins,aEckov,aTrTotal)  ;pTrTotal->SetMarkerStyle(30)          ;pTrTotal->SetMarkerColor(kYellow);
 189   TMultiGraph *pTrMG=new TMultiGraph("trans","Transmission;E_{#check{C}} [GeV]");  TLegend *pTrLe=new TLegend(0.2,0.4,0.2+kWidth,0.4+kHeight);
 190   pTrMG->Add(pQeCsI);       pTrLe->AddEntry(pQeCsI,   "CsI QE", "p");
 191   pTrMG->Add(pTrC6F14);     pTrLe->AddEntry(pTrC6F14, "C6F14" , "p");
 192   pTrMG->Add(pTrSiO2);      pTrLe->AddEntry(pTrSiO2,  "SiO2"  , "p");
 193   pTrMG->Add(pTrCH4);       pTrLe->AddEntry(pTrCH4,   "CH4"   , "p");
 194   pTrMG->Add(pTrTotal);     pTrLe->AddEntry(pTrTotal, "total" , "p");
 195
 196   TCanvas *pC=new TCanvas("c1","RICH optics to check",1100,900);  pC->Divide(2,2);
 197   pC->cd(1);                    pIdxMG ->Draw("AP");  pIdxLe ->Draw();      //ref idx
 198   pC->cd(2);  gPad->SetLogy();  pAbsMG ->Draw("AP");  pAbsLe ->Draw();      //absorption
 199   pC->cd(3);                    pCompMG->Draw("AP");  pCompLe->Draw();      //QE
 200   pC->cd(4);                    pTrMG  ->Draw("AP");  pTrLe  ->Draw();      //transmission
 201   }
 202 }//void AliRICH::CreateMaterials()
 203 //__________________________________________________________________________________________________
 204 void AliRICHv1::CreateGeometry()
 205 {
 206 //Creates detailed geometry simulation (currently GEANT volumes tree)
 207   AliDebug(1,"Start main.");
 208   if(!gMC->IsRootGeometrySupported()) return;
 209
 210   Double_t cm=1,mm=0.1*cm,mkm=0.001*mm,dx,dy,dz;//default is cm
 211
 212   TGeoVolume *pRich=gGeoManager->MakeBox("RICH",gGeoManager->GetMedium("RICH_CH4"),dx=(6*mm+1681*mm+6*mm)/2,  //main RICH volume
 213                                                                                    dy=(6*mm+1466*mm+6*mm)/2,
 214                                                                                    dz=(80*mm+40*mm)*2/2);     //x,y taken from 2033P1  z from p84 TDR
 215   const Double_t kAngHor=19.5; //  horizontal angle between chambers  19.5 grad
 216   const Double_t kAngVer=20;   //  vertical angle between chambers    20   grad
 217   const Double_t kAngCom=30;   //  common RICH rotation with respect to x axis  30   grad
 218   const Double_t trans[3]={490,0,0}; //center of the chamber is on window-gap surface
 219   for(Int_t iCh=1;iCh<=7;iCh++){//place 7 chambers
 220     TGeoHMatrix *pMatrix=new TGeoHMatrix;
 221     pMatrix->RotateY(90);           //rotate around y since initial position is in XY plane -> now in YZ plane
 222     pMatrix->SetTranslation(trans); //now plane in YZ is shifted along x
 223     switch(iCh){
 224       case 1:                pMatrix->RotateY(kAngHor);  pMatrix->RotateZ(-kAngVer);  break; //right and down
 225       case 2:                                            pMatrix->RotateZ(-kAngVer);  break; //down
 226       case 3:                pMatrix->RotateY(kAngHor);                               break; //right
 227       case 4:                                                                         break; //no rotation
 228       case 5:                pMatrix->RotateY(-kAngHor);                              break; //left
 229       case 6:                                            pMatrix->RotateZ(kAngVer);   break; //up
 230       case 7:                pMatrix->RotateY(-kAngHor); pMatrix->RotateZ(kAngVer);   break; //left and up
 231     }
 232     pMatrix->RotateZ(kAngCom);     //apply common rotation  in XY plane
 233     gGeoManager->GetVolume("ALIC")->AddNode(pRich,iCh,pMatrix);
 234   }
 235
 236   Float_t par[3];
 237   Int_t matrixIdReturn=0; //matrix id returned by AliMatrix
 238 //Pad Panel frame  6 sectors
 239   par[0]=648*mm/2;par[1]=  411*mm/2;par[2]=40  *mm/2;gMC->Gsvolu("Rppf"     ,"BOX ",(*fIdtmed)[kAl]  ,par,3);//PPF 2001P2 inner size of the slab by 1mm more
 240   par[0]=181*mm/2;par[1]=89.25*mm/2;par[2]=38.3*mm/2;gMC->Gsvolu("RppfLarge","BOX ",(*fIdtmed)[kAir] ,par,3);//large whole
 241   par[0]=114*mm/2;par[1]=89.25*mm/2;par[2]=38.3*mm/2;gMC->Gsvolu("RppfSmall","BOX ",(*fIdtmed)[kAir] ,par,3);//small whole
 242   par[0]=644*mm/2;par[1]=  407*mm/2;par[2]= 1.7*mm/2;gMC->Gsvolu("Rpc"      ,"BOX ",(*fIdtmed)[kCsI] ,par,3);//by 0.2 mm more then actual size (PCB 2006P1)
 243
 244   gMC->Gspos("Rppf",1,"RICH",    -335*mm,      -433*mm,  8*cm+20*mm,  0,"ONLY");//F1 2040P1 z p.84 TDR
 245   gMC->Gspos("Rppf",2,"RICH",    +335*mm,      -433*mm,  8*cm+20*mm,  0,"ONLY");
 246   gMC->Gspos("Rppf",3,"RICH",    -335*mm,         0*mm,  8*cm+20*mm,  0,"ONLY");
 247   gMC->Gspos("Rppf",4,"RICH",    +335*mm,         0*mm,  8*cm+20*mm,  0,"ONLY");
 248   gMC->Gspos("Rppf",5,"RICH",    -335*mm,      +433*mm,  8*cm+20*mm,  0,"ONLY");
 249   gMC->Gspos("Rppf",6,"RICH",    +335*mm,      +433*mm,  8*cm+20*mm,  0,"ONLY");
 250     gMC->Gspos("Rpc"      ,1,"Rppf",       0*mm,         0*mm,   -19.15*mm,  0,"ONLY");//PPF 2001P2
 251     gMC->Gspos("RppfLarge",1,"Rppf",  -224.5*mm,  -151.875*mm,     0.85*mm,  0,"ONLY");
 252     gMC->Gspos("RppfLarge",2,"Rppf",  -224.5*mm,  - 50.625*mm,     0.85*mm,  0,"ONLY");
 253     gMC->Gspos("RppfLarge",3,"Rppf",  -224.5*mm,  + 50.625*mm,     0.85*mm,  0,"ONLY");
 254     gMC->Gspos("RppfLarge",4,"Rppf",  -224.5*mm,  +151.875*mm,     0.85*mm,  0,"ONLY");
 255     gMC->Gspos("RppfSmall",1,"Rppf",  - 65.0*mm,  -151.875*mm,     0.85*mm,  0,"ONLY");
 256     gMC->Gspos("RppfSmall",2,"Rppf",  - 65.0*mm,  - 50.625*mm,     0.85*mm,  0,"ONLY");
 257     gMC->Gspos("RppfSmall",3,"Rppf",  - 65.0*mm,  + 50.625*mm,     0.85*mm,  0,"ONLY");
 258     gMC->Gspos("RppfSmall",4,"Rppf",  - 65.0*mm,  +151.875*mm,     0.85*mm,  0,"ONLY");
 259     gMC->Gspos("RppfSmall",5,"Rppf",  + 65.0*mm,  -151.875*mm,     0.85*mm,  0,"ONLY");
 260     gMC->Gspos("RppfSmall",6,"Rppf",  + 65.0*mm,  - 50.625*mm,     0.85*mm,  0,"ONLY");
 261     gMC->Gspos("RppfSmall",7,"Rppf",  + 65.0*mm,  + 50.625*mm,     0.85*mm,  0,"ONLY");
 262     gMC->Gspos("RppfSmall",8,"Rppf",  + 65.0*mm,  +151.875*mm,     0.85*mm,  0,"ONLY");
 263     gMC->Gspos("RppfLarge",5,"Rppf",  +224.5*mm,  -151.875*mm,     0.85*mm,  0,"ONLY");
 264     gMC->Gspos("RppfLarge",6,"Rppf",  +224.5*mm,  - 50.625*mm,     0.85*mm,  0,"ONLY");
 265     gMC->Gspos("RppfLarge",7,"Rppf",  +224.5*mm,  + 50.625*mm,     0.85*mm,  0,"ONLY");
 266     gMC->Gspos("RppfLarge",8,"Rppf",  +224.5*mm,  +151.875*mm,     0.85*mm,  0,"ONLY");
 267 //Gap - anod wires 6 copies to RICH
 268   par[0]=648*mm/2;par[1]=  411*mm/2 ;par[2]=4.45*mm/2;gMC->Gsvolu("Rgap","BOX ",(*fIdtmed)[kCH4] ,par,3);//xy as PPF 2001P2 z WP 2099P1
 269   par[0]=  0*mm  ;par[1]=  20*mkm/2 ;par[2]= 648*mm/2;gMC->Gsvolu("Rano","TUBE",(*fIdtmed)[kW]   ,par,3);//WP 2099P1 z = gap x PPF 2001P2
 270   AliMatrix(matrixIdReturn,180,0, 90,90, 90,0); //wires along x
 271
 272   gMC->Gspos("Rgap",1,"RICH",    -335*mm,      -433*mm,8*cm-2.225*mm, 0,"ONLY"); //F1 2040P1 z WP 2099P1
 273   gMC->Gspos("Rgap",2,"RICH",    +335*mm,      -433*mm,8*cm-2.225*mm, 0,"ONLY");
 274   gMC->Gspos("Rgap",3,"RICH",    -335*mm,         0*mm,8*cm-2.225*mm, 0,"ONLY");
 275   gMC->Gspos("Rgap",4,"RICH",    +335*mm,         0*mm,8*cm-2.225*mm, 0,"ONLY");
 276   gMC->Gspos("Rgap",5,"RICH",    -335*mm,      +433*mm,8*cm-2.225*mm, 0,"ONLY");
 277   gMC->Gspos("Rgap",6,"RICH",    +335*mm,      +433*mm,8*cm-2.225*mm, 0,"ONLY");
 278   for(int i=1;i<=96;i++)
 279     gMC->Gspos("Rano",i,"Rgap",     0*mm, -411/2*mm+i*4*mm, 0.185*mm, matrixIdReturn,"ONLY"); //WP 2099P1
 280 //Defines radiators geometry
 281   par[0]=1330*mm/2 ;par[1]= 413*mm/2  ;par[2]=  24*mm/2;  gMC->Gsvolu("Rrad"      ,"BOX ",(*fIdtmed)[kC6F14]     ,par,3); // Rad 2011P1
 282   par[0]=1330*mm/2 ;par[1]= 413*mm/2  ;par[2]=   4*mm/2;  gMC->Gsvolu("RradFront" ,"BOX ",(*fIdtmed)[kRoha]      ,par,3); //front
 283   par[0]=1330*mm/2 ;par[1]= 413*mm/2  ;par[2]=   5*mm/2;  gMC->Gsvolu("RradWin"   ,"BOX ",(*fIdtmed)[kSiO2]      ,par,3); //window
 284   par[0]=1330*mm/2 ;par[1]=   5*mm/2  ;par[2]=  15*mm/2;  gMC->Gsvolu("RradLong"  ,"BOX ",(*fIdtmed)[kRoha]      ,par,3); //long side
 285   par[0]=  10*mm/2 ;par[1]= 403*mm/2  ;par[2]=  15*mm/2;  gMC->Gsvolu("RradShort" ,"BOX ",(*fIdtmed)[kRoha]      ,par,3); //short side
 286   par[0]=   0      ;par[1]=  10*mm/2  ;par[2]=  15*mm/2;  gMC->Gsvolu("RradSpacer","TUBE",(*fIdtmed)[kSiO2]      ,par,3); //spacer
 287
 288   gMC->Gspos("Rrad",1,"RICH",   0*mm,-434*mm,   -12*mm,  0,"ONLY"); //3 radiators to RICH
 289   gMC->Gspos("Rrad",2,"RICH",   0*mm,   0*mm,   -12*mm,  0,"ONLY");
 290   gMC->Gspos("Rrad",3,"RICH",   0*mm,+434*mm,   -12*mm,  0,"ONLY");
 291     gMC->Gspos("RradFront",1,"Rrad",   0*mm,   0*mm, -10.0*mm,  0,"ONLY"); //front cover
 292     gMC->Gspos("RradWin"  ,1,"Rrad",   0*mm,   0*mm,   9.5*mm,  0,"ONLY"); //quartz window (back cover)
 293     gMC->Gspos("RradLong" ,1,"Rrad",   0*mm,-204*mm,  -0.5*mm,  0,"ONLY"); //long side
 294     gMC->Gspos("RradLong" ,2,"Rrad",   0*mm,+204*mm,  -0.5*mm,  0,"ONLY"); //long side
 295     gMC->Gspos("RradShort",1,"Rrad",-660*mm,   0*mm,  -0.5*mm,  0,"ONLY"); //short side
 296     gMC->Gspos("RradShort",2,"Rrad",+660*mm,   0*mm,  -0.5*mm,  0,"ONLY"); //short side
 297     for(int i=0;i<3;i++)
 298       for(int j=0;j<10;j++)
 299         gMC->Gspos("RradSpacer",10*i+j,"Rrad",-1330*mm/2+116*mm+j*122*mm,(i-1)*105*mm,-0.5*mm,0,"ONLY");//spacers
 300 //Defines SandBox geometry
 301   par[0]=1419*mm/2 ;par[1]=1378*mm/2;par[2]=50.5*mm/2; gMC->Gsvolu("Rsb"     ,"BOX ",(*fIdtmed)[kAir]  ,par,3);  //2072P1
 302   par[0]=1419*mm/2 ;par[1]=1378*mm/2;par[2]= 0.5*mm/2; gMC->Gsvolu("RsbCover","BOX ",(*fIdtmed)[kAl]   ,par,3);  //cover
 303   par[0]=1359*mm/2 ;par[1]=1318*mm/2;par[2]=49.5*mm/2; gMC->Gsvolu("RsbComb" ,"BOX ",(*fIdtmed)[kRoha] ,par,3);  //honeycomb structure
 304
 305   gMC->Gspos("Rsb",1,"RICH",   0*mm, 0*mm, -73.75*mm, 0,"ONLY"); //p.84 TDR sandbox to rich
 306     gMC->Gspos("RsbComb" ,1,"Rsb", 0*mm, 0*mm,      0*mm, 0,"ONLY"); //2072P1 honeycomv to sandbox
 307     gMC->Gspos("RsbCover",1,"Rsb", 0*mm, 0*mm,    +25*mm, 0,"ONLY"); //cover to sandbox
 308     gMC->Gspos("RsbCover",2,"Rsb", 0*mm, 0*mm,    -25*mm, 0,"ONLY"); //cover to sandbox
 309   AliDebug(1,"Stop v1. HMPID option");
 310 }//CreateGeometry()
 311 //__________________________________________________________________________________________________
 312 void AliRICHv1::Init()
 313 {
 314 // This methode defines ID for sensitive volumes, i.e. such geometry volumes for which there are if(gMC->CurrentVolID()==XXX) statements in StepManager()
 315 // Arguments: none
 316 //   Returns: none
 317   AliDebug(1,"Start v1 HMPID.");
 318   fIdRad     = gMC->VolId("Rrad");
 319   fIdWin     = gMC->VolId("RradWin");
 320   fIdPc      = gMC->VolId("Rpc");
 321   fIdAmpGap  = gMC->VolId("Rgap");
 322   fIdProxGap = gMC->VolId("Rgap");
 323   AliDebug(1,"Stop v1 HMPID.");
 324 }
 325 //__________________________________________________________________________________________________
 326 Bool_t AliRICHv1::IsLostByFresnel()
 327 {
 328 // Calculate probability for the photon to be lost by Fresnel reflection.
 329   TLorentzVector p4;
 330   Double_t mom[3],localMom[3];
 331   gMC->TrackMomentum(p4);   mom[0]=p4(1);   mom[1]=p4(2);   mom[2]=p4(3);
 332   localMom[0]=0; localMom[1]=0; localMom[2]=0;
 333   gMC->Gmtod(mom,localMom,2);
 334   Double_t localTc    = localMom[0]*localMom[0]+localMom[2]*localMom[2];
 335   Double_t localTheta = TMath::ATan2(TMath::Sqrt(localTc),localMom[1]);
 336   Double_t cotheta = TMath::Abs(TMath::Cos(localTheta));
 337   if(gMC->GetRandom()->Rndm() < Fresnel(p4.E()*1e9,cotheta,1)){
 338     AliDebug(1,"Photon lost");
 339     return kTRUE;
 340   }else
 341     return kFALSE;
 342 }//IsLostByFresnel()
 343 //__________________________________________________________________________________________________
 344 void AliRICHv1::GenFee(Int_t iChamber,Float_t eloss)
 345 {
 346 // Generate FeedBack photons for the current particle. To be invoked from StepManager().
 347 // eloss=0 means photon so only pulse height distribution is to be analysed. This one is done in AliRICHParam::TotQdc()
 348
 349   TLorentzVector x4;
 350   gMC->TrackPosition(x4);
 351   TVector2 x2=AliRICHParam::Instance()->Mars2Lors(iChamber,x4.Vect(),AliRICHParam::kPc);//hit position on photocathode plane
 352   TVector2 xspe=x2;
 353   Int_t sector=AliRICHParam::Loc2Sec(xspe);  if(sector==-1) return; //hit in dead zone, nothing to produce
 354   Int_t iTotQdc=AliRICHParam::TotQdc(x2,eloss);
 355   Int_t iNphotons=gMC->GetRandom()->Poisson(AliRICHParam::AlphaFeedback(iChamber,sector)*iTotQdc);
 356   AliDebug(1,Form("N photons=%i",iNphotons));
 357   Int_t j;
 358   Float_t cthf, phif, enfp = 0, sthf, e1[3], e2[3], e3[3], vmod, uswop,dir[3], phi,pol[3], mom[4];
 359 //Generate photons
 360   for(Int_t i=0;i<iNphotons;i++){//feedbacks loop
 361     Double_t ranf[2];
 362     gMC->GetRandom()->RndmArray(2,ranf);    //Sample direction
 363     cthf=ranf[0]*2-1.0;
 364     if(cthf<0) continue;
 365     sthf = TMath::Sqrt((1 - cthf) * (1 + cthf));
 366     phif = ranf[1] * 2 * TMath::Pi();
 367
 368     if(Double_t randomNumber=gMC->GetRandom()->Rndm()<=0.57)
 369       enfp = 7.5e-9;
 370     else if(randomNumber<=0.7)
 371       enfp = 6.4e-9;
 372     else
 373       enfp = 7.9e-9;
 374
 375
 376     dir[0] = sthf * TMath::Sin(phif);    dir[1] = cthf;    dir[2] = sthf * TMath::Cos(phif);
 377     gMC->Gdtom(dir, mom, 2);
 378     mom[0]*=enfp;    mom[1]*=enfp;    mom[2]*=enfp;
 379     mom[3] = TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]+mom[2]*mom[2]);
 380
 381     // Polarisation
 382     e1[0]=      0;    e1[1]=-dir[2];    e1[2]= dir[1];
 383     e2[0]=-dir[1];    e2[1]= dir[0];    e2[2]=      0;
 384     e3[0]= dir[1];    e3[1]=      0;    e3[2]=-dir[0];
 385
 386     vmod=0;
 387     for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
 388     if (!vmod) for(j=0;j<3;j++) {
 389       uswop=e1[j];
 390       e1[j]=e3[j];
 391       e3[j]=uswop;
 392     }
 393     vmod=0;
 394     for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
 395     if (!vmod) for(j=0;j<3;j++) {
 396       uswop=e2[j];
 397       e2[j]=e3[j];
 398       e3[j]=uswop;
 399     }
 400
 401     vmod=0;  for(j=0;j<3;j++) vmod+=e1[j]*e1[j];  vmod=TMath::Sqrt(1/vmod);  for(j=0;j<3;j++) e1[j]*=vmod;
 402     vmod=0;  for(j=0;j<3;j++) vmod+=e2[j]*e2[j];  vmod=TMath::Sqrt(1/vmod);  for(j=0;j<3;j++) e2[j]*=vmod;
 403
 404     phi = gMC->GetRandom()->Rndm()* 2 * TMath::Pi();
 405     for(j=0;j<3;j++) pol[j]=e1[j]*TMath::Sin(phi)+e2[j]*TMath::Cos(phi);
 406     gMC->Gdtom(pol, pol, 2);
 407     Int_t outputNtracksStored;
 408     gAlice->GetMCApp()->PushTrack(1,                             //transport
 409                      gAlice->GetMCApp()->GetCurrentTrackNumber(),//parent track
 410                      kFeedback,                                  //PID
 411                      mom[0],mom[1],mom[2],mom[3],                //track momentum
 412                      x4.X(),x4.Y(),x4.Z(),x4.T(),                //track origin
 413                      pol[0],pol[1],pol[2],                       //polarization
 414                      kPFeedBackPhoton,                           //process ID
 415                      outputNtracksStored,                        //on return how many new photons stored on stack
 416                      1.0);                                       //weight
 417   }//feedbacks loop
 418   AliDebug(1,"Stop.");
 419 }//GenerateFeedbacks()
 420 //__________________________________________________________________________________________________
 421 void AliRICHv1::Hits2SDigits()
 422 {
 423 // Create a list of sdigits corresponding to list of hits. Every hit generates one or more sdigits.
 424   AliDebug(1,"Start.");
 425   for(Int_t iEventN=0;iEventN<GetLoader()->GetRunLoader()->GetAliRun()->GetEventsPerRun();iEventN++){//events loop
 426     GetLoader()->GetRunLoader()->GetEvent(iEventN);//get next event
 427
 428     if(!GetLoader()->TreeH()) GetLoader()->LoadHits();    GetLoader()->GetRunLoader()->LoadHeader();
 429     if(!GetLoader()->GetRunLoader()->TreeK())             GetLoader()->GetRunLoader()->LoadKinematics();//from
 430     if(!GetLoader()->TreeS()) GetLoader()->MakeTree("S"); MakeBranch("S");//to
 431
 432     for(Int_t iPrimN=0;iPrimN<GetLoader()->TreeH()->GetEntries();iPrimN++){//prims loop
 433       GetLoader()->TreeH()->GetEntry(iPrimN);
 434       for(Int_t iHitN=0;iHitN<Hits()->GetEntries();iHitN++){//hits loop
 435         AliRICHHit *pHit=(AliRICHHit*)Hits()->At(iHitN);//get current hit
 436         TVector2 x2 = AliRICHParam::Instance()->Mars2Lors(pHit->C(),0.5*(pHit->InX3()+pHit->OutX3()),AliRICHParam::kAnod);//hit position in the anod plane
 437         Int_t iTotQdc=AliRICHParam::TotQdc(x2,pHit->Eloss());//total charge produced by hit, 0 if hit in dead zone
 438         if(iTotQdc==0) continue;
 439         //
 440         //need to quantize the anod....
 441         TVector padHit=AliRICHParam::Loc2Pad(x2);
 442         TVector2 padHitXY=AliRICHParam::Pad2Loc(padHit);
 443         TVector2 anod;
 444         if((x2.Y()-padHitXY.Y())>0) anod.Set(x2.X(),padHitXY.Y()+AliRICHParam::AnodPitch()/2);
 445         else anod.Set(x2.X(),padHitXY.Y()-AliRICHParam::AnodPitch()/2);
 446         //end to quantize anod
 447         //
 448         TVector area=AliRICHParam::Loc2Area(anod);//determine affected pads, dead zones analysed inside
 449         AliDebug(1,Form("hitanod(%6.2f,%6.2f)->area(%3.0f,%3.0f)-(%3.0f,%3.0f) QDC=%4i",anod.X(),anod.Y(),area[0],area[1],area[2],area[3],iTotQdc));
 450         TVector pad(2);
 451         for(pad[1]=area[1];pad[1]<=area[3];pad[1]++)//affected pads loop
 452           for(pad[0]=area[0];pad[0]<=area[2];pad[0]++){
 453             Double_t padQdc=iTotQdc*AliRICHParam::FracQdc(anod,pad);
 454             AliDebug(1,Form("current pad(%3.0f,%3.0f) with QDC  =%6.2f",pad[0],pad[1],padQdc));
 455             if(padQdc>0.1) SDigAdd(pHit->C(),pad,padQdc,GetLoader()->GetRunLoader()->Stack()->Particle(pHit->GetTrack())->GetPdgCode(),pHit->GetTrack());
 456           }//affected pads loop
 457       }//hits loop
 458     }//prims loop
 459     GetLoader()->TreeS()->Fill();
 460     GetLoader()->WriteSDigits("OVERWRITE");
 461     SDigReset();
 462   }//events loop
 463   GetLoader()->UnloadHits(); GetLoader()->GetRunLoader()->UnloadHeader(); GetLoader()->GetRunLoader()->UnloadKinematics();
 464   GetLoader()->UnloadSDigits();
 465   AliDebug(1,"Stop.");
 466 }//Hits2SDigits()
 467 //__________________________________________________________________________________________________
 468 void AliRICHv1::Digits2Raw()
 469 {
 470 //Creates raw data files in DDL format. Invoked by AliSimulation
 471 //loop over events is done outside in AliSimulation
 472 //Arguments: none
 473 //  Returns: none
 474   AliDebug(1,"Start.");
 475   GetLoader()->LoadDigits();
 476   GetLoader()->TreeD()->GetEntry(0);
 477
 478   ofstream file[AliRICHDigit::kNddls];   //output streams
 479   Int_t    cnt[AliRICHDigit::kNddls];        //data words counters for DDLs
 480   AliRawDataHeader header; //empty DDL header
 481   UInt_t w32=0;            //32 bits data word
 482
 483   for(Int_t i=0;i<AliRICHDigit::kNddls;i++){//open all 14 DDL in parallel
 484     file[i].open(Form("RICH_%4i.ddl",AliRICHDigit::kDdlOffset+i));  //first is 0x700
 485     file[i].write((char*)&header,sizeof(header));                //write dummy header as place holder, actual will be written later when total size of DDL is known
 486     cnt[i]=0; //reset counters
 487   }
 488
 489   for(Int_t iCh=0;iCh<fNcham;iCh++){ //digits are stored on chamber by chamber basis
 490     TClonesArray *pDigs=(TClonesArray*)fDig->UncheckedAt(iCh);
 491     for(Int_t iDig=0;iDig<pDigs->GetEntriesFast();iDig++){//digits loop for a given chamber
 492       AliRICHDigit *pDig=(AliRICHDigit*)pDigs->At(iDig);
 493       Int_t ddl=pDig->Dig2Raw(w32);  //ddl is 0..13
 494       file[ddl].write((char*)&w32,sizeof(w32));  cnt[ddl]++;//write formated digit to the propriate file (as decided in Dig2Raw) and increment corresponding counter
 495     }//digits loop for a given chamber
 496   }//chambers loop
 497   for(Int_t i=0;i<AliRICHDigit::kNddls;i++){
 498     header.fSize=sizeof(header)+cnt[i]*sizeof(w32); //now calculate total number of bytes for each DDL file
 499     header.SetAttribute(0);
 500     file[i].seekp(0); file[i].write((char*)&header,sizeof(header));//rewrite DDL header with fSize field properly set
 501     file[i].close();                                               //close DDL file
 502   }
 503   GetLoader()->UnloadDigits();
 504   AliDebug(1,"Stop.");
 505 }//Digits2Raw()
 506 //__________________________________________________________________________________________________
 507 Float_t AliRICHv1::Fresnel(Float_t ene,Float_t pdoti, Bool_t pola)
 508 {
 509 // Correction for Fresnel   ???????????
 510 // Arguments:   ene - photon energy [GeV],
 511 //              PDOTI=COS(INC.ANG.), PDOTR=COS(POL.PLANE ROT.ANG.)
 512 //   Returns:
 513     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,
 514                       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,
 515                       7.8,7.9,8.0,8.1,8.2,8.3,8.4,8.5};
 516     Float_t csin[36] = {2.14,2.21,2.33,2.48,2.76,2.97,2.99,2.59,2.81,3.05,
 517                         2.86,2.53,2.55,2.66,2.79,2.96,3.18,3.05,2.84,2.81,2.38,2.11,
 518                         2.01,2.13,2.39,2.73,3.08,3.15,2.95,2.73,2.56,2.41,2.12,1.95,
 519                         1.72,1.53};
 520     Float_t csik[36] = {0.,0.,0.,0.,0.,0.196,0.408,0.208,0.118,0.49,0.784,0.543,
 521                         0.424,0.404,0.371,0.514,0.922,1.102,1.139,1.376,1.461,1.253,0.878,
 522                         0.69,0.612,0.649,0.824,1.347,1.571,1.678,1.763,1.857,1.824,1.824,
 523                         1.714,1.498};
 524     Float_t xe=ene;
 525     Int_t  j=Int_t(xe*10)-49;
 526     Float_t cn=csin[j]+((csin[j+1]-csin[j])/0.1)*(xe-en[j]);
 527     Float_t ck=csik[j]+((csik[j+1]-csik[j])/0.1)*(xe-en[j]);
 528
 529     //FORMULAE FROM HANDBOOK OF OPTICS, 33.23 OR
 530     //W.R. HUNTER, J.O.S.A. 54 (1964),15 , J.O.S.A. 55(1965),1197
 531
 532     Float_t sinin=TMath::Sqrt(1-pdoti*pdoti);
 533     Float_t tanin=sinin/pdoti;
 534
 535     Float_t c1=cn*cn-ck*ck-sinin*sinin;
 536     Float_t c2=4*cn*cn*ck*ck;
 537     Float_t aO=TMath::Sqrt(0.5*(TMath::Sqrt(c1*c1+c2)+c1));
 538     Float_t b2=0.5*(TMath::Sqrt(c1*c1+c2)-c1);
 539
 540     Float_t rs=((aO-pdoti)*(aO-pdoti)+b2)/((aO+pdoti)*(aO+pdoti)+b2);
 541     Float_t rp=rs*((aO-sinin*tanin)*(aO-sinin*tanin)+b2)/((aO+sinin*tanin)*(aO+sinin*tanin)+b2);
 542
 543
 544     //CORRECTION FACTOR FOR SURFACE ROUGHNESS
 545     //B.J. STAGG  APPLIED OPTICS, 30(1991),4113
 546
 547     Float_t sigraf=18.;
 548     Float_t lamb=1240/ene;
 549     Float_t fresn;
 550
 551     Float_t  rO=TMath::Exp(-(4*TMath::Pi()*pdoti*sigraf/lamb)*(4*TMath::Pi()*pdoti*sigraf/lamb));
 552
 553     if(pola)
 554     {
 555         Float_t pdotr=0.8;                                 //DEGREE OF POLARIZATION : 1->P , -1->S
 556         fresn=0.5*(rp*(1+pdotr)+rs*(1-pdotr));
 557     }
 558     else
 559         fresn=0.5*(rp+rs);
 560
 561     fresn = fresn*rO;
 562     return fresn;
 563 }//Fresnel()
 564 //__________________________________________________________________________________________________
 565 void AliRICHv1::Print(Option_t *option)const
 566 {
 567 // Debug printout
 568   TObject::Print(option);
 569   Printf("Total number of MIP reached radiator        %9i",fCounters(kMipEnterRad));
 570   Printf("Total number of Ckov created                %9i",fCounters(kCkovNew));
 571   Printf("number of Ckov created in radiator          %9i",fCounters(kCkovNewRad));
 572   Printf("number of Ckov created in window            %9i",fCounters(kCkovNewWin));
 573   Printf("number of Ckov created in proximity gap     %9i",fCounters(kCkovNewProxGap));
 574   Printf("number of Ckov created in amplification gap %9i",fCounters(kCkovNewAmpGap));
 575   Printf("number of Ckov reached PC                   %9i",fCounters(kCkovEnterPc));
 576   Printf("number of photelectrons                     %9i",fCounters(kPhotoEle));
 577 }//void AliRICH::Print(Option_t *option)const
 578 //__________________________________________________________________________________________________
 579 void AliRICHv1::StepCount()
 580 {
 581 // Count number of ckovs created
 582   Int_t copy;
 583   if(gMC->TrackCharge()        &&gMC->CurrentVolID(copy)==fIdRad    &&gMC->IsTrackEntering()                ) fCounters(kMipEnterRad)++;
 584   if(gMC->TrackPid()==kCerenkov                                     &&gMC->IsNewTrack()                     ) fCounters(kCkovNew)++;
 585   if(gMC->TrackPid()==kCerenkov&&gMC->CurrentVolID(copy)==fIdRad    &&gMC->IsNewTrack()                     ) fCounters(kCkovNewRad)++;
 586   if(gMC->TrackPid()==kCerenkov&&gMC->CurrentVolID(copy)==fIdWin    &&gMC->IsNewTrack()                     ) fCounters(kCkovNewWin)++;
 587   if(gMC->TrackPid()==kCerenkov&&gMC->CurrentVolID(copy)==fIdProxGap&&gMC->IsNewTrack()                     ) fCounters(kCkovNewProxGap)++;
 588   if(gMC->TrackPid()==kCerenkov&&gMC->CurrentVolID(copy)==fIdAmpGap &&gMC->IsNewTrack()                     ) fCounters(kCkovNewAmpGap)++;
 589   if(gMC->TrackPid()==kCerenkov&&gMC->CurrentVolID(copy)==fIdPc     &&gMC->IsTrackEntering()                ) fCounters(kCkovEnterPc)++;
 590   if(gMC->TrackPid()==kCerenkov&&gMC->CurrentVolID(copy)==fIdPc     &&gMC->IsTrackEntering() &&gMC->Edep()>0) fCounters(kPhotoEle)++;
 591 }
 592 //__________________________________________________________________________________________________
 593 void AliRICHv1::StepHistory()
 594 {
 595 // This methode is invoked from StepManager() in order to print out
 596   static Int_t iStepN;
 597   const char *sParticle;
 598   switch(gMC->TrackPid()){
 599     case kProton:      sParticle="PROTON"    ;break;
 600     case kNeutron:     sParticle="neutron"   ;break;
 601     case kGamma:       sParticle="gamma"     ;break;
 602     case kCerenkov:    sParticle="CKOV"    ;break;
 603     case kPi0:         sParticle="Pi0"       ;break;
 604     case kPiPlus:      sParticle="Pi+"       ;break;
 605     case kPiMinus:     sParticle="Pi-"       ;break;
 606     case kElectron:    sParticle="electron"  ;break;
 607     default:           sParticle="not known" ;break;
 608   }
 609
 610   TString flag="fanny combination";
 611   if(gMC->IsTrackAlive())
 612       if(gMC->IsTrackEntering())      flag="enters to";
 613       else if(gMC->IsTrackExiting())  flag="exits from";
 614       else if(gMC->IsTrackInside())   flag="inside";
 615   else
 616       if(gMC->IsTrackStop())          flag="stoped in";
 617
 618   Int_t vid=0,copy=0;
 619   TString path=gMC->CurrentVolName(); path.Prepend("-");path.Prepend(gMC->CurrentVolOffName(1));//current volume and his mother are always there
 620   vid=gMC->CurrentVolOffID(2,copy);  if(vid) {path.Prepend("-");path.Prepend(gMC->VolName(vid));}
 621   vid=gMC->CurrentVolOffID(3,copy);  if(vid) {path.Prepend("-");path.Prepend(gMC->VolName(vid));}
 622
 623   Printf("Step %i: %s (%i) %s %s m=%.6f GeV q=%.1f dEdX=%.4f",iStepN,sParticle,gMC->TrackPid(),flag.Data(),path.Data(),gMC->TrackMass(),gMC->TrackCharge(),gMC->Edep()*1e9);
 624
 625   Printf("Step %i: tid=%i flags alive=%i disap=%i enter=%i exit=%i inside=%i out=%i stop=%i new=%i",
 626                             iStepN, gAlice->GetMCApp()->GetCurrentTrackNumber(),
 627                             gMC->IsTrackAlive(), gMC->IsTrackDisappeared(),gMC->IsTrackEntering(), gMC->IsTrackExiting(),
 628                             gMC->IsTrackInside(),gMC->IsTrackOut(),        gMC->IsTrackStop(),     gMC->IsNewTrack());
 629
 630   Float_t a,z,den,rad,abs; a=z=den=rad=abs=-1;
 631   Int_t mid=gMC->CurrentMaterial(a,z,den,rad,abs);
 632   Printf("Step %i: id=%i a=%7.2f z=%7.2f den=%9.4f rad=%9.2f abs=%9.2f\n\n",iStepN,mid,a,z,den,rad,abs);
 633   iStepN++;
 634 }//StepHistory()
 635 //__________________________________________________________________________________________________
 636 void AliRICHv1::StepManager()
 637 {
 638 // Full Step Manager.
 639 // Arguments: none
 640 //   Returns: none
 641 //  StepHistory(); return; //uncomment to print tracks history
 642 //  StepCount(); return;     //uncomment to count photons
 643
 644   Int_t          copy; //volume copy aka node
 645   static Int_t   iCham;
 646 //Treat photons
 647   static TLorentzVector cerX4;
 648   if((gMC->TrackPid()==kCerenkov||gMC->TrackPid()==kFeedback)&&gMC->CurrentVolID(copy)==fIdPc){//photon in PC
 649     if(gMC->Edep()>0){//photon survided QE test i.e. produces electron
 650       if(IsLostByFresnel()){ gMC->StopTrack(); return;} //photon lost due to fersnel reflection
 651       gMC->TrackPosition(cerX4); gMC->CurrentVolOffID(2,iCham);//RICH-Rppf-Rpc
 652       HitAdd(iCham,gMC->GetStack()->GetCurrentTrackNumber(),gMC->TrackPid(),cerX4.Vect(),cerX4.Vect());//HIT for PHOTON in conditions CF+CSI+DE
 653       GenFee(iCham);
 654     }//photon in PC and DE >0
 655   }//photon in PC
 656
 657 //Treat charged particles
 658   static Float_t eloss;
 659   static TLorentzVector mipInX4,mipOutX4;
 660   if(gMC->TrackCharge() && gMC->CurrentVolID(copy)==fIdAmpGap){//MIP in amplification gap
 661     gMC->CurrentVolOffID(1,iCham);//RICH-Rgap
 662     if(gMC->IsTrackEntering()||gMC->IsNewTrack()) {//MIP in GAP entering or newly created
 663       eloss=0;
 664       gMC->TrackPosition(mipInX4);
 665     }else if(gMC->IsTrackExiting()||gMC->IsTrackStop()||gMC->IsTrackDisappeared()){//MIP in GAP exiting or disappeared
 666       eloss+=gMC->Edep();//take into account last step dEdX
 667       gMC->TrackPosition(mipOutX4);
 668       HitAdd(iCham,gMC->GetStack()->GetCurrentTrackNumber(),gMC->TrackPid(),mipInX4.Vect(),mipOutX4.Vect(),eloss);//HIT for MIP: MIP in GAP Exiting
 669       GenFee(iCham,eloss);//MIP+GAP+Exit
 670     }else//MIP in GAP going inside
 671       eloss   += gMC->Edep();
 672   }//MIP in GAP
 673 }//StepManager()
 674 //__________________________________________________________________________________________________