RICH/AliRICHParam.h

   1 #ifndef AliRICHParam_h
   2 #define AliRICHParam_h
   3
   4 #include <TNamed.h>       //base class
   5 #include <TMath.h>        //QdcTot()
   6 #include <TVector.h>      //old style
   7 #include <TVector2.h>
   8 #include <TVector3.h>
   9 #include <TRandom.h>
  10 #include <TClonesArray.h>  //Hit2SDigs()
  11 #include <AliLog.h>
  12 #include <TGeoMatrix.h>    //Mars2Lors() Lors2Mars()
  13 #include <TF1.h>           //fields
  14 #include <TF2.h>           //fields
  15 #include "AliRICHDigit.h"  //Hit2Sdigs()
  16 #include <TGeoManager.h>   //Instance()
  17
  18 static const int kNchambers=7;     //number of RICH chambers
  19 static const int kNpadsX = 160;    //number of pads along X in single chamber
  20 static const int kNpadsY = 144;    //number of pads along Y in single chamber
  21 static const int kNsectors=6;      //number of sectors per chamber
  22
  23 static const int kCerenkov=50000050;  //??? go to something more general like TPDGCode
  24 static const int kFeedback=50000051;  //??? go to something more general like TPDGCode
  25
  26 // Class providing all the needed parametrised information
  27 // to construct the geometry, to define segmentation and to provide response model
  28 // In future will also provide all the staff needed for alignment and calibration
  29
  30
  31 class AliRICHParam :public TNamed
  32 {
  33 public:
  34 //ctor&dtor
  35   virtual        ~AliRICHParam()                                    {delete fIdxC6F14;fgInstance=0;}
  36 //test methodes
  37          void     Print(Option_t *opt="") const;                                         //print current parametrization
  38   static void     DrawAxis();
  39   static void     DrawSectors();
  40 //flags staff
  41   static inline AliRICHParam* Instance();                                //pointer to AliRICHParam singleton
  42   static        Int_t      Stack(Int_t evt=-1,Int_t tid=-1);              //Print stack info for event and tid
  43   static        Int_t      StackCount(Int_t pid,Int_t evt);               //Counts stack particles of given sort in given event
  44   static inline Double_t   ErrLoc                  (Double_t thetaC,Double_t phiC,Double_t thetaT,Double_t phiT,Double_t beta);
  45   static inline Double_t   ErrGeom                 (Double_t thetaC,Double_t phiC,Double_t thetaT,Double_t phiT,Double_t beta);
  46   static inline Double_t   ErrCrom                 (Double_t thetaC,Double_t phiC,Double_t thetaT,Double_t phiT,Double_t beta);
  47   static inline Double_t   SigmaSinglePhotonFormula(Double_t thetaC,Double_t phiC,Double_t thetaT,Double_t phiT,Double_t beta);
  48 //Geometrical properties
  49   static        Int_t      NpadsX      ()   {return kNpadsX;}                           //number of pads along X in chamber
  50   static        Int_t      NpadsY      ()   {return kNpadsY;}                           //number of pads along Y in chamber
  51   static        Int_t      NpadsXsec   ()   {return NpadsX()/2;}                        //number of pads along X in sector
  52   static        Int_t      NpadsYsec   ()   {return NpadsY()/3;}                        //number of pads along Y in sector
  53
  54   static        Double_t   AnodPitch   ()   {return PadSizeY()/2;}                      //cm between anode wires
  55   static        Double_t   AnodZ       ()   {return 7.806;}                             //Z positon of anod plane in LORS of the chamber, [cm]
  56   static        Double_t   CathPitch   ()   {return PadSizeY()/4;}                      //dist between cathode wires [cm]
  57   static        Double_t   CollPitch   ()   {return 0.5;}                               //dist between collection wires [cm]
  58   static        Double_t   DeadZone    ()   {return 2.6;}                               //dead zone thickness [cm]
  59   static        Double_t   PadSizeX    ()   {return 0.8;}                               //pad size x [cm]
  60   static        Double_t   PadSizeY    ()   {return 0.84;}                              //pad size y [cm]
  61   static        Double_t   PcSizeX     ()   {return NpadsX()*PadSizeX()+DeadZone();}    //PC size x [cm]
  62   static        Double_t   PcSizeY     ()   {return NpadsY()*PadSizeY()+2*DeadZone();}  //PC size y [cm]
  63   static        Double_t   Pc2Cath     ()   {return 0.445;}                             //dist between PC entrance plane and cathode wires plane [cm]
  64   static        Double_t   Pc2Win      ()   {return PcZ();}                             //dist between PC entrance plane and window exit plane [cm]
  65   static        Double_t   PcZ         ()   {return 8.0;  }                             //Z positon of PC entrance plane in LORS of the chamber [cm]
  66   static        Double_t   RadThick    ()   {return 1.5;}                               //radiator thickness [cm]
  67   static        Double_t   RadZ        ()   {return -2.0; }                             //Z positon of radiator entrance plane in LORS of the chamber [cm]
  68   static        Double_t   SecSizeX    ()   {return NpadsX()*PadSizeX()/2;}             //sector size x [cm]
  69   static        Double_t   SecSizeY    ()   {return NpadsY()*PadSizeY()/3;}             //sector size y [cm ]
  70   static        Double_t   WinThick    ()   {return 0.5;}                               //radiator window thickness [cm]
  71
  72
  73 //trasformation methodes
  74          inline TVector3   Lors2Mars     (Int_t c,Double_t x,Double_t y,Int_t p=kPc); //LORS->MARS transform of point [cm] for chamber c and plane p
  75          inline TVector3   Lors2MarsVec  (Int_t c,const TVector3 &p                ); //LORS->MARS transform of vector for chamber c
  76          inline TVector2   Mars2Lors     (Int_t c,const TVector3 &x    ,Int_t p=kPc); //MARS->LORS transform of point [cm] for chamber c and plane p
  77          inline TVector3   Mars2LorsVec  (Int_t c,const TVector3 &p                ); //MARS->LORS transform of vector for chamber c
  78
  79   static inline TVector3   Lors2MarsOld  (Int_t c,Double_t x,Double_t y,Int_t p); //LORS->MARS transform of position (cm) for chamber c and plane p
  80   static inline TVector2   Mars2LorsOld  (Int_t c,const TVector3 &x,Int_t p    ); //MARS->LORS transform of position (cm) for chamber c and plane p
  81   static inline TVector3   Center        (Int_t c,Int_t p                      ); //Center of plane p of chamber c in MARS (cm)
  82   static inline TVector3   Norm          (Int_t c                              ); //Norm vector to the chamber c in MARS (cm)
  83   static inline TGeoMatrix*Matrix        (Int_t iCh, Int_t iPlane              ); //TGeoMatrix for the given chamber plain
  84
  85   static        Int_t      Pad2Cha       (Int_t pad                      ){return pad/100000000;                     }//abs pad -> chamber
  86   static        Int_t      Pad2Sec       (Int_t pad                      ){return pad%100000000/1000000;             }//abs pad -> sector
  87   static        Int_t      Pad2PadX      (Int_t pad                      ){return pad%1000000/1000;                  }//abs pad -> pad x
  88   static        Int_t      Pad2PadY      (Int_t pad                      ){return pad%1000000%100;                   }//abs pad -> pad y
  89   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
  90   static inline TVector2   Pad2Loc       (Int_t pad                      );                                           //abs pad ->LORS
  91   static inline TVector2   Pad2Loc       (TVector pad                    );                                           //pad  -> LORS returns center of the pad
  92   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)
  93   static inline TVector    Loc2Area      (const TVector2 &x2             );                                           //pads area affected by hit x2. area is LeftDown-RightUp pad numbers
  94   static inline Int_t      Loc2Sec       (const TVector2 &x2             );                                           //LORS -> sector
  95   static        Int_t      Loc2Sec       (Double_t x,Double_t y          ){return Loc2Sec(TVector2(x,y));}            //LORS -> sector
  96   static inline TVector    Loc2Pad       (const TVector2 &x2             );                                           //LORS -> pad
  97   static        TVector    Loc2Pad       (Double_t x,Double_t y          ){return Loc2Pad(TVector2(x,y));}            //LORS -> pad
  98   static inline Int_t      Pad2Sec       (const TVector &pad             );                                           //pad  -> sector
  99   static inline Int_t      PadNeighbours (Int_t iPadX,Int_t iPadY,Int_t aListX[4],Int_t aListY[4]);                   //pad -> list of it neighbours
 100   static        Bool_t     IsAccepted    (const TVector2 &x2             ){return ( x2.X()>=0 && x2.X()<=PcSizeX() && x2.Y()>=0 && x2.Y()<=PcSizeY() ) ? kTRUE:kFALSE;}
 101 //optical properties methodes
 102   static        Float_t    EckovMean     (                               ){return 6.766e-9;}                          //mean Ckov energy according to the total trasmission curve
 103   static        Float_t    EckovMin      (                               ){return 5.5e-9;}                            //min photon energy [GeV] defined in optical curves
 104   static        Float_t    EckovMax      (                               ){return 8.5e-9;}                            //min photon energy [GeV] defined in optical curves
 105
 106                 Float_t    IdxC6F14      (Float_t gev                    ){return fIdxC6F14->Eval(gev,fIdxC6F14->GetUniqueID());}  //n=f(Eckov) [GeV] for C6H14 used as radiator
 107   static        Float_t    IdxSiO2       (Float_t gev                    ){return TMath::Sqrt(1+46.411/(10.666*10.666-gev*gev*1e18)+228.71/(18.125*18.125-gev*gev*1e18));} //n=f(Eckov) [GeV] for SiO2 used as window TDR p.35
 108   static        Float_t    IdxCH4        (Float_t gev                    ){return 1+0.12489e-6/(2.62e-4 - TMath::Power(1239.84e-9/gev,-2));}              //n=f(Eckov) [GeV] for CF4
 109   static        Float_t    AbsCH4        (Float_t gev                    );                                                                          //abs len=f(Eckov) [GeV] for CF4
 110
 111                 void       CdbRead   (Int_t run,Int_t version           );                                           //read all calibration information for requested run
 112
 113   static Double_t IonisationPotential()      {return 26.0e-9;}                            //for CH4 in GeV taken from ????
 114   static TVector2 MathiesonDelta()           {return TVector2(5*0.18,5*0.18);}            //area of 5 sigmas of Mathieson distribution (cm)
 115   static Int_t    MaxQdc()                   {return 4095;}                               //QDC number of channels
 116
 117
 118   static Int_t    QthMIP()                   {return 100;}
 119   static Double_t DmatchMIP()                {return 1.;}
 120   static Double_t PmodMax()                  {return 6.5;}
 121   static Int_t    HV(Int_t sector)           {if (sector>=1 && sector <=6) return 2050;  else return -1;} //high voltage for this sector
 122 //charge response methodes
 123   inline static Double_t Mathieson(Double_t x1,Double_t x2,Double_t y1,Double_t y2);               //Mathienson integral over given limits
 124
 125   inline static Double_t GainSag(Double_t x,Int_t sector);                                         //gain variations in %
 126          static Double_t Gain(const TVector2 &x2,Bool_t isSag=kTRUE){//gives chamber gain in terms of QDC channels for given point in local ref system
 127                           if(isSag) return QdcSlope(Loc2Sec(x2))*(1+GainSag(x2.X(),Loc2Sec(x2))/100);
 128                           else      return QdcSlope(Loc2Sec(x2));}
 129   inline static Double_t FracQdc(const TVector2 &x2,const TVector &pad);                           //charge fraction to pad from hit
 130   inline static Int_t    TotQdc(TVector2 x2,Double_t e    );                                       //total charge for Eloss (GeV) 0 for photons
 131          static Double_t QdcSlope(Int_t sec){switch(sec){case -1: return 0;  default:   return 33;}} //weight of electon in QDC channels
 132
 133   static inline Int_t    Lors2Pad      (Double_t x,Double_t y                         );                                            //LORS (x,y) [cm] -> abs pad number
 134   static        Double_t IonPot        (                                              ){return 26.0e-9;}                            //for CH4 in GeV taken from ????
 135   static inline Int_t    QdcTot        (Int_t iPad,Double_t e                         );                                            //total QDC generated by Eloss or Etot [GeV]
 136   static inline Double_t QdcSag        (Int_t iPad                                    );                                            //mean QDC variation due to sagita [0,1]
 137   static        Double_t QdcEle        (Int_t iPad,Bool_t isSag=kTRUE                 ){return isSag?33*(1+QdcSag(iPad)):33;}       //mean QDC per electron
 138   static inline Int_t    Hit2SDigs     (Int_t iPad,  Double_t e,TClonesArray* pSDigLst);                                            //hit->sdigits, returns Qtot
 139   static inline Int_t    Hit2SDigs     (TVector2 hit,Double_t e,TClonesArray* pSDigLst);                                            //hit->sdigits, returns Qtot, old style
 140   static        void     TestHit2SDigs (Double_t x,Double_t y,Double_t e,Bool_t isNew=kFALSE);                                      //test hit->sdigits
 141
 142   inline static Bool_t   IsOverTh(Int_t c,TVector pad,Double_t q);                                 //is QDC of the pad registered by FEE
 143          static Int_t    NsigmaTh()                    {return 4;}                        //
 144          static Float_t  SigmaThMean()                 {return 1.132;}                    //QDC electronic noise mean
 145          static Float_t  SigmaThSpread()               {return 0.035;}                    //QDC electronic noise width
 146
 147          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
 148                                                     -2.66575e-3*TMath::Sin(4*TMath::Pi()/PadSizeX()*x)
 149                                                     +2.80553e-3*TMath::Sin(6*TMath::Pi()/PadSizeX()*x)+0.0070;}
 150
 151          TVector3 ForwardTracing(TVector3 entranceTrackPoint,TVector3 vectorTrack, Double_t thetaC, Double_t phiC); //it traces foward a photon from Emission Point to PC
 152   static TVector3 PlaneIntersect(const TVector3 &lineDir,const TVector3 &linePoint,const TVector3 &planeNorm,const TVector3 &planePoint); //intersection between line and plane
 153   static Double_t SnellAngle(Float_t n1, Float_t n2, Float_t theta1);                              // Snell law
 154   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
 155   static        Double_t AlphaFeedback(Int_t c,Int_t s) {c++;s++; return 0.02;} //for sector s of chamber c
 156 //test part
 157   static void     Test()                            {TestSeg();TestTrans();TestResp();}  //test all groups of methodes
 158   static void     TestResp();                                                            //test the response group of methodes
 159   static void     TestSeg();                                                             //test the segmentation group of methodes
 160   static void     TestTrans();                                                           //test the transform group of methodes
 161
 162   static Double_t fgMass[5];                                // mass array
 163   enum EPlaneId {kCenter,kPc,kRad,kAnod,kNch=7};            //4 planes in chamber and total number of chambers
 164 protected:
 165          AliRICHParam();             //default ctor is protected to enforce it to be singleton
 166   static AliRICHParam *fgInstance;   //static pointer  to instance of AliRICHParam singleton
 167   TF2         *fIdxC6F14;            //n=f(Ephot,T) [GeV] for radiator freon   C6F14
 168   TGeoHMatrix *fMatrix[kNchambers];  //poiners to matrices defining RICH chambers rotations-translations
 169   ClassDef(AliRICHParam,0)           //RICH main parameters class
 170 };
 171
 172 AliRICHParam* AliRICHParam::Instance()
 173 {
 174 // Return pointer to the AliRICHParam singleton.
 175 // Arguments: none
 176 //   Returns: pointer to the instance of AliRICHParam or 0 if no geometry
 177   if(!fgInstance&&gGeoManager) new AliRICHParam;
 178   else if(!gGeoManager)                        Printf("No geometry imported");
 179   return fgInstance;
 180 }//Instance()
 181 //__________________________________________________________________________________________________
 182 Int_t AliRICHParam::PadNeighbours(Int_t iPadX,Int_t iPadY,Int_t listX[4],Int_t listY[4])
 183 {
 184 //Determines all the neighbouring pads for the given one (iPadX,iPadY). Returns total number of these pads.
 185 //Dead zones are taken into account, meaning pads from different sector are not taken.
 186 //   1
 187 // 2   3
 188 //   4
 189   Int_t nPads=0;
 190   if(iPadY!=NpadsY()&&iPadY!=2*NpadsYsec()&&iPadY!=NpadsYsec()){listX[nPads]=iPadX;   listY[nPads]=iPadY+1; nPads++;}       //1
 191   if(iPadX!=1&&iPadX!=NpadsXsec()+1)                           {listX[nPads]=iPadX-1; listY[nPads]=iPadY;   nPads++;}       //2
 192   if(iPadX!=NpadsXsec()&&iPadX!=NpadsX())                      {listX[nPads]=iPadX+1; listY[nPads]=iPadY;   nPads++;}       //3
 193   if(iPadY!=1&&iPadY!=NpadsYsec()+1&&2*NpadsYsec()+1)          {listX[nPads]=iPadX;   listY[nPads]=iPadY-1; nPads++;}       //4
 194
 195   return nPads;
 196 }//Pad2ClosePads()
 197 //__________________________________________________________________________________________________
 198 Int_t AliRICHParam::Loc2Sec(const TVector2 &v2)
 199 {
 200 // Determines sector containing the given point. y ^  5 6
 201 //                                                 |  3 4
 202 //                                                 |  1 2
 203 //                                                  -------> x
 204 // Arguments: v2- LORS position [cm]
 205 //   Returns: sector code
 206   Double_t x0=0; Double_t x1=SecSizeX(); Double_t x2=SecSizeX()+DeadZone(); Double_t x3=PcSizeX();
 207   Double_t y0=0; Double_t y1=SecSizeY(); Double_t y2=SecSizeY()+DeadZone(); Double_t y3=2*SecSizeY()+DeadZone();
 208   Double_t y4=PcSizeY()-SecSizeY();      Double_t y5=PcSizeY();
 209
 210   Int_t sector=-1;
 211   if     (v2.X() >= x0 && v2.X() <= x1 )  sector=1;
 212   else if(v2.X() >= x2 && v2.X() <= x3 )  sector=2;
 213   else                                    return -1;
 214
 215   if     (v2.Y() >= y0 && v2.Y() <= y1 )  ;                    //sectors 1 or 2
 216   else if(v2.Y() >= y2 && v2.Y() <= y3 )  sector+=2;           //sectors 3 or 4
 217   else if(v2.Y() >= y4 && v2.Y() <= y5 )  sector+=4;           //sectors 5 or 6
 218   else                                    return -1;
 219   return sector;
 220 }//Loc2Sec(Double_t x, Double_t y)
 221 //__________________________________________________________________________________________________
 222 TVector AliRICHParam::Loc2Pad(const TVector2 &loc)
 223 {
 224 //Determines pad number TVector(padx,pady) containing the given point x2 defined in the chamber RS.
 225 //Pad count starts in lower left corner from 1,1 to 144,160 in upper right corner of a chamber.
 226 //y ^  5 6
 227 //  |  3 4
 228 //  |  1 2
 229 //   -------> x
 230   TVector pad(2);
 231   Int_t sec=Loc2Sec(loc);//trasforms x2 to sector reference system
 232   if(sec==-1) {pad[0]=pad[1]=-1; return pad;}
 233 //first we deal with x
 234   if(sec==1||sec==3||sec==5)    pad[0]=           Int_t(            loc.X()   / PadSizeX() )+1; //sector 1 or 3 or 5
 235   else                          pad[0]=NpadsX() - Int_t( (PcSizeX()-loc.X())  / PadSizeX() )  ; //sector 2 or 4 or 6
 236 //second deal with y
 237        if(sec==1||sec==2)       pad[1]=Int_t(             loc.Y()                / PadSizeY())+1;               //sector 1 or 2
 238   else if(sec==3||sec==4)       pad[1]=Int_t( (loc.Y()-SecSizeY()-DeadZone()) / PadSizeY())+NpadsYsec()+1;  //sector 3 or 4
 239   else                          pad[1]=NpadsY() - Int_t( (PcSizeY()-loc.Y())     / PadSizeY());                //sector 5 or 6
 240   return pad;
 241 }
 242 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 243 Int_t AliRICHParam::Lors2Pad(Double_t x,Double_t y)
 244 {
 245 // Determines abs pad number containing the given point (x,y) defined in the chamber RS.
 246 // Pad count starts in lower left corner from 1,1 to 144,160 in upper right corner of a chamber.
 247 // y ^  5 6
 248 //   |  3 4
 249 //   |  1 2
 250 //    -------> x
 251   Int_t padx,pady;
 252   if     (x>=          0           && x<=  SecSizeX()            )  padx=  1              + Int_t(            x /PadSizeX() ); //sector 1 or 3 or 5
 253   else if(x>=SecSizeX()+DeadZone() && x<=  PcSizeX()             )  padx=     NpadsX()    - Int_t( (PcSizeX()-x)/PadSizeX() ); //sector 2 or 4 or 6
 254   else                                                              return -1;                                             //dead zone or out of chamber
 255
 256
 257   if     (y>=          0           && y<=  SecSizeY()            )  pady= 1 +               Int_t(            y /PadSizeY()              ); //sector 1 or 2
 258   else if(y>=SecSizeY()+DeadZone() && y<=2*SecSizeY()+DeadZone() )  pady= 1 + NpadsYsec() + Int_t( (y-SecSizeY()-DeadZone()) / PadSizeY()); //sector 3 or 4
 259   else if(y>= PcSizeY()-SecSizeY() && y<=  PcSizeY()             )  pady=     NpadsY()    - Int_t( (PcSizeY()-y)/PadSizeY()              ); //sector 5 or 6
 260   else                                                              return -1;                                               //dead zone or out of chamber
 261
 262   return AliRICHDigit::P2A(0,padx,pady);
 263 }//Lors2Pad()
 264 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 265 Int_t AliRICHParam::Pad2Sec(const TVector &pad)
 266 {
 267 //Determines sector containing the given pad.
 268   Int_t sector=-1;
 269   if     (pad[0] >= 1           && pad[0] <=   NpadsXsec() )    {sector=1;}
 270   else if(pad[0] >  NpadsXsec() && pad[0] <=   NpadsX()    )    {sector=2;}
 271   else                                                         AliDebugClass(1,Form("Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]));
 272
 273   if     (pad[1] >= 1             && pad[1] <=   NpadsYsec() )    {}
 274   else if(pad[1] >  NpadsYsec()   && pad[1] <= 2*NpadsYsec() )    {sector+=2;}
 275   else if(pad[1] >  2*NpadsYsec() && pad[1] <=   NpadsY()    )    {sector+=4;}
 276   else                                                         AliDebugClass(1,Form("Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]));
 277
 278   return sector;
 279 }//Pad2Sec()
 280 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 281 TVector2 AliRICHParam::Pad2Loc(TVector pad)
 282 {
 283 //Returns position of the center of the given pad in local system of the chamber (cm)
 284 // y ^  5 6
 285 //   |  3 4        sector numbers
 286 //   |  1 2
 287 //    -------> x
 288   Double_t x=-1,y=-1;
 289   if(pad[0] > 0 && pad[0] <= NpadsXsec())//it's 1 or 3 or 5
 290     x=(pad[0]-0.5)*PadSizeX();
 291   else if(pad[0] > NpadsXsec() && pad[0] <= NpadsX())//it's 2 or 4 or 6
 292     x=(pad[0]-0.5)*PadSizeX()+DeadZone();
 293   else
 294     AliDebugClass(1,Form("Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]));
 295
 296   if(pad[1] > 0 && pad[1] <= NpadsYsec())//it's 1 or 2
 297     y=(pad[1]-0.5)*PadSizeY();
 298   else if(pad[1] > NpadsYsec() && pad[1] <= 2*NpadsYsec())//it's 3 or 4
 299     y=(pad[1]-0.5)*PadSizeY()+DeadZone();
 300   else if(pad[1] > 2*NpadsYsec() && pad[1]<= NpadsY())//it's 5 or 6
 301     y=(pad[1]-0.5)*PadSizeY()+2*DeadZone();
 302   else
 303     AliDebugClass(1,Form("Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]));
 304
 305   return TVector2(x,y);
 306 }
 307 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 308 TVector2 AliRICHParam::Pad2Loc(Int_t pad)
 309 {
 310 // Converts absolute pad number to local position in LORS
 311 // LORS is a chamber  reference system with origin in left-down coner looking from IP
 312 // Arguments: pad- absolute pad number
 313 //   Returns: pad center position as TVector2 in PCRS
 314   TVector2 pos;
 315   pos.Set((Pad2PadX(pad)-0.5)*PadSizeX() , (Pad2PadY(pad)-0.5)*PadSizeY());//set to sector LORS
 316   return pos;
 317 }
 318 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 319 Double_t AliRICHParam::GainSag(Double_t x,Int_t sector)
 320 {
 321 //Returns % of gain variation due to wire sagita.
 322 //All curves are parametrized as per sector basis, so x must be apriory transformed to the Sector RS.
 323 //Here x is a distance along wires.
 324   x-=SecSizeX()/2;
 325   if(x>SecSizeX()) x-=SecSizeX();
 326   switch(HV(sector)){
 327     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;//%
 328     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;
 329     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;
 330     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;
 331     default:   return 0;
 332   }
 333 }
 334 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 335 Double_t AliRICHParam::QdcSag(Int_t iPad)
 336 {
 337 // It was observed at BNL that wires are affected by gravitation field providing a significant sagita leading to the local electric field variation
 338 // which means that different pads produce different signals.
 339 // Arguments:  iPad- absolute pad number
 340 //   Returns:  gain variation due to wire sagita 0 < QdcSag < 1.
 341 // Curves are parametrised in terms of distance x (cm) along wires having 0 on the left edge of the photocathode
 342   Double_t x=AliRICHDigit::P2X(iPad)*PadSizeX()-0.5*PadSizeX(); //center of the padx (count from 1)
 343   switch(HV(iPad)){
 344     case 2150: return 0.01*(9e-6*TMath::Power(x,4)+2e-7*TMath::Power(x,3)-0.0316*TMath::Power(x,2)-3e-4*x+25.367);//function is a fit in % so multiply by 0.01
 345     case 2100: return 0.01*(8e-6*TMath::Power(x,4)+2e-7*TMath::Power(x,3)-0.0283*TMath::Power(x,2)-2e-4*x+23.015);
 346     case 2050: return 0.01*(7e-6*TMath::Power(x,4)+1e-7*TMath::Power(x,3)-0.0254*TMath::Power(x,2)-2e-4*x+20.888);
 347     case 2000: return 0.01*(6e-6*TMath::Power(x,4)+8e-8*TMath::Power(x,3)-0.0227*TMath::Power(x,2)-1e-4*x+18.961);
 348     default:   return 0;
 349   }
 350 }//QdcSag()
 351 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 352 Int_t AliRICHParam::QdcTot(Int_t iPad,Double_t e)
 353 {
 354 // Calculates the total charge produced by the hit. Method:
 355 // 1. number of electrons is calculated as energy lost in amp gas divided by ionisation potential (for photon only one electron as Etot is always less then ionization potential)
 356 // 2. each electron imposes a charge distributed as Poisson with QdcEle() mean. Different pads produce different means. See QdcEle().
 357 // Arguments:  iPad- absolute pad number contaning the hit;
 358 //                e- Eloss for mip in amplification gas or Etot for photon
 359 //   Returns:  charge parametrised in QDC channels.
 360   Int_t iNele=Int_t(e/IonPot()); if(iNele==0) iNele=1;//e < ion. pot. means it's photoelectron
 361   Double_t dQdc=0;
 362   for(Int_t i=1;i<=iNele;i++) dQdc+=-QdcEle(iPad)*TMath::Log(gRandom->Rndm());
 363   return Int_t(dQdc);
 364 }//QdcTot()
 365 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 366 Int_t AliRICHParam::TotQdc(TVector2 x2,Double_t eloss)
 367 {
 368 //Calculates the total charge produced by the eloss in point x2 (Chamber RS).
 369 //Returns this change parametrised in QDC channels, or 0 if the hit in the dead zone.
 370 //eloss=0 means photon which produces 1 electron only eloss > 0 for Mip
 371   if(Loc2Sec(x2)==-1) return 0; //hit in the dead zone
 372   Int_t iNelectrons=Int_t(eloss/IonisationPotential()); if(iNelectrons==0) iNelectrons=1;
 373   Double_t qdc=0;
 374   for(Int_t i=1;i<=iNelectrons;i++) qdc+=-Gain(x2)*TMath::Log(gRandom->Rndm());
 375   return Int_t(qdc);
 376 }
 377 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 378 Double_t AliRICHParam::FracQdc(const TVector2 &x2,const TVector &pad)
 379 {
 380 //Calculates the charge fraction induced to given pad by the hit from the given point.
 381 //Integrated Mathieson distribution is used.
 382   TVector2 center2=Pad2Loc(pad);//gives center of requested pad
 383   Double_t normXmin=(x2.X()-center2.X()-PadSizeX()/2)  /Pc2Cath();//parametrise for Mathienson
 384   Double_t normXmax=(x2.X()-center2.X()+PadSizeX()/2)  /Pc2Cath();
 385   Double_t normYmin=(x2.Y()-center2.Y()-PadSizeY()/2)  /Pc2Cath();
 386   Double_t normYmax=(x2.Y()-center2.Y()+PadSizeY()/2)  /Pc2Cath();
 387
 388 //requested pad might not belong to the sector of the given hit position, hence the check:
 389   return (Loc2Sec(x2)!=Pad2Sec(pad)) ? 0:Mathieson(normXmin, normYmin, normXmax, normYmax);
 390 }
 391 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 392 Double_t AliRICHParam::Mathieson(Double_t x1,Double_t y1,Double_t x2,Double_t y2)
 393 {
 394 // This is the answer to electrostatic problem of charge distrubution in MWPC described elsewhere. (NIM A370(1988)602-603)
 395 // Arguments: x1- diff between center of distribution which is a hit position and left edge of interested pad divided by anod-cathode distance
 396 //            x2- right edge of the pad
 397 //            y1- up edge of the pad
 398 //            y2- bottom edge of the pad
 399 //  Returns: a charge fraction [0-1] imposed into the pad
 400   const Double_t kSqrtKx3=0.77459667;const Double_t kX2=0.962;const Double_t kX4=0.379;
 401   const Double_t kSqrtKy3=0.77459667;const Double_t kY2=0.962;const Double_t kY4=0.379;
 402
 403   Double_t ux1=kSqrtKx3*TMath::TanH(kX2*x1);
 404   Double_t ux2=kSqrtKx3*TMath::TanH(kX2*x2);
 405   Double_t uy1=kSqrtKy3*TMath::TanH(kY2*y1);
 406   Double_t uy2=kSqrtKy3*TMath::TanH(kY2*y2);
 407   return 4*kX4*(TMath::ATan(ux2)-TMath::ATan(ux1))*kY4*(TMath::ATan(uy2)-TMath::ATan(uy1));
 408 }
 409 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 410 TVector AliRICHParam::Loc2Area(const TVector2 &x2)
 411 {
 412 // Calculates the area of disintegration for a given point. It's assumed here that this points lays on anode wire.
 413 // Area is a rectangulare set of pads defined by its left-down and right-up coners.
 414   TVector area(4);
 415   TVector pad=Loc2Pad(x2);
 416   area[0]=area[2]=pad[0]; area[1]=area[3]=pad[1];//area is just a pad fired
 417   if(pad[0]!=1           && pad[0]!= NpadsXsec()+1                            ) area[0]--; //left down coner X
 418   if(pad[1]!=1           && pad[1]!= NpadsYsec()+1 && pad[1]!= 2*NpadsYsec()+1) area[1]--; //left down coner Y
 419   if(pad[0]!=NpadsXsec() && pad[0]!= NpadsX()                                 ) area[2]++; //right up coner X
 420   if(pad[1]!=NpadsYsec() && pad[1]!= 2*NpadsYsec() && pad[1]!= NpadsY()       ) area[3]++; //right up coner Y
 421   return area;
 422 }
 423 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 424 Bool_t AliRICHParam::IsOverTh(Int_t ,TVector ,Double_t qdc)
 425 {
 426 // Checks if the current QDC is over threshold and FEE will save this value to data concentrator.
 427 // This is done on pad by pad level, so the pad pedestal map is to be used. ??????????????
 428 // Arguments:
 429 //   Returns:  true if QDC over treshold
 430   return (qdc>NsigmaTh()*(SigmaThMean()+(1.-2*gRandom->Rndm())*SigmaThSpread())); //??????????? to be change to real values
 431 }
 432 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 433 TGeoMatrix* AliRICHParam::Matrix(Int_t iChamN,Int_t iPlane)
 434 {
 435   TGeoHMatrix *pMatrix=new TGeoHMatrix;
 436
 437   const Double_t kAngHor=19.5; //  horizontal angle between chambers  19.5 grad
 438   const Double_t kAngVer=20;   //  vertical angle between chambers    20   grad
 439   const Double_t kAngCom=30;   //  common RICH rotation with respect to x axis  30   grad
 440
 441   pMatrix->RotateY(90);        //rotate around y since initial position is in XY plane -> now in YZ plane
 442   Double_t trans[3]={490,0,0}; //center of the chamber is on window-gap surface
 443
 444   switch(iPlane){
 445     case kCenter:                          break;
 446     case kPc    :   trans[0]+=PcZ();       break;
 447     case kRad   :   trans[0]+=RadZ();      break;
 448     case kAnod  :   trans[0]+=AnodZ();     break;
 449     default:               return 0;       break;
 450   }
 451   pMatrix->SetTranslation(trans); //now plane in YZ is shifted along x
 452
 453   switch(iChamN){
 454     case 1:                pMatrix->RotateY(kAngHor);  pMatrix->RotateZ(-kAngVer);  break; //right and down
 455     case 2:                                            pMatrix->RotateZ(-kAngVer);  break; //down
 456     case 3:                pMatrix->RotateY(kAngHor);                               break; //right
 457     case 4:                                                                         break; //no rotation
 458     case 5:                pMatrix->RotateY(-kAngHor);                              break; //left
 459     case 6:                                            pMatrix->RotateZ(kAngVer);   break; //up
 460     case 7:                pMatrix->RotateY(-kAngHor); pMatrix->RotateZ(kAngVer);   break; //left and up
 461     default:               return 0;                                  break;
 462   }//switch(iChamber)
 463   pMatrix->RotateZ(kAngCom);     //apply common rotation  in XY plane
 464   return pMatrix;
 465 }//Matrix()
 466 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 467 TVector3  AliRICHParam::Lors2Mars(Int_t iChId,Double_t x,Double_t y,Int_t iPlnId)
 468 {
 469 // Trasform from LORS to MARS
 470 // Arguments: iChId  - chamber code 1..7
 471 //            x,y    - point in LORS
 472 //            iPlnN  - chamber plane code  might be kPc kRad kCenter kAnod
 473   Double_t z=0;
 474   switch(iPlnId){
 475     case kPc    : z=PcZ()  ; break;
 476     case kAnod  : z=AnodZ(); break;
 477     case kCenter: z=0      ; break;
 478     case kRad   : z=RadZ() ; break;
 479   }
 480   Double_t lors[3]={x-0.5*PcSizeX(),y-0.5*PcSizeY(),z},  mars[3];
 481   fMatrix[iChId-1]->LocalToMaster(lors,mars);
 482   return TVector3(mars);
 483 }
 484 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 485 TVector3  AliRICHParam::Lors2MarsVec(Int_t iCh,const TVector3 &p)
 486 {
 487   Double_t mars[3], lors[3]; p.GetXYZ(lors);
 488   fMatrix[iCh-1]->LocalToMasterVect(lors,mars);
 489   return TVector3(mars);
 490 }
 491 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 492 TVector2  AliRICHParam::Mars2Lors(Int_t iChId,const TVector3 &x,Int_t iPlnId)
 493 {
 494 // Trasform from MARS to LORS
 495 // Arguments: iChId  - chamber code 1..7
 496 //            mars   - point in MARS
 497 //            iPlnN  - chamber plane code  might be kPc kRad kCenter kAnod
 498   Double_t z=0;
 499   switch(iPlnId){
 500     case kPc    : z=PcZ()  ; break;
 501     case kAnod  : z=AnodZ(); break;
 502     case kCenter: z=0      ; break;
 503     case kRad   : z=RadZ() ; break;
 504   }
 505   Double_t lors[3],mars[3];
 506   x.GetXYZ(mars);
 507   fMatrix[iChId-1]->MasterToLocal(mars,lors);
 508   return TVector2(lors[0]+0.5*PcSizeX(),lors[1]+0.5*PcSizeY());
 509 }
 510 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 511 TVector3  AliRICHParam::Mars2LorsVec(Int_t iCh,const TVector3 &p)
 512 {
 513   Double_t mars[3], lors[3]; p.GetXYZ(mars);
 514   fMatrix[iCh-1]->MasterToLocalVect(mars,lors);
 515   return TVector3(lors);
 516 }
 517 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 518 TVector3  AliRICHParam::Lors2MarsOld(Int_t iChId,Double_t x,Double_t y,Int_t iPlnId)
 519 {
 520 // Trasform from LORS to MARS
 521 // Arguments: iChId  - chamber code 0..6
 522 //            x,y    - point in LORS
 523 //            iPlnN  - chamber plane code  might be kPc kRad kCenter kAnod
 524   TGeoMatrix *pMatrix=Matrix(iChId,iPlnId);
 525   Double_t lors[3]={x-0.5*PcSizeX(),y-0.5*PcSizeY(),0},  mars[3];  pMatrix->LocalToMaster(lors,mars);  delete pMatrix;
 526   return TVector3(mars);
 527 }
 528 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 529 TVector2  AliRICHParam::Mars2LorsOld(Int_t iChamN,const TVector3 &x,Int_t iPlaneN)
 530 {
 531   TGeoMatrix *pMatrix=Matrix(iChamN,iPlaneN);
 532   Double_t mars[3]={x.X(),x.Y(),x.Z()}  , lors[3];  pMatrix->MasterToLocal(mars,lors);  delete pMatrix;
 533   return TVector2(lors[0]+0.5*PcSizeX(),lors[1]+0.5*PcSizeY());
 534 }
 535 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 536 TVector3  AliRICHParam::Center(Int_t iChamN,Int_t iPlaneN)
 537 {
 538   TGeoMatrix *pMatrix=Matrix(iChamN,iPlaneN);
 539   Double_t mars[3]  , lors[3]={0,0,0};  pMatrix->LocalToMaster(lors,mars);  delete pMatrix;
 540   return TVector3(mars);
 541 }
 542 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 543 TVector3  AliRICHParam::Norm(Int_t iChamN)
 544 {
 545   TGeoMatrix *pMatrix=Matrix(iChamN,kPc);
 546   Double_t mars[3] , lors[3]={0,0,1};  pMatrix->LocalToMasterVect(lors,mars);  delete pMatrix;
 547   return TVector3(mars);
 548 }
 549 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 550 Int_t AliRICHParam::Hit2SDigs(Int_t iHitPad,Double_t e,TClonesArray *pSDigLst)
 551 {
 552 // Determines a number of pads affected by the hit and calculates the charge induced to each pad.
 553 // Integrated Mathieson distribution is used. Invoked from AliRICHvX::Hits2SDigits()
 554 // Arguments: iHitPad  - hit pad absolute number
 555 //            e        - energy (GeV) of this hit (Eloss for mip or Etot for photon)
 556 //            pSDigLst - pointer to clones array to store in calculated sdigits
 557 //   Returns:          total QDC for this hit
 558   Int_t iQtot=QdcTot(iHitPad,e);                    //total QDC value collected for this hit
 559   Int_t a=1;                                        //analise current pad +- a pads in both directions
 560   Int_t iLeftX=0,iBotY=0,iRightX=0,iTopY=0;         //area of disintegration for cluster formation, shifts to hit pad, not pad numbers
 561   if(AliRICHDigit::P2X(iHitPad) >                         a) iLeftX =-a;//determine area of disintegration as hit pad +- parametrised number
 562   if(AliRICHDigit::P2X(iHitPad) < AliRICHDigit::kPadsSecX-a) iRightX= a;//of pads. this number is determined by5 sigmas of Mathieson shape
 563   if(AliRICHDigit::P2Y(iHitPad) >                         a) iBotY  =-a;//see RICH TDR page 29
 564   if(AliRICHDigit::P2Y(iHitPad) < AliRICHDigit::kPadsSecY-a) iTopY  = a;//also boundary conditions are checked (edge of sector aka PC)
 565   Int_t iPadsCnt=0;
 566   for(Int_t iShiftX=iLeftX;iShiftX<=iRightX;iShiftX++){//affected pads loop iShiftX is a distance (in pads) between hit pad and pad under analisys
 567     for(Int_t iShiftY=iBotY;iShiftY<=iTopY;iShiftY++){//affected pads loop
 568       iHitPad+=AliRICHDigit::kPadAbsX*iShiftX+iShiftY;
 569       Double_t x1=PadSizeX()/Pc2Cath()*(iShiftX-0.5);//parametrise for Mathienson
 570       Double_t x2=PadSizeX()/Pc2Cath()*(iShiftX+0.5);//parametrise for Mathienson
 571       Double_t y1=PadSizeY()/Pc2Cath()*(iShiftY-0.5);//parametrise for Mathienson
 572       Double_t y2=PadSizeY()/Pc2Cath()*(iShiftY+0.5);//parametrise for Mathienson
 573       (*pSDigLst)[iPadsCnt++]= new AliRICHDigit(iHitPad,iQtot*Mathieson(x1,y1,x2,y2));
 574     }//Y loop
 575   }//X loop
 576   return iQtot;
 577 }//Hit2SDigs() for abs pad
 578 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 579 Int_t AliRICHParam::Hit2SDigs(TVector2 hitX2,Double_t e,TClonesArray *pSDigLst)
 580 {
 581 // Determines a number of pads affected by the hit and calculates the charge induced to each pad.
 582 // Integrated Mathieson distribution is used. Invoked from AliRICHvX::Hits2SDigits()
 583 // Arguments: hitX2    - hit position in LORS, cm
 584 //            e        - energy (GeV) of this hit (Eloss for mip or Etot for photon)
 585 //            pSDigLst - pointer to clones array to store in calculated sdigits
 586 //   Returns:          total QDC for this hit
 587   Int_t iQtot=TotQdc(hitX2,e);//total charge produced by hit, 0 if hit in dead zone
 588   if(iQtot==0) return 0;
 589
 590   TVector  hitPad=Loc2Pad(hitX2);  TVector2 padCenterX2=Pad2Loc(hitPad);     //shift the hit position to the nearest anod wire
 591   TVector2 anod;
 592   if((hitX2.Y()-padCenterX2.Y())>0) anod.Set(hitX2.X(),padCenterX2.Y()+AnodPitch()/2); //upper part of the pad: shift to upper anod wire
 593   else                              anod.Set(hitX2.X(),padCenterX2.Y()-AnodPitch()/2); //lower part of the pad: shift to lower anod wire
 594
 595   TVector area=Loc2Area(anod);//determine affected pads, dead zones analysed inside
 596   TVector pad(2); //current pad
 597   Int_t iPadsCnt=0;
 598   for(pad[1]=area[1];pad[1]<=area[3];pad[1]++){//affected pads loop
 599     for(pad[0]=area[0];pad[0]<=area[2];pad[0]++){
 600       Double_t dQpad=iQtot*FracQdc(anod,pad);
 601       if(dQpad>0.1) (*pSDigLst)[iPadsCnt++]= new AliRICHDigit(pad,dQpad);//make sdigit if Qpad is large enough, meaning after merging there is a chance to go above threshold
 602     }//X loop
 603   }//Y loop
 604   return iQtot;
 605 }//Hit2SDigs() for TVector2
 606 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 607 Double_t AliRICHParam::SigmaSinglePhotonFormula(Double_t thetaCer, Double_t phiCer, Double_t theta, Double_t phi, Double_t beta)
 608 {
 609   TVector3 v(-999,-999,-999);
 610
 611   v.SetX(ErrLoc(thetaCer,phiCer,theta,phi,beta));
 612   v.SetY(ErrGeom(thetaCer,phiCer,theta,phi,beta));
 613   v.SetZ(ErrCrom(thetaCer,phiCer,theta,phi,beta));
 614
 615   return v.Mag2();
 616 }
 617 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 618 Double_t AliRICHParam::ErrLoc(Double_t thetaC, Double_t phiC, Double_t Ptheta, Double_t Pphi, Double_t beta)
 619 {
 620 //par->RefIdxC6F14(par->MeanCkovEnergy())
 621 //(Float_t)1.29337525367736816e+00
 622 Double_t RefC6F14m = 1.29337;
 623   Double_t Hgap = Pc2Win();
 624   Double_t dphi = phiC - Pphi;
 625
 626   Double_t alpha =TMath::Cos(Ptheta)-TMath::Tan(thetaC)*TMath::Cos(dphi)*TMath::Sin(Ptheta);
 627   Double_t k = 1.-RefC6F14m*RefC6F14m+alpha*alpha/(beta*beta);
 628
 629   Double_t mu = TMath::Sin(Ptheta)*TMath::Sin(Pphi) + TMath::Tan(thetaC)*(TMath::Cos(Ptheta)*TMath::Cos(dphi)*TMath::Sin(Pphi)
 630 + TMath::Sin(dphi)*TMath::Cos(Pphi));
 631
 632   Double_t e = TMath::Sin(Ptheta)*TMath::Cos(Pphi)+TMath::Tan(thetaC)*(TMath::Cos(Ptheta)*TMath::Cos(dphi)*TMath::Cos(Pphi) -TMath::Sin(dphi)*TMath::Sin(Pphi));
 633
 634   Double_t kk = beta*TMath::Sqrt(k)/(Hgap*alpha);
 635   Double_t dtdxc = kk*(k*(TMath::Cos(dphi)*TMath::Cos(Pphi) - TMath::Cos(Ptheta)*TMath::Sin(dphi)*TMath::Sin(Pphi)) - ( alpha*
 636   mu/(beta*beta) )*TMath::Sin(Ptheta)*TMath::Sin(dphi));
 637
 638   Double_t dtdyc = kk*(k*(TMath::Cos(dphi)*TMath::Sin(Pphi) + TMath::Cos(Ptheta)*TMath::Sin(dphi)*TMath::Cos(Pphi)) + ( alpha*
 639   e/(beta*beta) )* TMath::Sin(Ptheta)*TMath::Sin(dphi));
 640
 641   return  TMath::Sqrt(0.2*0.2*dtdxc*dtdxc + 0.25*0.25*dtdyc*dtdyc);
 642 }
 643 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 644 Double_t AliRICHParam::ErrCrom(Double_t thetaC, Double_t phiC, Double_t Ptheta, Double_t Pphi, Double_t beta)
 645 {
 646   Double_t dphi = phiC - Pphi;
 647   Double_t RefC6F14m = 1.29337;
 648   Double_t alpha =TMath::Cos(Ptheta)-TMath::Tan(thetaC)*TMath::Cos(dphi)*TMath::Sin(Ptheta);
 649
 650   Double_t dtdn = TMath::Cos(Ptheta)*RefC6F14m*beta*beta/(alpha*TMath::Tan(thetaC));
 651
 652   Double_t f = 0.00928*(7.75-5.635)/TMath::Sqrt(12.);
 653
 654   return f*dtdn;
 655 }
 656 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 657 Double_t AliRICHParam::ErrGeom(Double_t thetaC, Double_t phiC, Double_t Ptheta, Double_t Pphi, Double_t beta )
 658 {
 659
 660   Double_t Tr = RadThick();
 661   Double_t Xep = 0.5*Tr;
 662
 663   Double_t dphi = phiC - Pphi;
 664   Double_t RefC6F14m = 1.29337;
 665   Double_t alpha =TMath::Cos(Ptheta)-TMath::Tan(thetaC)*TMath::Cos(dphi)*TMath::Sin(Ptheta);
 666
 667   Double_t k = 1.-RefC6F14m*RefC6F14m+alpha*alpha/(beta*beta);
 668
 669   Double_t Hgap = Pc2Win();
 670
 671
 672   Double_t eTr = (Tr - Xep)*beta*TMath::Sqrt(k)/(Hgap*alpha);
 673   Double_t lambda = 1.-TMath::Sin(Ptheta)*TMath::Sin(Ptheta)*TMath::Sin(phiC)*TMath::Sin(phiC);
 674
 675   Double_t c = 1./(1.+ eTr*k/(alpha*alpha*TMath::Cos(thetaC)*TMath::Cos(thetaC)));
 676   Double_t I = beta*TMath::Tan(thetaC)*lambda*TMath::Power(k,1.5);
 677   Double_t II = 1.+eTr*beta*I;
 678
 679   Double_t err = c * (I/(alpha*alpha*Hgap) +  II* (1.-lambda) / ( alpha*alpha*Hgap*beta*(1.+eTr)) );
 680   Double_t TrErr = Tr/(TMath::Sqrt(12.)*TMath::Cos(Ptheta));
 681
 682   return TrErr*err;
 683 }//ErrGeom()
 684
 685 #endif //AliRICHParam_h