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