Raw data reconstruction (K.Shileev)
[u/mrichter/AliRoot.git] / RICH / AliRICHParam.h
1 #ifndef AliRICHParam_h
2 #define AliRICHParam_h
3
4 #include <TError.h>
5 #include <TMath.h>
6 #include <TObjArray.h>
7 #include <TObject.h>
8 #include <TMath.h>
9 #include <TRandom.h>
10 #include <TVector.h>
11 #include <TVector2.h>
12 #include <TVector3.h>
13 #include <TRandom.h>
14 #include <TError.h>
15 #include <TObjArray.h>
16 #include <AliLog.h>
17 #include <TClass.h>
18
19 static const int kNchambers=7;     //number of RICH chambers 
20 static const int kNpadsX = 160;    //number of pads along X in single chamber
21 static const int kNpadsY = 144;    //number of pads along Y in single chamber
22 static const int kNsectors=6;      //number of sectors per chamber
23
24 static const int kCerenkov=50000050;  //??? go to something more general like TPDGCode
25 static const int kFeedback=50000051;  //??? go to something more general like TPDGCode
26
27 class AliRICHChamber;
28
29 // Class providing all the needed parametrised information
30 // to construct the geometry, to define segmentation and to provide response model
31 // In future will also provide all the staff needed for alignment and calibration
32
33
34 class AliRICHParam :public TObject  
35 {
36 public:
37 //ctor&dtor    
38                   AliRICHParam():TObject(),fpChambers(0)  {CreateChambers();}
39   virtual        ~AliRICHParam()                          {delete fpChambers;}
40 //test methodes  
41          void     Print(Option_t *opt="") const;                                         //print current parametrization
42          void     Test()                            {TestSeg();TestTrans();TestResp();}  //test all groups of methodes
43          void     TestResp();                                                            //test the response group of methodes
44          void     TestSeg();                                                             //test the segmentation group of methodes
45          void     TestTrans();                                                           //test the transform group of methodes
46   static void     DrawAxis();
47   static void     DrawSectors();
48 //flags staff         
49   static           void     SetAerogel(Bool_t a)                   {fgIsAerogel=a;}
50   static           Bool_t   IsAerogel()                            {return fgIsAerogel;}
51   static           void     SetRadioSrc(Bool_t a)                   {fgIsRadioSrc=a;}
52   static           Bool_t   IsRadioSrc()                            {return fgIsRadioSrc;}
53   static              void     SetTestBeam(Bool_t a)                {fgIsTestBeam=a;}
54   static              Bool_t   IsTestBeam()                         {return fgIsTestBeam;}
55   static                void     SetWireSag(Bool_t a)               {fgIsWireSag=a;}
56   static                Bool_t   IsWireSag()                        {return fgIsWireSag;}
57   static                   void     SetResolveClusters(Bool_t a)    {fgIsResolveClusters=a;}
58   static                   Bool_t   IsResolveClusters()             {return fgIsResolveClusters;}
59 //Chambers manipulation  methodes 
60   void            CreateChambers();                                                      //form chamber structure  
61   AliRICHChamber* C(Int_t i)                 {return (AliRICHChamber*)fpChambers->UncheckedAt(i-1);}      //returns pointer to chamber i
62   Int_t           Nchambers()                {return fpChambers->GetEntriesFast();}      //returns number of chambers 
63 //Geometrical properties  
64   static        Int_t      NpadsX()                   {return kNpadsX;}                           //pads along X in chamber
65   static        Int_t      NpadsY()                   {return kNpadsY;}                           //pads along Y in chamber
66   static        Int_t      NpadsXsec()                {return NpadsX()/2;}                        //pads along X in sector
67   static        Int_t      NpadsYsec()                {return NpadsY()/3;}                        //pads along Y in sector
68   static        Double_t   DeadZone()                 {return 2.6;}                               //dead zone size in cm  
69   static        Double_t   SectorSizeX()              {return NpadsX()*PadSizeX()/2;}             //sector size x, cm
70   static        Double_t   SectorSizeY()              {return NpadsY()*PadSizeY()/3;}             //sector size y, cm 
71   static        Double_t   PcSizeX()                  {return NpadsX()*PadSizeX()+DeadZone();}    //PC size x, cm
72   static        Double_t   PcSizeY()                  {return NpadsY()*PadSizeY()+2*DeadZone();}  //PC size y, cm
73   static        Double_t   Zfreon()                   {return 1.5;}                               //freon thinkness, cm
74   static        Double_t   Zwin()                     {return 0.5;}                               //radiator quartz window, cm   
75   static        Double_t   Pc2Win()                   {return 8.0;}                               //cm between CsI PC and radiator quartz window
76   static        Double_t   Pc2Coll()                  {return 7.0;}                               //cm between CsI PC and third wire grid (collection wires)     
77   static        Double_t   Pc2Anod()                  {return 0.204;}                             //cm between CsI PC and first wire grid (anod wires)     
78   static        Double_t   Pc2Cath()                  {return 0.445;}                             //cm between CsI PC and second wire grid (cathode wires)
79   static        Double_t   Freon2Pc()                 {return Zfreon()+Zwin()+Pc2Win();}          //cm between CsI PC and entrance to freon
80   static        Double_t   PitchAnod()                {return PadSizeY()/2;}                      //cm between anode wires
81   static        Double_t   PitchCath()                {return PadSizeY()/4;}                      //cm between cathode wires
82   static        Double_t   PitchColl()                {return 0.5;}                               //cm between collection wires
83   static        Double_t   PadSizeX()                 {return 0.8;}                               //pad size x, cm 
84   static        Double_t   PadSizeY      (                               ){return 0.84;}                              //pad size y, cm   
85 //trasformation methodes
86   static        Int_t      Pad2Cha       (Int_t pad                      ){return pad/100000000;                     }//abs pad -> chamber
87   static        Int_t      Pad2Sec       (Int_t pad                      ){return pad%100000000/1000000;             }//abs pad -> sector
88   static        Int_t      Pad2PadX      (Int_t pad                      ){return pad%1000000/1000;                  }//abs pad -> pad x 
89   static        Int_t      Pad2PadY      (Int_t pad                      ){return pad%1000000%100;                   }//abs pad -> pad y
90   static        Int_t      PadAbs        (Int_t c,Int_t s,Int_t x,Int_t y){return 100000000*c+1000000*s+1000*x+y;    }//(c,s,x,y) -> abs pad
91   static inline TVector2   Pad2Loc       (Int_t pad                      );                                           //abs pad ->LORS
92   static inline TVector2   Pad2Loc       (TVector pad                    );                                           //pad  -> LORS returns center of the pad
93   static        TVector2   Pad2Loc       (Int_t x,Int_t y                ){TVector pad(2);pad[0]=x;pad[1]=y;return Pad2Loc(pad);}//return center of the pad (x,y)
94   static inline TVector    Loc2Area      (const TVector2 &x2             );                                           //pads area affected by hit x2. area is LeftDown-RightUp pad numbers
95   static inline Int_t      Loc2Sec       (const TVector2 &x2             );                                           //LORS -> sector
96   static        Int_t      Loc2Sec       (Double_t x,Double_t y          ){return Loc2Sec(TVector2(x,y));}            //LORS -> sector
97   static inline TVector    Loc2Pad       (const TVector2 &x2             );                                           //LORS -> pad
98   static        TVector    Loc2Pad       (Double_t x,Double_t y          ){return Loc2Pad(TVector2(x,y));}            //LORS -> pad
99   static inline Int_t      Pad2Sec       (const TVector &pad             );                                           //pad  -> sector
100   static inline Int_t      PadNeighbours (Int_t iPadX,Int_t iPadY,Int_t aListX[4],Int_t aListY[4]);                   //pad -> list of it neighbours
101   static        Bool_t     IsAccepted    (const TVector2 &x2             ){return ( x2.X()>=0 && x2.X()<=PcSizeX() && x2.Y()>=0 && x2.Y()<=PcSizeY() ) ? kTRUE:kFALSE;}
102 //optical properties methodes  
103   static        Double_t   MeanCkovEnergy(                               ){return 6.766;}                 //mean Ckov energy according to the total trasmission curve
104   static        Float_t    PhotonEnergy  (Int_t i                        ){return 0.1*i+5.5;}                         //photon energy (eV) for i-th point
105   static        Float_t    AbsCH4        (Float_t ev                     );                                           //CH4 abs len (cm) 
106   static        Float_t    AbsGel        (Float_t                        ){return 500;}                               //Aerogel abs len (cm)
107   static        Float_t    RefIdxC6F14   (Float_t eV                     ){return eV*0.0172+1.177;}                   //Freon ref idx
108   static        Float_t    RefIdxCH4     (Float_t                        ){return 1.000444;}                          //Methane ref idx 
109   static        Float_t    RefIdxSiO2    (Float_t eV                     ){Float_t e1=10.666,e2=18.125,f1=46.411,f2= 228.71; return TMath::Sqrt(1.+f1/(e1*e1-eV*eV)+f2/(e2*e2-eV*eV));}//Quartz window ref index from TDR p.35
110   static        Float_t    RefIdxGel     (Float_t                        ){return 1.05;}                              //aerogel ref index 
111   static        Float_t    DenGel        (                               ){return (RefIdxGel(0)-1)/0.21;}             //aerogel density gr/cm^3 parametrization by E.Nappi
112
113   
114   static Double_t IonisationPotential()      {return 26.0e-9;}                            //for CH4 in GeV taken from ????
115   static TVector2 MathiesonDelta()           {return TVector2(5*0.18,5*0.18);}            //area of 5 sigmas of Mathieson distribution (cm)
116   static Int_t    MaxQdc()                   {return 4095;}                               //QDC number of channels          
117   
118   static Int_t    HV(Int_t sector)           {if (sector>=1 && sector <=6) return fgHV[sector-1];  else return -1;} //high voltage for this sector
119   static void     SetHV(Int_t sector,Int_t hv){fgHV[sector-1]=hv;}  
120 //charge response methodes  
121   inline static Double_t Mathieson(Double_t x1,Double_t x2,Double_t y1,Double_t y2);               //Mathienson integral over given limits
122   inline static Double_t GainSag(Double_t x,Int_t sector);                                         //gain variations in %
123          static Double_t QdcSlope(Int_t sec){switch(sec){case -1: return 0;  default:   return 33;}} //weight of electon in QDC channels
124          static Double_t Gain(const TVector2 &x2){//gives chamber gain in terms of QDC channels for given point in local ref system
125                           if(fgIsWireSag) return QdcSlope(Loc2Sec(x2))*(1+GainSag(x2.X(),Loc2Sec(x2))/100);
126                           else            return QdcSlope(Loc2Sec(x2));}
127   inline static Double_t FracQdc(const TVector2 &x2,const TVector &pad);                           //charge fraction to pad from hit
128   inline static Int_t    TotQdc(TVector2 x2,Double_t eloss);                                       //total charge for Eloss (GeV) 0 for photons
129   inline static Bool_t   IsOverTh(Int_t c,TVector pad,Double_t q);                                 //is QDC of the pad registered by FEE  
130          static Int_t    NsigmaTh()                    {return fgNsigmaTh;}                        //
131          static Float_t  SigmaThMean()                 {return fgSigmaThMean;}                     //QDC electronic noise mean
132          static Float_t  SigmaThSpread()               {return fgSigmaThSpread;}                   //QDC electronic noise width
133                 
134          static Double_t CogCorr(Double_t x) {return 3.31267e-2*TMath::Sin(2*TMath::Pi()/PadSizeX()*x) //correction of cluster CoG due to sinoidal
135                                                     -2.66575e-3*TMath::Sin(4*TMath::Pi()/PadSizeX()*x)
136                                                     +2.80553e-3*TMath::Sin(6*TMath::Pi()/PadSizeX()*x)+0.0070;}
137          static void     ReadErrFiles();                                                                  //Read Err file parameters
138          static TVector3 SigmaSinglePhoton(Int_t Npart, Double_t mom, Double_t theta, Double_t phi);      //Find Sigma for single photon
139          static Double_t Interpolate(Double_t par[4][330],Double_t x, Double_t y, Double_t phi);          //Find the error value from interpolation
140          
141          static TVector3 ForwardTracing(TVector3 entranceTrackPoint,TVector3 vectorTrack, Double_t thetaC, Double_t phiC); //it traces foward a photon from Emission Point to PC
142          static TVector3 PlaneIntersect(TVector3 vstart,TVector3 p0,TVector3 n,TVector3 v0);              //it finds intersection between straight track and plane
143          static Double_t SnellAngle(Float_t n1, Float_t n2, Float_t theta1);                              // Snell law
144          static void     AnglesInDRS(Double_t trackTheta,Double_t trackPhi,Double_t thetaCerenkov,Double_t phiCerenkov,Double_t &tout,Double_t &pout);//It finds photon angles in 
145                                                                                                                                                       //Detector Reference System
146
147   static Bool_t     fgIsAerogel;                            //aerogel geometry instead of normal RICH flag
148 protected:
149   static Bool_t     fgIsRadioSrc;                           //radioactive source instead of radiators flag
150   static Bool_t     fgIsTestBeam;                           //test beam geometry instead of normal RICH flag
151   static Bool_t     fgIsWireSag;                            //wire sagitta ON/OFF flag
152   static Bool_t     fgIsResolveClusters;                    //declustering ON/OFF flag
153   static Bool_t     fgIsFeedback;                           //generate feedback photon?
154
155          TObjArray *fpChambers;                             //list of chambers    
156   static Int_t      fgHV[6];                                //HV applied to anod wires
157   static Int_t      fgNsigmaTh;                             //n. of sigmas to cut for zero suppression
158   static Float_t    fgSigmaThMean;                          //sigma threshold value
159   static Float_t    fgSigmaThSpread;                        //spread of sigma
160   
161   static Double_t fgErrChrom[4][330];                       //
162   static Double_t fgErrGeom[4][330];                        //
163   static Double_t fgErrLoc[4][330];                         //Chromatic, Geometric and Localization array to parametrize SigmaCerenkov
164   
165   ClassDef(AliRICHParam,6)                                  //RICH main parameters class
166 };
167 //__________________________________________________________________________________________________
168 Int_t AliRICHParam::PadNeighbours(Int_t iPadX,Int_t iPadY,Int_t listX[4],Int_t listY[4])
169 {
170 //Determines all the neighbouring pads for the given one (iPadX,iPadY). Returns total number of these pads.
171 //Dead zones are taken into account, meaning pads from different sector are not taken. 
172 //   1  
173 // 2   3
174 //   4     
175   Int_t nPads=0;
176   if(iPadY!=NpadsY()&&iPadY!=2*NpadsYsec()&&iPadY!=NpadsYsec()){listX[nPads]=iPadX;   listY[nPads]=iPadY+1; nPads++;}       //1
177   if(iPadX!=1&&iPadX!=NpadsXsec()+1)                           {listX[nPads]=iPadX-1; listY[nPads]=iPadY;   nPads++;}       //2
178   if(iPadX!=NpadsXsec()&&iPadX!=NpadsX())                      {listX[nPads]=iPadX+1; listY[nPads]=iPadY;   nPads++;}       //3
179   if(iPadY!=1&&iPadY!=NpadsYsec()+1&&2*NpadsYsec()+1)          {listX[nPads]=iPadX;   listY[nPads]=iPadY-1; nPads++;}       //4
180
181   return nPads;
182 }//Pad2ClosePads()
183 //__________________________________________________________________________________________________
184 Int_t AliRICHParam::Loc2Sec(const TVector2 &v2)
185 {
186 //Determines sector containing the given point.
187 //Returns sector code:                       
188 //y ^  5 6
189 //  |  3 4
190 //  |  1 2
191 //   -------> x  
192   Double_t x0=0; Double_t x1=SectorSizeX(); Double_t x2=SectorSizeX()+DeadZone(); Double_t x3=PcSizeX();
193   Double_t y0=0; Double_t y1=SectorSizeY(); Double_t y2=SectorSizeY()+DeadZone(); Double_t y3=2*SectorSizeY()+DeadZone(); 
194   Double_t y4=PcSizeY()-SectorSizeY();      Double_t y5=PcSizeY();
195   
196   Int_t sector=-1;  
197   if     (v2.X() >= x0 && v2.X() <= x1 )  sector=1;
198   else if(v2.X() >= x2 && v2.X() <= x3 )  sector=2;
199   else                                    return -1;
200   
201   if     (v2.Y() >= y0 && v2.Y() <= y1 )  ;                    //sectors 1 or 2 
202   else if(v2.Y() >= y2 && v2.Y() <= y3 )  sector+=2;           //sectors 3 or 4
203   else if(v2.Y() >= y4 && v2.Y() <= y5 )  sector+=4;           //sectors 5 or 6
204   else                                    return -1;
205   return sector;
206 }//Loc2Sec(Double_t x, Double_t y)
207 //__________________________________________________________________________________________________
208 TVector AliRICHParam::Loc2Pad(const TVector2 &loc)
209 {
210 //Determines pad number TVector(padx,pady) containing the given point x2 defined in the chamber RS.
211 //Pad count starts in lower left corner from 1,1 to 144,160 in upper right corner of a chamber.
212 //y ^  5 6
213 //  |  3 4
214 //  |  1 2
215 //   -------> x  
216   TVector pad(2);
217   Int_t sec=Loc2Sec(loc);//trasforms x2 to sector reference system
218   if(sec==-1) {pad[0]=pad[1]=-1; return pad;}
219 //first we deal with x  
220   if(sec==1||sec==3||sec==5)    pad[0]=           Int_t(            loc.X()   / PadSizeX() )+1; //sector 1 or 3 or 5
221   else                          pad[0]=NpadsX() - Int_t( (PcSizeX()-loc.X())  / PadSizeX() )  ; //sector 2 or 4 or 6
222 //second deal with y
223        if(sec==1||sec==2)       pad[1]=Int_t(             loc.Y()                / PadSizeY())+1;               //sector 1 or 2 
224   else if(sec==3||sec==4)       pad[1]=Int_t( (loc.Y()-SectorSizeY()-DeadZone()) / PadSizeY())+NpadsYsec()+1;  //sector 3 or 4    
225   else                          pad[1]=NpadsY() - Int_t( (PcSizeY()-loc.Y())     / PadSizeY());                //sector 5 or 6        
226   return pad;
227 }
228 //__________________________________________________________________________________________________
229 Int_t AliRICHParam::Pad2Sec(const TVector &pad)
230 {
231 //Determines sector containing the given pad.
232   Int_t sector=-1;      
233   if     (pad[0] >= 1           && pad[0] <=   NpadsXsec() )    {sector=1;}
234   else if(pad[0] >  NpadsXsec() && pad[0] <=   NpadsX()    )    {sector=2;} 
235   else                                                         AliDebugClass(1,Form("Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]));
236     
237   if     (pad[1] >= 1             && pad[1] <=   NpadsYsec() )    {}
238   else if(pad[1] >  NpadsYsec()   && pad[1] <= 2*NpadsYsec() )    {sector+=2;}
239   else if(pad[1] >  2*NpadsYsec() && pad[1] <=   NpadsY()    )    {sector+=4;}
240   else                                                         AliDebugClass(1,Form("Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]));
241
242   return sector;
243 }//Pad2Sec()
244 //__________________________________________________________________________________________________
245 TVector2 AliRICHParam::Pad2Loc(TVector pad)
246 {
247 //Returns position of the center of the given pad in local system of the chamber (cm)    
248 // y ^  5 6
249 //   |  3 4        sector numbers
250 //   |  1 2
251 //    -------> x  
252   Double_t x=-1,y=-1;
253   if(pad[0] > 0 && pad[0] <= NpadsXsec())//it's 1 or 3 or 5
254     x=(pad[0]-0.5)*PadSizeX();
255   else if(pad[0] > NpadsXsec() && pad[0] <= NpadsX())//it's 2 or 4 or 6
256     x=(pad[0]-0.5)*PadSizeX()+DeadZone();
257   else
258     AliDebugClass(1,Form("Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]));
259   
260   if(pad[1] > 0 && pad[1] <= NpadsYsec())//it's 1 or 2
261     y=(pad[1]-0.5)*PadSizeY();
262   else if(pad[1] > NpadsYsec() && pad[1] <= 2*NpadsYsec())//it's 3 or 4
263     y=(pad[1]-0.5)*PadSizeY()+DeadZone();
264   else if(pad[1] > 2*NpadsYsec() && pad[1]<= NpadsY())//it's 5 or 6
265     y=(pad[1]-0.5)*PadSizeY()+2*DeadZone();
266   else
267     AliDebugClass(1,Form("Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]));
268     
269   return TVector2(x,y);
270 }
271 //__________________________________________________________________________________________________
272 TVector2 AliRICHParam::Pad2Loc(Int_t pad)
273 {
274 //Converts absolute pad number to local position in LORS
275 //LORS is a chamber  reference system with origin in left-down coner looking from IP
276 //Arguments: pad- absolute pad number
277 //  Returns: pad center position as TVector2 in PCRS  
278   TVector2 pos;
279   pos.Set((Pad2PadX(pad)-0.5)*PadSizeX() , (Pad2PadY(pad)-0.5)*PadSizeY());//set to sector LORS
280   return pos;
281 }
282 //__________________________________________________________________________________________________
283 Double_t AliRICHParam::GainSag(Double_t x,Int_t sector)
284 {
285 //Returns % of gain variation due to wire sagita.
286 //All curves are parametrized as per sector basis, so x must be apriory transformed to the Sector RS.    
287 //Here x is a distance along wires.  
288   x-=SectorSizeX()/2;
289   if(x>SectorSizeX()) x-=SectorSizeX(); 
290   switch(HV(sector)){
291     case 2150: return 9e-6*TMath::Power(x,4)+2e-7*TMath::Power(x,3)-0.0316*TMath::Power(x,2)-3e-4*x+25.367;//%
292     case 2100: return 8e-6*TMath::Power(x,4)+2e-7*TMath::Power(x,3)-0.0283*TMath::Power(x,2)-2e-4*x+23.015;
293     case 2050: return 7e-6*TMath::Power(x,4)+1e-7*TMath::Power(x,3)-0.0254*TMath::Power(x,2)-2e-4*x+20.888;
294     case 2000: return 6e-6*TMath::Power(x,4)+8e-8*TMath::Power(x,3)-0.0227*TMath::Power(x,2)-1e-4*x+18.961;
295     default:   return 0;
296   }
297 }
298 //__________________________________________________________________________________________________
299 Int_t AliRICHParam::TotQdc(TVector2 x2,Double_t eloss)
300 {
301 //Calculates the total charge produced by the eloss in point x2 (Chamber RS).
302 //Returns this change parametrised in QDC channels, or 0 if the hit in the dead zone.
303 //eloss=0 means photon which produces 1 electron only eloss > 0 for Mip
304   if(Loc2Sec(x2)==-1) return 0; //hit in the dead zone     
305   Int_t iNelectrons=Int_t(eloss/IonisationPotential()); if(iNelectrons==0) iNelectrons=1;
306   Double_t qdc=0;
307   for(Int_t i=1;i<=iNelectrons;i++) qdc+=-Gain(x2)*TMath::Log(gRandom->Rndm());
308   return Int_t(qdc);
309 }
310 //__________________________________________________________________________________________________
311 Double_t AliRICHParam::FracQdc(const TVector2 &x2,const TVector &pad)
312 {
313 //Calculates the charge fraction induced to given pad by the hit from the given point.
314 //Integrated Mathieson distribution is used.  
315   TVector2 center2=Pad2Loc(pad);//gives center of requested pad
316   Double_t normXmin=(x2.X()-center2.X()-PadSizeX()/2)  /Pc2Cath();//parametrise for Mathienson
317   Double_t normXmax=(x2.X()-center2.X()+PadSizeX()/2)  /Pc2Cath();
318   Double_t normYmin=(x2.Y()-center2.Y()-PadSizeY()/2)  /Pc2Cath();
319   Double_t normYmax=(x2.Y()-center2.Y()+PadSizeY()/2)  /Pc2Cath();
320  
321 //requested pad might not belong to the sector of the given hit position, hence the check:
322   return (Loc2Sec(x2)!=Pad2Sec(pad)) ? 0:Mathieson(normXmin, normYmin, normXmax, normYmax);
323 }
324 //__________________________________________________________________________________________________
325 Double_t AliRICHParam::Mathieson(Double_t x1,Double_t y1,Double_t x2,Double_t y2)
326 {
327 //This is the answer to electrostatic problem of charge distrubution in MWPC described elsewhere. (NIM A370(1988)602-603)
328 //Arguments: x1- diff between center of distribution and left margin of interested pad divided by anod-cathode distance
329 //           x2,y1,y2- analogically  
330 //  Returns: a charge fraction [0-1].
331   const Double_t kSqrtKx3=0.77459667;const Double_t kX2=0.962;const Double_t kX4=0.379;
332   const Double_t kSqrtKy3=0.77459667;const Double_t kY2=0.962;const Double_t kY4=0.379;
333
334   Double_t ux1=kSqrtKx3*TMath::TanH(kX2*x1);
335   Double_t ux2=kSqrtKx3*TMath::TanH(kX2*x2);
336   Double_t uy1=kSqrtKy3*TMath::TanH(kY2*y1);
337   Double_t uy2=kSqrtKy3*TMath::TanH(kY2*y2);
338   return 4*kX4*(TMath::ATan(ux2)-TMath::ATan(ux1))*kY4*(TMath::ATan(uy2)-TMath::ATan(uy1));
339 }
340 //__________________________________________________________________________________________________
341 TVector AliRICHParam::Loc2Area(const TVector2 &x2)
342 {
343 //Calculates the area of disintegration for a given point. It's assumed here that this points lays on anode wire.
344 //Area is a rectangulare set of pads defined by its left-down and right-up coners.
345   TVector area(4);
346   TVector pad=Loc2Pad(x2); 
347   area[0]=area[2]=pad[0]; area[1]=area[3]=pad[1];//area is just a pad fired  
348   if(pad[0]!=1           && pad[0]!= NpadsXsec()+1                            ) area[0]--; //left down coner X
349   if(pad[1]!=1           && pad[1]!= NpadsYsec()+1 && pad[1]!= 2*NpadsYsec()+1) area[1]--; //left down coner Y 
350   if(pad[0]!=NpadsXsec() && pad[0]!= NpadsX()                                 ) area[2]++; //right up coner X
351   if(pad[1]!=NpadsYsec() && pad[1]!= 2*NpadsYsec() && pad[1]!= NpadsY()       ) area[3]++; //right up coner Y
352   return area;          
353 }
354 //__________________________________________________________________________________________________
355 Bool_t AliRICHParam::IsOverTh(Int_t ,TVector ,Double_t q)
356 {
357 //Checks if the current q is over threshold and FEE will save this value to data concentrator.
358   return (q>NsigmaTh()*(SigmaThMean()+(1.-2*gRandom->Rndm())*SigmaThSpread()));
359 }
360 //__________________________________________________________________________________________________
361 #endif //AliRICHParam_h