LCOV - code coverage report
Current view: top level - STEER/STEERBase - AliHMPIDPIDResponse.cxx (source / functions) Hit Total Coverage
Test: coverage.info Lines: 42 267 15.7 %
Date: 2016-06-14 17:26:59 Functions: 5 28 17.9 %

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

Generated by: LCOV version 1.11