PID probability calculated using refractive index evaluated from OCDB information...
[u/mrichter/AliRoot.git] / 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(thetaCerTh,pTrk);
76     
77     if(sigmaRing==0) {
78       for(Int_t jPart=0;jPart<nsp;jPart++) prob[jPart]=1.0/(Float_t)nsp;
79       delete pPid ; pPid=0x0; delete [] h; return;
80     } 
81     
82     if(TMath::Abs(thetaCerExp-thetaCerTh)<4*sigmaRing) desert = kFALSE;                                                                //   
83     h[iPart] =TMath::Gaus(thetaCerTh,thetaCerExp,sigmaRing,kTRUE);
84     hTot    +=h[iPart]; //total height of all theoretical heights for normalization
85     
86   }//species loop
87
88   for(Int_t iPart=0;iPart<nsp;iPart++) {//species loop to assign probabilities
89     
90     if(!desert) prob[iPart]=h[iPart]/hTot;
91     else prob[iPart]=1.0/(Float_t)nsp;            //all theoretical values are far away from experemental one
92     
93   }
94   
95   delete [] h;
96   delete pPid ; pPid=0x0;
97 }
98 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
99 Double_t AliHMPIDPid::Resolution(Double_t thetaCerTh, AliESDtrack *pTrk)
100 {
101   AliHMPIDParam *pParam = AliHMPIDParam::Instance();
102       
103   AliHMPIDRecon rec;
104   Float_t xRa,yRa,thRa,phRa;
105   pTrk->GetHMPIDtrk(xRa,yRa,thRa,phRa);
106   rec.SetTrack(xRa,yRa,thRa,phRa);
107   Double_t thetaMax = TMath::ACos(1./pParam->MeanIdxRad());
108   Int_t nPhots = (Int_t)(21.*TMath::Sin(thetaCerTh)*TMath::Sin(thetaCerTh)/(TMath::Sin(thetaMax)*TMath::Sin(thetaMax))+0.01);
109
110   Double_t sigmatot = 0;
111   Int_t nTrks = 20;
112   for(Int_t iTrk=0;iTrk<nTrks;iTrk++) {
113     Double_t invSigma = 0;
114     Int_t nPhotsAcc = 0;
115     for(Int_t j=0;j<nPhots;j++){
116       Double_t phi = gRandom->Rndm()*TMath::TwoPi();
117       TVector2 pos; pos=rec.TracePhot(thetaCerTh,phi);
118       if(!pParam->IsInside(pos.X(),pos.Y())) continue;
119       if(pParam->IsInDead(pos.X(),pos.Y())) continue;
120       Double_t sigma2 = pParam->Sigma2(thRa,phRa,thetaCerTh,phi);//photon candidate sigma^2
121       if(sigma2!=0) {
122         invSigma += 1./sigma2;
123         nPhotsAcc++;
124       }
125     }      
126     if(invSigma!=0) sigmatot += 1./TMath::Sqrt(invSigma);  
127   }
128   return sigmatot/nTrks;
129 }