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