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