Updating angular resolution calculation and correction factor as a function of chambe...
[u/mrichter/AliRoot.git] / HMPID / AliHMPIDPid.cxx
CommitLineData
3278403b 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
ac66a50d 26#include "AliHMPIDRecon.h" //class header
9eb30f08 27#include <AliESDtrack.h> //FindPid()
ac66a50d 28#include <TRandom.h> //Resolution()
3278403b 29
30//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
f21fc003 31AliHMPIDPid::AliHMPIDPid():TNamed("HMPIDrec","HMPIDPid")
3278403b 32{
33//..
34//init of data members
35//..
36}
37//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
b7a2d22d 38void AliHMPIDPid::FindPid(AliESDtrack *pTrk,Double_t nmean,Int_t nsp,Double_t *prob)
3278403b 39{
ac66a50d 40// Calculates probability to be a electron-muon-pion-kaon-proton with the "amplitude" method
3278403b 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
ac66a50d 44 AliPID *pPid = new AliPID();
ac66a50d 45
6ba13ae0 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
dac53a45 51 for(Int_t iPart=0;iPart<nsp;iPart++) prob[iPart]=1.0/(Float_t)nsp;
89e82ce1 52 delete pPid ; pPid=0x0; return;
9eb30f08 53 }
122305ea 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;
89e82ce1 60 delete pPid ; pPid=0x0; return;
122305ea 61 }
62
ac66a50d 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
9eb30f08 65
ac66a50d 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
b7a2d22d 72 Double_t cosThetaTh = TMath::Sqrt(mass*mass+pmod*pmod)/(nmean*pmod); // evaluate the theor. Theta Cherenkov
ac66a50d 73 if(cosThetaTh>1) continue; // no light emitted, zero height
74 Double_t thetaCerTh = TMath::ACos(cosThetaTh); // theoretical Theta Cherenkov
94d2d456 75 Double_t sigmaRing = Resolution(iPart,thetaCerTh,pTrk);
ac66a50d 76
77 if(sigmaRing==0) {
78 for(Int_t jPart=0;jPart<nsp;jPart++) prob[jPart]=1.0/(Float_t)nsp;
aa8db9d6 79 delete pPid ; pPid=0x0; delete [] h; return;
ac66a50d 80 }
81
82 if(TMath::Abs(thetaCerExp-thetaCerTh)<4*sigmaRing) desert = kFALSE; //
83 h[iPart] =TMath::Gaus(thetaCerTh,thetaCerExp,sigmaRing,kTRUE);
9eb30f08 84 hTot +=h[iPart]; //total height of all theoretical heights for normalization
ac66a50d 85
9eb30f08 86 }//species loop
87
dac53a45 88 for(Int_t iPart=0;iPart<nsp;iPart++) {//species loop to assign probabilities
ac66a50d 89
90 if(!desert) prob[iPart]=h[iPart]/hTot;
dac53a45 91 else prob[iPart]=1.0/(Float_t)nsp; //all theoretical values are far away from experemental one
ac66a50d 92
9eb30f08 93 }
ac66a50d 94
dac53a45 95 delete [] h;
ac66a50d 96 delete pPid ; pPid=0x0;
3278403b 97}
98//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
94d2d456 99Double_t AliHMPIDPid::Resolution(Int_t iPart, Double_t thetaCerTh, AliESDtrack *pTrk)
ac66a50d 100{
101 AliHMPIDParam *pParam = AliHMPIDParam::Instance();
102
103 AliHMPIDRecon rec;
b018f09d 104 Float_t xPc,yPc,thRa,phRa;
94d2d456 105 Float_t x, y;
106 Int_t q, nph;
b018f09d 107 pTrk->GetHMPIDtrk(xPc,yPc,thRa,phRa);
94d2d456 108 pTrk->GetHMPIDmip(x,y,q,nph);
b018f09d 109
110 Double_t xRa = xPc - (pParam->RadThick()+pParam->WinThick()+pParam->GapThick())*TMath::Cos(phRa)*TMath::Tan(thRa); //just linear extrapolation back to RAD
111 Double_t yRa = yPc - (pParam->RadThick()+pParam->WinThick()+pParam->GapThick())*TMath::Sin(phRa)*TMath::Tan(thRa); //just linear extrapolation back to RAD
ac66a50d 112 rec.SetTrack(xRa,yRa,thRa,phRa);
b018f09d 113
114 Double_t occupancy = pTrk->GetHMPIDoccupancy();
ac66a50d 115 Double_t thetaMax = TMath::ACos(1./pParam->MeanIdxRad());
94d2d456 116 Int_t nPhotsTh = (Int_t)(12.*TMath::Sin(thetaCerTh)*TMath::Sin(thetaCerTh)/(TMath::Sin(thetaMax)*TMath::Sin(thetaMax))+0.01);
ac66a50d 117
118 Double_t sigmatot = 0;
119 Int_t nTrks = 20;
94d2d456 120
ac66a50d 121 for(Int_t iTrk=0;iTrk<nTrks;iTrk++) {
122 Double_t invSigma = 0;
123 Int_t nPhotsAcc = 0;
94d2d456 124
125 Int_t nPhots = 0;
126 if(nph<nPhotsTh+TMath::Sqrt(nPhotsTh) && nph>nPhotsTh-TMath::Sqrt(nPhotsTh)) nPhots = nph;
127 else nPhots = gRandom->Poisson(nPhotsTh);
128
ac66a50d 129 for(Int_t j=0;j<nPhots;j++){
130 Double_t phi = gRandom->Rndm()*TMath::TwoPi();
131 TVector2 pos; pos=rec.TracePhot(thetaCerTh,phi);
132 if(!pParam->IsInside(pos.X(),pos.Y())) continue;
94d2d456 133 if(pParam->IsInDead(pos.X(),pos.Y())) continue;
ac66a50d 134 Double_t sigma2 = pParam->Sigma2(thRa,phRa,thetaCerTh,phi);//photon candidate sigma^2
135 if(sigma2!=0) {
136 invSigma += 1./sigma2;
137 nPhotsAcc++;
138 }
139 }
140 if(invSigma!=0) sigmatot += 1./TMath::Sqrt(invSigma);
141 }
94d2d456 142
143 return (sigmatot/nTrks)*pParam->SigmaCorrFact(iPart,occupancy);
ac66a50d 144}