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