HMPID/AliHMPIDPid.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 //                                                                      //
  18 // AliHMPIDPid                                                          //
  19 //                                                                      //
  20 // HMPID class to perfom particle identification                        //
  21 //                                                                      //
  22 //////////////////////////////////////////////////////////////////////////
  23
  24 #include "AliHMPIDPid.h"       //class header
  25 #include "AliHMPIDParam.h"     //class header
  26 #include "AliHMPIDRecon.h"     //class header
  27 #include <AliESDtrack.h>       //FindPid()
  28 #include <TRandom.h>           //Resolution()
  29
  30 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  31 AliHMPIDPid::AliHMPIDPid():TNamed("HMPIDrec","HMPIDPid")
  32 {
  33 //..
  34 //init of data members
  35 //..
  36 }
  37 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  38 void AliHMPIDPid::FindPid(AliESDtrack *pTrk,Double_t nmean,Int_t nsp,Double_t *prob)
  39 {
  40 // Calculates probability to be a electron-muon-pion-kaon-proton with the "amplitude" method
  41 // from the given Cerenkov angle and momentum assuming no initial particle composition
  42 // (i.e. apriory probability to be the particle of the given sort is the same for all sorts)
  43
  44   AliPID *pPid = new AliPID();
  45
  46   Double_t thetaCerExp = -999.;
  47   if(pTrk->GetHMPIDsignal()<=0) thetaCerExp = pTrk->GetHMPIDsignal();
  48   else                          thetaCerExp = pTrk->GetHMPIDsignal() - (Int_t)pTrk->GetHMPIDsignal();     //  measured thetaCherenkov
  49
  50   if(thetaCerExp<=0){                                                                                     //HMPID does not find anything reasonable for this track, assign 0.2 for all species
  51     for(Int_t iPart=0;iPart<nsp;iPart++) prob[iPart]=1.0/(Float_t)nsp;
  52     delete pPid ; pPid=0x0; return;
  53   }
  54
  55   Double_t p[3] = {0}, pmod = 0;
  56   if(pTrk->GetOuterHmpPxPyPz(p))  pmod = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);  // Momentum of the charged particle
  57
  58   else {
  59     for(Int_t iPart=0;iPart<nsp;iPart++) prob[iPart]=1.0/(Float_t)nsp;
  60     delete pPid ; pPid=0x0; return;
  61   }
  62
  63   Double_t hTot=0;                               // Initialize the total height of the amplitude method
  64   Double_t *h = new Double_t [nsp];              // number of charged particles to be considered
  65
  66   Bool_t desert = kTRUE;                                                                                                     //  Flag to evaluate if ThetaC is far ("desert") from the given Gaussians
  67
  68   for(Int_t iPart=0;iPart<nsp;iPart++){                                                                                      //  for each particle
  69
  70     h[iPart] = 0;                                                                                                            //  reset the height
  71     Double_t mass = pPid->ParticleMass(iPart);                                                                             //  with the given mass
  72     Double_t cosThetaTh = TMath::Sqrt(mass*mass+pmod*pmod)/(nmean*pmod);                   //  evaluate the theor. Theta Cherenkov
  73     if(cosThetaTh>1) continue;                                                                                               //  no light emitted, zero height
  74     Double_t thetaCerTh = TMath::ACos(cosThetaTh);                                                                           //  theoretical Theta Cherenkov
  75     Double_t sigmaRing = Resolution(iPart,thetaCerTh,pTrk);
  76
  77     if(sigmaRing==0) continue;
  78
  79     if(TMath::Abs(thetaCerExp-thetaCerTh)<4*sigmaRing) desert = kFALSE;                                                                //
  80     h[iPart] =TMath::Gaus(thetaCerTh,thetaCerExp,sigmaRing,kTRUE);
  81     hTot    +=h[iPart]; //total height of all theoretical heights for normalization
  82
  83   }//species loop
  84
  85   for(Int_t iPart=0;iPart<nsp;iPart++) {//species loop to assign probabilities
  86
  87     if(!desert) prob[iPart]=h[iPart]/hTot;
  88     else prob[iPart]=1.0/(Float_t)nsp;            //all theoretical values are far away from experemental one
  89
  90   }
  91
  92   delete [] h;
  93   delete pPid ; pPid=0x0;
  94 }
  95 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  96 Double_t AliHMPIDPid::Resolution(Int_t iPart, Double_t thetaCerTh, AliESDtrack *pTrk)
  97 {
  98   AliHMPIDParam *pParam = AliHMPIDParam::Instance();
  99
 100   AliHMPIDRecon rec;
 101   Float_t xPc,yPc,thRa,phRa;
 102   Float_t x, y;
 103   Int_t q, nph;
 104   pTrk->GetHMPIDtrk(xPc,yPc,thRa,phRa);
 105   pTrk->GetHMPIDmip(x,y,q,nph);
 106
 107   Double_t xRa = xPc - (pParam->RadThick()+pParam->WinThick()+pParam->GapThick())*TMath::Cos(phRa)*TMath::Tan(thRa); //just linear extrapolation back to RAD
 108   Double_t yRa = yPc - (pParam->RadThick()+pParam->WinThick()+pParam->GapThick())*TMath::Sin(phRa)*TMath::Tan(thRa); //just linear extrapolation back to RAD
 109   rec.SetTrack(xRa,yRa,thRa,phRa);
 110
 111   Double_t occupancy = pTrk->GetHMPIDoccupancy();
 112   Double_t thetaMax = TMath::ACos(1./pParam->MeanIdxRad());
 113   Int_t    nPhotsTh = (Int_t)(12.*TMath::Sin(thetaCerTh)*TMath::Sin(thetaCerTh)/(TMath::Sin(thetaMax)*TMath::Sin(thetaMax))+0.01);
 114
 115   Double_t sigmatot = 0;
 116   Int_t nTrks = 20;
 117
 118   for(Int_t iTrk=0;iTrk<nTrks;iTrk++) {
 119     Double_t invSigma = 0;
 120     Int_t nPhotsAcc = 0;
 121
 122     Int_t nPhots = 0;
 123     if(nph<nPhotsTh+TMath::Sqrt(nPhotsTh) && nph>nPhotsTh-TMath::Sqrt(nPhotsTh)) nPhots = nph;
 124     else nPhots = gRandom->Poisson(nPhotsTh);
 125
 126     for(Int_t j=0;j<nPhots;j++){
 127       Double_t phi = gRandom->Rndm()*TMath::TwoPi();
 128       TVector2 pos; pos=rec.TracePhot(thetaCerTh,phi);
 129       if(!pParam->IsInside(pos.X(),pos.Y())) continue;
 130       if(pParam->IsInDead(pos.X(),pos.Y()))  continue;
 131       Double_t sigma2 = pParam->Sigma2(thRa,phRa,thetaCerTh,phi);//photon candidate sigma^2
 132       if(sigma2!=0) {
 133         invSigma += 1./sigma2;
 134         nPhotsAcc++;
 135       }
 136     }
 137     if(invSigma!=0) sigmatot += 1./TMath::Sqrt(invSigma);
 138   }
 139
 140   return (sigmatot/nTrks)*pParam->SigmaCorrFact(iPart,occupancy);
 141 }