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 : // //
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 8 : AliHMPIDPid::AliHMPIDPid():TNamed("HMPIDrec","HMPIDPid")
32 40 : {
33 : //..
34 : //init of data members
35 : //..
36 16 : }
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 16 : AliPID *pPid = new AliPID();
45 :
46 : Double_t thetaCerExp = -999.;
47 8 : if(pTrk->GetHMPIDsignal()<=0) thetaCerExp = pTrk->GetHMPIDsignal();
48 8 : else thetaCerExp = pTrk->GetHMPIDsignal() - (Int_t)pTrk->GetHMPIDsignal(); // measured thetaCherenkov
49 :
50 8 : if(thetaCerExp<=0){ //HMPID does not find anything reasonable for this track, assign 0.2 for all species
51 0 : for(Int_t iPart=0;iPart<nsp;iPart++) prob[iPart]=1.0/(Float_t)nsp;
52 0 : delete pPid ; pPid=0x0; return;
53 : }
54 :
55 8 : Double_t p[3] = {0}, pmod = 0;
56 16 : 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 0 : for(Int_t iPart=0;iPart<nsp;iPart++) prob[iPart]=1.0/(Float_t)nsp;
60 0 : delete pPid ; pPid=0x0; return;
61 : }
62 :
63 : Double_t hTot=0; // Initialize the total height of the amplitude method
64 8 : 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 96 : for(Int_t iPart=0;iPart<nsp;iPart++){ // for each particle
69 :
70 40 : h[iPart] = 0; // reset the height
71 40 : Double_t mass = pPid->ParticleMass(iPart); // with the given mass
72 40 : Double_t cosThetaTh = TMath::Sqrt(mass*mass+pmod*pmod)/(nmean*pmod); // evaluate the theor. Theta Cherenkov
73 44 : if(cosThetaTh>1) continue; // no light emitted, zero height
74 36 : Double_t thetaCerTh = TMath::ACos(cosThetaTh); // theoretical Theta Cherenkov
75 36 : Double_t sigmaRing = Resolution(iPart,thetaCerTh,pTrk);
76 :
77 36 : if(sigmaRing==0) continue;
78 :
79 60 : if(TMath::Abs(thetaCerExp-thetaCerTh)<4*sigmaRing) desert = kFALSE; //
80 36 : h[iPart] =TMath::Gaus(thetaCerTh,thetaCerExp,sigmaRing,kTRUE);
81 36 : hTot +=h[iPart]; //total height of all theoretical heights for normalization
82 :
83 36 : }//species loop
84 :
85 96 : for(Int_t iPart=0;iPart<nsp;iPart++) {//species loop to assign probabilities
86 :
87 80 : if(!desert) prob[iPart]=h[iPart]/hTot;
88 0 : else prob[iPart]=1.0/(Float_t)nsp; //all theoretical values are far away from experemental one
89 :
90 : }
91 :
92 16 : delete [] h;
93 16 : delete pPid ; pPid=0x0;
94 24 : }
95 : //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
96 : Double_t AliHMPIDPid::Resolution(Int_t iPart, Double_t thetaCerTh, AliESDtrack *pTrk)
97 : {
98 72 : AliHMPIDParam *pParam = AliHMPIDParam::Instance();
99 :
100 36 : AliHMPIDRecon rec;
101 36 : Float_t xPc,yPc,thRa,phRa;
102 36 : Float_t x, y;
103 36 : Int_t q, nph;
104 36 : pTrk->GetHMPIDtrk(xPc,yPc,thRa,phRa);
105 36 : pTrk->GetHMPIDmip(x,y,q,nph);
106 :
107 36 : Double_t xRa = xPc - (pParam->RadThick()+pParam->WinThick()+pParam->GapThick())*TMath::Cos(phRa)*TMath::Tan(thRa); //just linear extrapolation back to RAD
108 36 : Double_t yRa = yPc - (pParam->RadThick()+pParam->WinThick()+pParam->GapThick())*TMath::Sin(phRa)*TMath::Tan(thRa); //just linear extrapolation back to RAD
109 36 : rec.SetTrack(xRa,yRa,thRa,phRa);
110 :
111 36 : Double_t occupancy = pTrk->GetHMPIDoccupancy();
112 72 : Double_t thetaMax = TMath::ACos(1./pParam->MeanIdxRad());
113 36 : 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 1512 : 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 1840 : if(nph<nPhotsTh+TMath::Sqrt(nPhotsTh) && nph>nPhotsTh-TMath::Sqrt(nPhotsTh)) nPhots = nph;
124 220 : else nPhots = gRandom->Poisson(nPhotsTh);
125 :
126 14748 : for(Int_t j=0;j<nPhots;j++){
127 13308 : Double_t phi = gRandom->Rndm()*TMath::TwoPi();
128 19962 : TVector2 pos; pos=rec.TracePhot(thetaCerTh,phi);
129 7723 : if(!pParam->IsInside(pos.X(),pos.Y())) continue;
130 5730 : if(pParam->IsInDead(pos.X(),pos.Y())) continue;
131 5440 : Double_t sigma2 = pParam->Sigma2(thRa,phRa,thetaCerTh,phi);//photon candidate sigma^2
132 5440 : if(sigma2!=0) {
133 5440 : invSigma += 1./sigma2;
134 5440 : nPhotsAcc++;
135 5440 : }
136 12094 : }
137 1427 : if(invSigma!=0) sigmatot += 1./TMath::Sqrt(invSigma);
138 : }
139 :
140 72 : return (sigmatot/nTrks)*pParam->SigmaCorrFact(iPart,occupancy);
141 36 : }
|
