]> git.uio.no Git - u/mrichter/AliRoot.git/blame - RICH/AliRICHv1.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Check the Minuit status after each fit
[u/mrichter/AliRoot.git] / RICH / AliRICHv1.cxx
CommitLineData
53fd478b 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// **************************************************************************
c28632f0 15
c28632f0 16
0422a446 17#include "AliRICHv1.h" //class header
53fd478b 18#include "AliRICHParam.h"
db910db9 19#include <TParticle.h> //Hits2SDigits()
88cb7938 20#include <TRandom.h>
db910db9 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()
c60862bf 26#include <AliConst.h>
27#include <AliPDG.h>
db910db9 28#include <AliMC.h> //StepManager()
0422a446 29#include <AliRawDataHeader.h> //Digits2Raw()
db910db9 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
d128c9d1 41ClassImp(AliRICHv1)
db910db9 42
43void 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()
3582c1f9 203//__________________________________________________________________________________________________
db910db9 204void AliRICHv1::CreateGeometry()
53fd478b 205{
db910db9 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);
e42a7b46 234 }
db910db9 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)
ed3ceb24 243
db910db9 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//__________________________________________________________________________________________________
312void 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}
3582c1f9 325//__________________________________________________________________________________________________
326Bool_t AliRICHv1::IsLostByFresnel()
ed3ceb24 327{
e42a7b46 328// Calculate probability for the photon to be lost by Fresnel reflection.
3582c1f9 329 TLorentzVector p4;
330 Double_t mom[3],localMom[3];
ab296411 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;
3582c1f9 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)){
998b831f 338 AliDebug(1,"Photon lost");
3582c1f9 339 return kTRUE;
340 }else
341 return kFALSE;
342}//IsLostByFresnel()
343//__________________________________________________________________________________________________
db910db9 344void AliRICHv1::GenFee(Int_t iChamber,Float_t eloss)
e42a7b46 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);
db910db9 351 TVector2 x2=AliRICHParam::Instance()->Mars2Lors(iChamber,x4.Vect(),AliRICHParam::kPc);//hit position on photocathode plane
998b831f 352 TVector2 xspe=x2;
db910db9 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);
998b831f 356 AliDebug(1,Form("N photons=%i",iNphotons));
e42a7b46 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;
db910db9 408 gAlice->GetMCApp()->PushTrack(1, //transport
e42a7b46 409 gAlice->GetMCApp()->GetCurrentTrackNumber(),//parent track
db910db9 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
e42a7b46 417 }//feedbacks loop
998b831f 418 AliDebug(1,"Stop.");
e42a7b46 419}//GenerateFeedbacks()
0422a446 420//__________________________________________________________________________________________________
421void 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
db910db9 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
0422a446 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;
db910db9 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);
0422a446 446 //end to quantize anod
447 //
db910db9 448 TVector area=AliRICHParam::Loc2Area(anod);//determine affected pads, dead zones analysed inside
0422a446 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]++){
db910db9 453 Double_t padQdc=iTotQdc*AliRICHParam::FracQdc(anod,pad);
0422a446 454 AliDebug(1,Form("current pad(%3.0f,%3.0f) with QDC =%6.2f",pad[0],pad[1],padQdc));
db910db9 455 if(padQdc>0.1) SDigAdd(pHit->C(),pad,padQdc,GetLoader()->GetRunLoader()->Stack()->Particle(pHit->GetTrack())->GetPdgCode(),pHit->GetTrack());
0422a446 456 }//affected pads loop
457 }//hits loop
458 }//prims loop
459 GetLoader()->TreeS()->Fill();
460 GetLoader()->WriteSDigits("OVERWRITE");
db910db9 461 SDigReset();
0422a446 462 }//events loop
463 GetLoader()->UnloadHits(); GetLoader()->GetRunLoader()->UnloadHeader(); GetLoader()->GetRunLoader()->UnloadKinematics();
464 GetLoader()->UnloadSDigits();
465 AliDebug(1,"Stop.");
466}//Hits2SDigits()
467//__________________________________________________________________________________________________
468void 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
db910db9 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);
0422a446 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//__________________________________________________________________________________________________
db910db9 507Float_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//__________________________________________________________________________________________________
565void 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//__________________________________________________________________________________________________
579void 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//__________________________________________________________________________________________________
593void 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//__________________________________________________________________________________________________
636void 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//__________________________________________________________________________________________________