STEER/STEERBase/AliHMPIDPIDResponse.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //***********************************************************
  17 // Class AliHMPIDPIDResponse
  18 //
  19 // HMPID class to perfom particle identification
  20 //
  21 // Author: G. Volpe, giacomo.volpe@cern.ch
  22 //***********************************************************
  23
  24
  25 #include "AliHMPIDPIDResponse.h"  //class header
  26 #include "AliPID.h"               //FindPid()
  27 #include "AliVTrack.h"            //FindPid()
  28 #include "AliLog.h"               //general
  29 #include <TRandom.h>              //Resolution()
  30 #include <TVector2.h>             //Resolution()
  31 #include <TRotation.h>
  32 #include <TF1.h>
  33 #include <TGeoManager.h>          //Instance()
  34 #include <TGeoMatrix.h>           //Instance()
  35 #include <TGeoPhysicalNode.h>     //ctor
  36 #include <TGeoBBox.h>
  37 #include <TObjArray.h>
  38
  39 Float_t AliHMPIDPIDResponse::fgkMinPcX[]={0.,0.,0.,0.,0.,0.};
  40 Float_t AliHMPIDPIDResponse::fgkMaxPcX[]={0.,0.,0.,0.,0.,0.};
  41 Float_t AliHMPIDPIDResponse::fgkMinPcY[]={0.,0.,0.,0.,0.,0.};
  42 Float_t AliHMPIDPIDResponse::fgkMaxPcY[]={0.,0.,0.,0.,0.,0.};
  43
  44 Float_t AliHMPIDPIDResponse::fgCellX=0.;
  45 Float_t AliHMPIDPIDResponse::fgCellY=0.;
  46
  47 Float_t AliHMPIDPIDResponse::fgPcX=0;
  48 Float_t AliHMPIDPIDResponse::fgPcY=0;
  49
  50 Float_t AliHMPIDPIDResponse::fgAllX=0;
  51 Float_t AliHMPIDPIDResponse::fgAllY=0;
  52
  53 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  54 AliHMPIDPIDResponse::AliHMPIDPIDResponse():
  55   TNamed("HMPIDPIDResponseRec","HMPIDPIDResponsePid"),
  56   fRefIdx(1.28947),
  57   fTrkDir(0,0,1),
  58   fTrkPos(30,40),
  59   fRefIndexArray(0x0)
  60 {
  61   //
  62   // ctor
  63   //
  64
  65   Float_t dead=2.6;// cm of the dead zones between PCs-> See 2CRC2099P1
  66
  67   fgCellX=0.8; fgCellY=0.84;
  68
  69   fgPcX  = 80.*fgCellX;      fgPcY = 48.*fgCellY;
  70   fgAllX = 2.*fgPcX+dead;
  71   fgAllY = 3.*fgPcY+2.*dead;
  72
  73   fgkMinPcX[1]=fgPcX+dead; fgkMinPcX[3]=fgkMinPcX[1];  fgkMinPcX[5]=fgkMinPcX[3];
  74   fgkMaxPcX[0]=fgPcX;      fgkMaxPcX[2]=fgkMaxPcX[0];  fgkMaxPcX[4]=fgkMaxPcX[2];
  75   fgkMaxPcX[1]=fgAllX;     fgkMaxPcX[3]=fgkMaxPcX[1];  fgkMaxPcX[5]=fgkMaxPcX[3];
  76
  77   fgkMinPcY[2]=fgPcY+dead;       fgkMinPcY[3]=fgkMinPcY[2];
  78   fgkMinPcY[4]=2.*fgPcY+2.*dead; fgkMinPcY[5]=fgkMinPcY[4];
  79   fgkMaxPcY[0]=fgPcY;            fgkMaxPcY[1]=fgkMaxPcY[0];
  80   fgkMaxPcY[2]=2.*fgPcY+dead;    fgkMaxPcY[3]=fgkMaxPcY[2];
  81   fgkMaxPcY[4]=fgAllY;           fgkMaxPcY[5]=fgkMaxPcY[4];
  82
  83   for(Int_t i=kMinCh;i<=kMaxCh;i++)
  84     if(gGeoManager && gGeoManager->IsClosed()) {
  85       TGeoPNEntry* pne = gGeoManager->GetAlignableEntry(Form("/HMPID/Chamber%i",i));
  86       if (!pne) {
  87         AliErrorClass(Form("The symbolic volume %s does not correspond to any physical entry!",Form("HMPID_%i",i)));
  88         fM[i]=new TGeoHMatrix;
  89         IdealPosition(i,fM[i]);
  90       } else {
  91         TGeoPhysicalNode *pnode = pne->GetPhysicalNode();
  92         if(pnode) fM[i]=new TGeoHMatrix(*(pnode->GetMatrix()));
  93         else {
  94           fM[i]=new TGeoHMatrix;
  95           IdealPosition(i,fM[i]);
  96         }
  97       }
  98     } else{
  99       fM[i]=new TGeoHMatrix;
 100       IdealPosition(i,fM[i]);
 101     }
 102
 103 }//ctor
 104 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 105 AliHMPIDPIDResponse::AliHMPIDPIDResponse(const AliHMPIDPIDResponse& c):
 106   TNamed(c),
 107   fRefIdx(c.fRefIdx),
 108   fTrkDir(c.fTrkDir),
 109   fTrkPos(c.fTrkPos),
 110   fRefIndexArray(c.fRefIndexArray)
 111  {
 112    //
 113    // copy ctor
 114    //
 115
 116    for(Int_t i=0; i<6; i++) {
 117
 118       fgkMinPcX[i] = c.fgkMinPcX[i];
 119       fgkMinPcY[i] = c.fgkMinPcY[i];
 120       fgkMaxPcX[i] = c.fgkMaxPcX[i];
 121       fgkMaxPcY[i] = c.fgkMaxPcY[i];
 122      }
 123
 124    for(Int_t i=0; i<7; i++) fM[i] = c.fM[i] ? new TGeoHMatrix(*c.fM[i]) : 0;
 125  }
 126 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 127 AliHMPIDPIDResponse::~AliHMPIDPIDResponse()
 128 {
 129   // d-tor
 130   for (int i=7;i--;) delete fM[i];
 131 }
 132
 133 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 134 AliHMPIDPIDResponse& AliHMPIDPIDResponse::operator=(const AliHMPIDPIDResponse& c) {
 135
 136    //
 137    // assignment operator
 138    //
 139   if(this!=&c){
 140      TNamed::operator=(c);
 141      fgCellX = c.fgCellX;
 142      fgCellY = c.fgCellY;
 143      fgPcX   = c.fgPcX;
 144      fgPcY   = c.fgPcY;
 145      fgAllX  = c.fgAllX;
 146      fgAllY  = c.fgAllY;
 147      fRefIdx = c.fRefIdx;
 148      fTrkDir = c.fTrkDir;
 149      fTrkPos = c.fTrkPos;
 150      fRefIndexArray = c.fRefIndexArray;
 151      for(Int_t i=0; i<6; i++) {
 152       fgkMinPcX[i] = c.fgkMinPcX[i];
 153       fgkMinPcY[i] = c.fgkMinPcY[i];
 154       fgkMaxPcX[i] = c.fgkMaxPcX[i];
 155       fgkMaxPcY[i] = c.fgkMaxPcY[i];
 156      }
 157      for(Int_t i=0; i<7; i++) fM[i] = c.fM[i] ? new TGeoHMatrix(*c.fM[i]) : 0;
 158     }
 159
 160   return *this;
 161 }
 162 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 163 void AliHMPIDPIDResponse::IdealPosition(Int_t iCh, TGeoHMatrix *pMatrix) {
 164
 165 // Construct ideal position matrix for a given chamber
 166 // Arguments: iCh- chamber ID; pMatrix- pointer to precreated unity matrix where to store the results
 167 // Returns: none
 168
 169   const Double_t kAngHor=19.5;        //  horizontal angle between chambers  19.5 grad
 170   const Double_t kAngVer=20;          //  vertical angle between chambers    20   grad
 171   const Double_t kAngCom=30;          //  common HMPID rotation with respect to x axis  30   grad
 172   const Double_t kTrans[3]={490,0,0}; //  center of the chamber is on window-gap surface
 173   pMatrix->RotateY(90);               //  rotate around y since initial position is in XY plane -> now in YZ plane
 174   pMatrix->SetTranslation(kTrans);    //  now plane in YZ is shifted along x
 175   switch(iCh){
 176     case 0:                pMatrix->RotateY(kAngHor);  pMatrix->RotateZ(-kAngVer);  break; //right and down
 177     case 1:                                            pMatrix->RotateZ(-kAngVer);  break; //down
 178     case 2:                pMatrix->RotateY(kAngHor);                               break; //right
 179     case 3:                                                                         break; //no rotation
 180     case 4:                pMatrix->RotateY(-kAngHor);                              break; //left
 181     case 5:                                            pMatrix->RotateZ(kAngVer);   break; //up
 182     case 6:                pMatrix->RotateY(-kAngHor); pMatrix->RotateZ(kAngVer);   break; //left and up
 183   }
 184   pMatrix->RotateZ(kAngCom);     //apply common rotation  in XY plane
 185
 186 }
 187 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 188 Double_t AliHMPIDPIDResponse::GetExpectedSignal(const AliVTrack *vTrk, AliPID::EParticleType specie) const {
 189
 190   // expected Cherenkov angle calculation
 191
 192   const Double_t nmean = GetNMean(vTrk);
 193   return ExpectedSignal(vTrk,nmean,specie);
 194 }
 195 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 196 Double_t AliHMPIDPIDResponse::GetExpectedSigma(const AliVTrack *vTrk, AliPID::EParticleType specie) const {
 197
 198   // expected resolution calculation
 199
 200   const Double_t nmean = GetNMean(vTrk);
 201   return ExpectedSigma(vTrk,nmean,specie);
 202 }
 203 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 204 Double_t AliHMPIDPIDResponse::ExpectedSignal(const AliVTrack *vTrk, Double_t nmean, AliPID::EParticleType specie) const {
 205
 206   // expected Cherenkov angle calculation
 207
 208   Double_t thetaTheor = -999.;
 209
 210   Double_t p[3] = {0}, mom = 0;
 211   if(vTrk->GetOuterHmpPxPyPz(p))  mom = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);  // Momentum of the charged particle
 212   else return thetaTheor;
 213
 214   const Double_t mass = AliPID::ParticleMass(specie);
 215   const Double_t cosTheta = TMath::Sqrt(mass*mass+mom*mom)/(nmean*mom);
 216
 217   if(cosTheta>1) return thetaTheor;
 218
 219   else thetaTheor = TMath::ACos(cosTheta);
 220
 221   return thetaTheor;                                                                                          //  evaluate the theor. Theta Cherenkov
 222 }
 223 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 224 Double_t AliHMPIDPIDResponse::ExpectedSigma(const AliVTrack *vTrk, Double_t nmean, AliPID::EParticleType specie) const {
 225
 226   // expected resolution calculation
 227
 228   Float_t x=0., y=0.;
 229   Int_t q=0, nph=0;
 230   Float_t xPc=0.,yPc=0.,thRa=0.,phRa=0.;
 231
 232   vTrk->GetHMPIDmip(x,y,q,nph);
 233   vTrk->GetHMPIDtrk(xPc,yPc,thRa,phRa);
 234
 235   const Double_t xRa = xPc - (RadThick()+WinThick()+GapThick())*TMath::Cos(phRa)*TMath::Tan(thRa);  //just linear extrapolation back to RAD
 236   const Double_t yRa = yPc - (RadThick()+WinThick()+GapThick())*TMath::Sin(phRa)*TMath::Tan(thRa);  //just linear extrapolation back to RAD
 237
 238   const Double_t thetaCerTh = ExpectedSignal(vTrk,nmean,specie);
 239   const Double_t occupancy  = vTrk->GetHMPIDoccupancy();
 240   const Double_t thetaMax   = TMath::ACos(1./nmean);
 241   const Int_t    nPhotsTh   = (Int_t)(12.*TMath::Sin(thetaCerTh)*TMath::Sin(thetaCerTh)/(TMath::Sin(thetaMax)*TMath::Sin(thetaMax))+0.01);
 242
 243   Double_t sigmatot = 0;
 244   Int_t nTrks = 20;
 245   for(Int_t iTrk=0;iTrk<nTrks;iTrk++) {
 246     Double_t invSigma = 0;
 247     Int_t nPhotsAcc = 0;
 248
 249     Int_t nPhots = 0;
 250     if(nph<nPhotsTh+TMath::Sqrt(nPhotsTh) && nph>nPhotsTh-TMath::Sqrt(nPhotsTh)) nPhots = nph;
 251     else nPhots = gRandom->Poisson(nPhotsTh);
 252
 253     for(Int_t j=0;j<nPhots;j++){
 254       Double_t phi = gRandom->Rndm()*TMath::TwoPi();
 255       TVector2 pos; pos = TracePhot(xRa,yRa,thRa,phRa,thetaCerTh,phi);
 256       if(!IsInside(pos.X(),pos.Y())) continue;
 257       if(IsInDead(pos.X(),pos.Y()))  continue;
 258       Double_t sigma2 = Sigma2(thRa,phRa,thetaCerTh,phi); //photon candidate sigma^2
 259
 260       if(sigma2!=0) {
 261         invSigma += 1./sigma2;
 262         nPhotsAcc++;
 263       }
 264     }
 265     if(invSigma!=0) sigmatot += 1./TMath::Sqrt(invSigma);
 266   }
 267
 268   return (sigmatot/nTrks)*SigmaCorrFact(specie,occupancy);
 269 }
 270 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 271 Double_t AliHMPIDPIDResponse::GetNumberOfSigmas(const AliVTrack *vTrk, AliPID::EParticleType specie) const {
 272
 273   // Number of sigmas calculation
 274
 275   Double_t nSigmas = -999.;
 276
 277   if(vTrk->GetHMPIDsignal()<0.) return nSigmas;
 278
 279   const Double_t nmean = GetNMean(vTrk);
 280
 281   const Double_t expSigma = ExpectedSigma(vTrk, nmean, specie);
 282
 283   if(expSigma > 0.) nSigmas = (vTrk->GetHMPIDsignal() - ExpectedSignal(vTrk,nmean,specie))/expSigma;
 284
 285   return nSigmas;
 286
 287 }
 288 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 289 void AliHMPIDPIDResponse::GetProbability(const AliVTrack *vTrk,Int_t nSpecies,Double_t *prob) const {
 290
 291 // Calculates probability to be a electron-muon-pion-kaon-proton with the "amplitude" method
 292 // from the given Cerenkov angle and momentum assuming no initial particle composition
 293
 294   const Double_t thetaCerExp = vTrk->GetHMPIDsignal();
 295
 296   const Double_t nmean = GetNMean(vTrk);
 297
 298   if(thetaCerExp<=0){                                                                                     // HMPID does not find anything reasonable for this track, assign 0.2 for all species
 299     for(Int_t iPart=0;iPart<nSpecies;iPart++) prob[iPart]=1.0/(Float_t)nSpecies;
 300     return;
 301   }
 302
 303   Double_t p[3] = {0,0,0};
 304
 305   if(!(vTrk->GetOuterHmpPxPyPz(p))) for(Int_t iPart=0;iPart<nSpecies;iPart++) prob[iPart]=1.0/(Float_t)nSpecies;
 306
 307   Double_t hTot=0;                                                                                        // Initialize the total height of the amplitude method
 308   Double_t *h = new Double_t [nSpecies];                                                                  // number of charged particles to be considered
 309
 310   Bool_t desert = kTRUE;                                                                                  // Flag to evaluate if ThetaC is far ("desert") from the given Gaussians
 311
 312   for(Int_t iPart=0;iPart<nSpecies;iPart++){                                                              // for each particle
 313
 314
 315     h[iPart] = 0;                                                                                         // reset the height
 316     Double_t thetaCerTh = ExpectedSignal(vTrk,nmean,(AliPID::EParticleType)iPart);                        // theoretical Theta Cherenkov
 317     if(thetaCerTh>900.) continue;                                                                         // no light emitted, zero height
 318     Double_t sigmaRing = ExpectedSigma(vTrk,nmean,(AliPID::EParticleType)iPart);
 319
 320     if(sigmaRing==0) continue;
 321
 322     if(TMath::Abs(thetaCerExp-thetaCerTh)<4*sigmaRing) desert = kFALSE;
 323     h[iPart] =TMath::Gaus(thetaCerTh,thetaCerExp,sigmaRing,kTRUE);
 324     hTot    +=h[iPart];                                                                                   // total height of all theoretical heights for normalization
 325
 326   }//species loop
 327
 328   for(Int_t iPart=0;iPart<nSpecies;iPart++) {                                                             // species loop to assign probabilities
 329
 330     if(!desert) prob[iPart]=h[iPart]/hTot;
 331     else        prob[iPart]=1.0/(Float_t)nSpecies;                                                        // all theoretical values are far away from experemental one
 332
 333   }
 334
 335   delete [] h;
 336 }
 337 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 338 Double_t AliHMPIDPIDResponse::GetSignalDelta(const AliVTrack *vTrk, AliPID::EParticleType specie, Bool_t ratio/*=kFALSE*/) const {
 339
 340   //
 341   // calculation of Experimental Cherenkov angle - Theoretical Cherenkov angle
 342   //
 343   const Double_t signal    = vTrk->GetHMPIDsignal();
 344   const Double_t expSignal = GetExpectedSignal(vTrk,specie);
 345
 346   Double_t delta = -9999.;
 347   if (!ratio) delta=signal-expSignal;
 348   else if (expSignal>1.e-20) delta=signal/expSignal;
 349
 350   return delta;
 351 }
 352 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 353 TVector2 AliHMPIDPIDResponse::TracePhot(Double_t xRa, Double_t yRa, Double_t thRa, Double_t phRa, Double_t ckovThe,Double_t ckovPhi) const {
 354
 355 // Trace a single Ckov photon from emission point somewhere in radiator up to photocathode taking into account ref indexes of materials it travereses
 356 // Returns: distance between photon point on PC and track projection
 357
 358   Double_t theta=0.,phi=0.;
 359   TVector3  dirTRS,dirLORS;
 360   dirTRS.SetMagThetaPhi(1,ckovThe,ckovPhi);                     //photon in TRS
 361   Trs2Lors(thRa,phRa,dirTRS,theta,phi);
 362   dirLORS.SetMagThetaPhi(1,theta,phi);                          //photon in LORS
 363   return TraceForward(xRa,yRa,dirLORS);                                 //now foward tracing
 364 }//TracePhot()
 365 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 366 TVector2 AliHMPIDPIDResponse::TraceForward(Double_t xRa, Double_t yRa, TVector3 dirCkov) const {
 367
 368 // Trace forward a photon from (x,y) up to PC
 369 // Returns: pos of traced photon at PC
 370
 371   TVector2 pos(-999,-999);
 372   Double_t thetaCer = dirCkov.Theta();
 373   if(thetaCer > TMath::ASin(1./GetRefIdx())) return pos;          //total refraction on WIN-GAP boundary
 374   Double_t zRad= -0.5*RadThick()-0.5*WinThick();          //z position of middle of RAD
 375   TVector3  posCkov(xRa,yRa,zRad);                        //RAD: photon position is track position @ middle of RAD
 376   Propagate(dirCkov,posCkov,           -0.5*WinThick());          //go to RAD-WIN boundary
 377   Refract  (dirCkov,         GetRefIdx(),WinIdx());       //RAD-WIN refraction
 378   Propagate(dirCkov,posCkov,            0.5*WinThick());          //go to WIN-GAP boundary
 379   Refract  (dirCkov,         WinIdx(),GapIdx());          //WIN-GAP refraction
 380   Propagate(dirCkov,posCkov,0.5*WinThick()+GapThick());   //go to PC
 381   pos.Set(posCkov.X(),posCkov.Y());
 382   return pos;
 383 }//TraceForward()
 384 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 385 void AliHMPIDPIDResponse::Propagate(const TVector3 dir,TVector3 &pos,Double_t z) const {
 386
 387 // Finds an intersection point between a line and XY plane shifted along Z.
 388 // Arguments:  dir,pos   - vector along the line and any point of the line
 389 //             z         - z coordinate of plain
 390 // Returns:  none
 391 // On exit:  pos is the position if this intesection if any
 392
 393   static TVector3 nrm(0,0,1);
 394          TVector3 pnt(0,0,z);
 395
 396   TVector3 diff=pnt-pos;
 397   Double_t sint=(nrm*diff)/(nrm*dir);
 398   pos+=sint*dir;
 399 }//Propagate()
 400 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 401 void AliHMPIDPIDResponse::Refract(TVector3 &dir,Double_t n1,Double_t n2) const {
 402
 403 // Refract direction vector according to Snell law
 404 // Arguments:
 405 //            n1 - ref idx of first substance
 406 //            n2 - ref idx of second substance
 407 //   Returns: none
 408 //   On exit: dir is new direction
 409
 410   Double_t sinref=(n1/n2)*TMath::Sin(dir.Theta());
 411   if(TMath::Abs(sinref)>1.) dir.SetXYZ(-999,-999,-999);
 412   else             dir.SetTheta(TMath::ASin(sinref));
 413 }//Refract()
 414 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 415 void AliHMPIDPIDResponse::Trs2Lors(Double_t thRa, Double_t phRa, TVector3 dirCkov,Double_t &thetaCer,Double_t &phiCer) const {
 416
 417   // Theta Cerenkov reconstruction
 418   // Returns: thetaCer of photon in LORS
 419   //          phiCer of photon in LORS
 420
 421   TRotation mtheta;   mtheta.RotateY(thRa);
 422   TRotation mphi;       mphi.RotateZ(phRa);
 423   TRotation mrot=mphi*mtheta;
 424   TVector3 dirCkovLORS;
 425   dirCkovLORS=mrot*dirCkov;
 426   phiCer  = dirCkovLORS.Phi();                                          //actual value of the phi of the photon
 427   thetaCer= dirCkovLORS.Theta();                                        //actual value of thetaCerenkov of the photon
 428 }
 429 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 430 Bool_t AliHMPIDPIDResponse::IsInDead(Float_t x,Float_t y)  {
 431
 432 // Check is the current point is outside of sensitive area or in dead zones
 433 // Arguments: x,y -position
 434 // Returns: 1 if not in sensitive zone
 435
 436   for(Int_t iPc=0;iPc<6;iPc++)
 437     if(x>=fgkMinPcX[iPc] && x<=fgkMaxPcX[iPc] && y>=fgkMinPcY[iPc] && y<=fgkMaxPcY [iPc]) return kFALSE; //in current pc
 438
 439   return kTRUE;
 440 }
 441 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 442 Double_t AliHMPIDPIDResponse::Sigma2(Double_t trkTheta,Double_t trkPhi,Double_t ckovTh, Double_t ckovPh) const {
 443
 444 // Analithical calculation of total error (as a sum of localization, geometrical and chromatic errors) on Cerenkov angle for a given Cerenkov photon
 445 // created by a given MIP. Formules according to CERN-EP-2000-058
 446 // Arguments: Cerenkov and azimuthal angles for Cerenkov photon, [radians]
 447 //            dip and azimuthal angles for MIP taken at the entrance to radiator, [radians]
 448 //            MIP beta
 449 // Returns: absolute error on Cerenkov angle, [radians]
 450
 451   TVector3 v(-999,-999,-999);
 452   Double_t trkBeta = 1./(TMath::Cos(ckovTh)*GetRefIdx());
 453
 454   if(trkBeta > 1) trkBeta = 1;                 //protection against bad measured thetaCer
 455   if(trkBeta < 0) trkBeta = 0.0001;            //
 456
 457   v.SetX(SigLoc (trkTheta,trkPhi,ckovTh,ckovPh,trkBeta));
 458   v.SetY(SigGeom(trkTheta,trkPhi,ckovTh,ckovPh,trkBeta));
 459   v.SetZ(SigCrom(trkTheta,ckovTh,ckovPh,trkBeta));
 460
 461   return v.Mag2();
 462 }
 463 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 464 Double_t AliHMPIDPIDResponse::SigLoc(Double_t trkTheta,Double_t trkPhi,Double_t thetaC, Double_t phiC,Double_t betaM) const {
 465
 466 // Analitical calculation of localization error (due to finite segmentation of PC) on Cerenkov angle for a given Cerenkov photon
 467 // created by a given MIP. Fromulae according to CERN-EP-2000-058
 468 // Arguments: Cerenkov and azimuthal angles for Cerenkov photon, [radians]
 469 //            dip and azimuthal angles for MIP taken at the entrance to radiator, [radians]
 470 //            MIP beta
 471 // Returns: absolute error on Cerenkov angle, [radians]
 472
 473   Double_t phiDelta = phiC;
 474
 475   Double_t sint     = TMath::Sin(trkTheta);
 476   Double_t cost     = TMath::Cos(trkTheta);
 477   Double_t sinf     = TMath::Sin(trkPhi);
 478   Double_t cosf     = TMath::Cos(trkPhi);
 479   Double_t sinfd    = TMath::Sin(phiDelta);
 480   Double_t cosfd    = TMath::Cos(phiDelta);
 481   Double_t tantheta = TMath::Tan(thetaC);
 482
 483   Double_t alpha =cost-tantheta*cosfd*sint;                                                    // formula (11)
 484   Double_t k = 1.-GetRefIdx()*GetRefIdx()+alpha*alpha/(betaM*betaM);                           // formula (after 8 in the text)
 485   if (k<0) return 1e10;
 486   Double_t mu =sint*sinf+tantheta*(cost*cosfd*sinf+sinfd*cosf);                                // formula (10)
 487   Double_t e  =sint*cosf+tantheta*(cost*cosfd*cosf-sinfd*sinf);                                // formula (9)
 488
 489   Double_t kk = betaM*TMath::Sqrt(k)/(GapThick()*alpha);                                       // formula (6) and (7)
 490   Double_t dtdxc = kk*(k*(cosfd*cosf-cost*sinfd*sinf)-(alpha*mu/(betaM*betaM))*sint*sinfd);    // formula (6)
 491   Double_t dtdyc = kk*(k*(cosfd*sinf+cost*sinfd*cosf)+(alpha* e/(betaM*betaM))*sint*sinfd);    // formula (7)            pag.4
 492
 493   Double_t errX = 0.2,errY=0.25;                                                                //end of page 7
 494   return  TMath::Sqrt(errX*errX*dtdxc*dtdxc + errY*errY*dtdyc*dtdyc);
 495 }
 496 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 497 Double_t AliHMPIDPIDResponse::SigCrom(Double_t trkTheta,Double_t thetaC, Double_t phiC,Double_t betaM) const {
 498
 499 // Analitical calculation of chromatic error (due to lack of knowledge of Cerenkov photon energy) on Cerenkov angle for a given Cerenkov photon
 500 // created by a given MIP. Fromulae according to CERN-EP-2000-058
 501 // Arguments: Cerenkov and azimuthal angles for Cerenkov photon, [radians]
 502 //            dip and azimuthal angles for MIP taken at the entrance to radiator, [radians]
 503 //            MIP beta
 504 //   Returns: absolute error on Cerenkov angle, [radians]
 505
 506   Double_t phiDelta = phiC;
 507
 508   Double_t sint     = TMath::Sin(trkTheta);
 509   Double_t cost     = TMath::Cos(trkTheta);
 510   Double_t cosfd    = TMath::Cos(phiDelta);
 511   Double_t tantheta = TMath::Tan(thetaC);
 512
 513   Double_t alpha =cost-tantheta*cosfd*sint;                                         // formula (11)
 514   Double_t dtdn = cost*GetRefIdx()*betaM*betaM/(alpha*tantheta);                    // formula (12)
 515
 516   Double_t f = 0.0172*(7.75-5.635)/TMath::Sqrt(24.);
 517
 518   return f*dtdn;
 519 }//SigCrom()
 520 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 521 Double_t AliHMPIDPIDResponse::SigGeom(Double_t trkTheta,Double_t trkPhi,Double_t thetaC, Double_t phiC,Double_t betaM) const {
 522
 523 // Analitical calculation of geometric error (due to lack of knowledge of creation point in radiator) on Cerenkov angle for a given Cerenkov photon
 524 // created by a given MIP. Formulae according to CERN-EP-2000-058
 525 // Arguments: Cerenkov and azimuthal angles for Cerenkov photon, [radians]
 526 //            dip and azimuthal angles for MIP taken at the entrance to radiator, [radians]
 527 //            MIP beta
 528 //   Returns: absolute error on Cerenkov angle, [radians]
 529
 530   Double_t phiDelta = phiC;
 531
 532   Double_t sint     = TMath::Sin(trkTheta);
 533   Double_t cost     = TMath::Cos(trkTheta);
 534   Double_t sinf     = TMath::Sin(trkPhi);
 535   Double_t cosfd    = TMath::Cos(phiDelta);
 536   Double_t costheta = TMath::Cos(thetaC);
 537   Double_t tantheta = TMath::Tan(thetaC);
 538
 539   Double_t alpha =cost-tantheta*cosfd*sint;                                                // formula (11)
 540
 541   Double_t k = 1.-GetRefIdx()*GetRefIdx()+alpha*alpha/(betaM*betaM);                       // formula (after 8 in the text)
 542   if (k<0) return 1e10;
 543
 544   Double_t eTr = 0.5*RadThick()*betaM*TMath::Sqrt(k)/(GapThick()*alpha);                    // formula (14)
 545   Double_t lambda = (1.-sint*sinf)*(1.+sint*sinf);                                          // formula (15)
 546
 547   Double_t c1 = 1./(1.+ eTr*k/(alpha*alpha*costheta*costheta));                             // formula (13.a)
 548   Double_t c2 = betaM*TMath::Power(k,1.5)*tantheta*lambda/(GapThick()*alpha*alpha);         // formula (13.b)
 549   Double_t c3 = (1.+eTr*k*betaM*betaM)/((1+eTr)*alpha*alpha);                               // formula (13.c)
 550   Double_t c4 = TMath::Sqrt(k)*tantheta*(1-lambda)/(GapThick()*betaM);                      // formula (13.d)
 551   Double_t dtdT = c1 * (c2+c3*c4);
 552   Double_t trErr = RadThick()/(TMath::Sqrt(12.)*cost);
 553
 554   return trErr*dtdT;
 555 }//SigGeom()
 556 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 557 Double_t AliHMPIDPIDResponse::GetNMean(const AliVTrack *vTrk) const {
 558
 559   //
 560   // mean refractive index calculation
 561   //
 562   Double_t nmean = -999.;
 563
 564   Float_t xPc=0.,yPc=0.,thRa=0.,phRa=0.;
 565   vTrk->GetHMPIDtrk(xPc,yPc,thRa,phRa);
 566
 567   const Int_t ch = vTrk->GetHMPIDcluIdx()/1000000;
 568
 569   const Double_t yRa = yPc - (RadThick()+WinThick()+GapThick())*TMath::Sin(phRa)*TMath::Tan(thRa);  //just linear extrapolation back to RAD
 570
 571   TF1 *RefIndex=0x0;
 572
 573   if(GetRefIndexArray()) RefIndex = (TF1*)(GetRefIndexArray()->At(ch));
 574   else return nmean;
 575
 576   if(RefIndex) nmean = RefIndex->Eval(yRa);
 577   else return nmean;
 578
 579   return nmean;
 580 }
 581 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 582 Double_t AliHMPIDPIDResponse::SigmaCorrFact  (Int_t iPart, Double_t occupancy) {
 583
 584 // calculation of sigma correction factor
 585
 586   Double_t corr = 1.0;
 587
 588   switch(iPart) {
 589     case 0: corr = 0.115*occupancy + 1.166; break;
 590     case 1: corr = 0.115*occupancy + 1.166; break;
 591     case 2: corr = 0.115*occupancy + 1.166; break;
 592     case 3: corr = 0.065*occupancy + 1.137; break;
 593     case 4: corr = 0.048*occupancy + 1.202; break;
 594   }
 595  return corr;
 596 }
 597