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