added protection in the destructor, checking if OCDB cache is active
[u/mrichter/AliRoot.git] / HMPID / AliHMPIDRecon.cxx
CommitLineData
d3da6dc4 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// AliHMPIDRecon //
19// //
20// HMPID class to perfom pattern recognition based on Hough transfrom //
21// for single chamber //
22//////////////////////////////////////////////////////////////////////////
23
a591e55f 24#include "AliHMPIDRecon.h" //class header
d3da6dc4 25#include "AliHMPIDCluster.h" //CkovAngle()
a591e55f 26#include <TRotation.h> //TracePhot()
27#include <TH1D.h> //HoughResponse()
28#include <TClonesArray.h> //CkovAngle()
29#include <AliESDtrack.h> //CkovAngle()
d3da6dc4 30
d3da6dc4 31//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
ffb1ac19 32AliHMPIDRecon::AliHMPIDRecon():TTask("RichRec","RichPat")
d3da6dc4 33{
ffb1ac19 34//..
35//init of data members
36//..
37
38 fPhotCnt = -1;
39 fPhotFlag = 0x0;
40 fPhotCkov = 0x0;
41 fPhotPhi = 0x0;
42 fPhotWei = 0x0;
43 fCkovSigma2 = 0;
44 fIsWEIGHT = kFALSE;
45 fDTheta = 0.001;
46 fWindowWidth = 0.045;
47 fTrkDir = TVector3(0,0,1); // init just for test
48 fTrkPos = TVector2(30,40); // init just for test
49
50 AliHMPIDParam *pParam=AliHMPIDParam::Instance();
51 fParam = pParam;
52
53 fParam->SetRefIdx(fParam->MeanIdxRad()); // initialization of ref index to a default one
54}
55//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
56void AliHMPIDRecon::InitVars(Int_t n)
57{
58//..
59//Init some variables
60//..
61 if(n<0) return;
62 fPhotFlag = new Int_t[n];
63 fPhotCkov = new Double_t[n];
64 fPhotPhi = new Double_t[n];
65 fPhotWei = new Double_t[n];
66//
67}
68//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
69void AliHMPIDRecon::DeleteVars()
70{
71//..
72//Delete variables
73//..
74 delete fPhotFlag;
75 delete fPhotCkov;
76 delete fPhotPhi;
77 delete fPhotWei;
d3da6dc4 78}
79//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
afe12692 80void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t nmean,Double_t qthre)
d3da6dc4 81{
82// Pattern recognition method based on Hough transform
59280a5a 83// Arguments: pTrk - track for which Ckov angle is to be found
84// pCluLst - list of clusters for this chamber
85// Returns: - track ckov angle, [rad],
a591e55f 86
ffb1ac19 87 Int_t nClusTot = pCluLst->GetEntries();
88 if(nClusTot>fParam->MultCut()) fIsWEIGHT = kTRUE; // offset to take into account bkg in reconstruction
89 else fIsWEIGHT = kFALSE;
d3da6dc4 90
ffb1ac19 91 InitVars(nClusTot);
92
611e810d 93 Float_t xRa,yRa,th,ph;
a591e55f 94 pTrk->GetHMPIDtrk(xRa,yRa,th,ph); //initialize this track: th and ph angles at middle of RAD
a591e55f 95 SetTrack(xRa,yRa,th,ph);
611e810d 96
ffb1ac19 97 fParam->SetRefIdx(nmean);
d3da6dc4 98
59280a5a 99 Float_t dMin=999,mipX=-1,mipY=-1;Int_t chId=-1,mipId=-1,mipQ=-1;
d3da6dc4 100 fPhotCnt=0;
101 for (Int_t iClu=0; iClu<pCluLst->GetEntriesFast();iClu++){//clusters loop
102 AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(iClu); //get pointer to current cluster
59280a5a 103 chId=pClu->Ch();
afe12692 104 if(pClu->Q()>qthre){ //charge compartible with MIP clusters
a591e55f 105 Float_t dX=fPc.X()-pClu->X(),dY=fPc.Y()-pClu->Y(),d =TMath::Sqrt(dX*dX+dY*dY); //distance between current cluster and intersection point
106 if( d < dMin) {mipId=iClu; dMin=d;mipX=pClu->X();mipY=pClu->Y();mipQ=(Int_t)pClu->Q();} //current cluster is closer, overwrite data for min cluster
107 }else{ //charge compatible with photon cluster
108 Double_t thetaCer,phiCer;
109 if(FindPhotCkov(pClu->X(),pClu->Y(),thetaCer,phiCer)){ //find ckov angle for this photon candidate
110 fPhotCkov[fPhotCnt]=thetaCer; //actual theta Cerenkov (in TRS)
b4ad85e9 111 fPhotPhi [fPhotCnt]=phiCer; //actual phi Cerenkov (in TRS): -pi to come back to "unusual" ref system (X,Y,-Z)
28500fe1 112 //PH Printf("photon n. %i reconstructed theta = %f",fPhotCnt,fPhotCkov[fPhotCnt]);
a591e55f 113 fPhotCnt++; //increment counter of photon candidates
114 }
59280a5a 115 }
d3da6dc4 116 }//clusters loop
4598109f 117 fMipPos.Set(mipX,mipY);
76fd1a96 118 if(fPhotCnt<=3) pTrk->SetHMPIDsignal(kNoPhotAccept); //no reconstruction with <=3 photon candidates
a591e55f 119 Int_t iNacc=FlagPhot(HoughResponse()); //flag photons according to individual theta ckov with respect to most probable
120 pTrk->SetHMPIDmip(mipX,mipY,mipQ,iNacc); //store mip info
59280a5a 121
a591e55f 122 if(mipId==-1) {pTrk->SetHMPIDsignal(kMipQdcCut); return;} //no clusters with QDC more the threshold at all
ffb1ac19 123 if(dMin>fParam->DistCut()) {pTrk->SetHMPIDsignal(kMipDistCut); return;} //closest cluster with enough charge is still too far from intersection
a591e55f 124 pTrk->SetHMPIDcluIdx(chId,mipId); //set index of cluster
76fd1a96 125 if(iNacc<1){
126 pTrk->SetHMPIDsignal(kNoPhotAccept); //no photon candidates is accepted
127 }
128 else {
129 pTrk->SetHMPIDsignal(FindRingCkov(pCluLst->GetEntries())); //find best Theta ckov for ring i.e. track
130 pTrk->SetHMPIDchi2(fCkovSigma2); //errors squared
131 }
d3da6dc4 132
ffb1ac19 133 DeleteVars();
43400d2d 134}//CkovAngle()
d3da6dc4 135//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
a591e55f 136Bool_t AliHMPIDRecon::FindPhotCkov(Double_t cluX,Double_t cluY,Double_t &thetaCer,Double_t &phiCer)
d3da6dc4 137{
138// Finds Cerenkov angle for this photon candidate
139// Arguments: cluX,cluY - position of cadidate's cluster
a591e55f 140// Returns: Cerenkov angle
d3da6dc4 141
a591e55f 142 TVector3 dirCkov;
143
ffb1ac19 144 Double_t zRad= -0.5*fParam->RadThick()-0.5*fParam->WinThick(); //z position of middle of RAD
67a1c24c 145 TVector3 rad(fTrkPos.X(),fTrkPos.Y(),zRad); //impact point at middle of RAD
ffb1ac19 146 TVector3 pc(cluX,cluY,0.5*fParam->WinThick()+fParam->GapIdx()); //mip at PC
a591e55f 147 Double_t cluR = TMath::Sqrt((cluX-fTrkPos.X())*(cluX-fTrkPos.X())+
148 (cluY-fTrkPos.Y())*(cluY-fTrkPos.Y()));//ref. distance impact RAD-CLUSTER
67a1c24c 149 Double_t phi=(pc-rad).Phi(); //phi of photon
a591e55f 150
b4ad85e9 151 Double_t ckov1=0;
67a1c24c 152 Double_t ckov2=0.75+fTrkDir.Theta(); //start to find theta cerenkov in DRS
b4ad85e9 153 const Double_t kTol=0.01;
d3da6dc4 154 Int_t iIterCnt = 0;
155 while(1){
a591e55f 156 if(iIterCnt>=50) return kFALSE;
d3da6dc4 157 Double_t ckov=0.5*(ckov1+ckov2);
67a1c24c 158 dirCkov.SetMagThetaPhi(1,ckov,phi);
a591e55f 159 TVector2 posC=TraceForward(dirCkov); //trace photon with actual angles
160 Double_t dist=cluR-(posC-fTrkPos).Mod(); //get distance between trial point and cluster position
161 if(posC.X()==-999) dist = - 999; //total reflection problem
162 iIterCnt++; //counter step
b4ad85e9 163 if (dist> kTol) ckov1=ckov; //cluster @ larger ckov
d3da6dc4 164 else if(dist<-kTol) ckov2=ckov; //cluster @ smaller ckov
a591e55f 165 else{ //precision achived: ckov in DRS found
166 dirCkov.SetMagThetaPhi(1,ckov,phi); //
167 RecPhot(dirCkov,thetaCer,phiCer); //find ckov (in TRS:the effective Cherenkov angle!)
168 return kTRUE;
169 }
d3da6dc4 170 }
171}//FindPhotTheta()
172//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
a591e55f 173TVector2 AliHMPIDRecon::TraceForward(TVector3 dirCkov)const
d3da6dc4 174{
a591e55f 175 //Trace forward a photon from (x,y) up to PC
176 // Arguments: dirCkov photon vector in LORS
177 // Returns: pos of traced photon at PC
ffb1ac19 178
a591e55f 179 TVector2 pos(-999,-999);
67a1c24c 180 Double_t thetaCer = dirCkov.Theta();
ffb1ac19 181 if(thetaCer > TMath::ASin(1./fParam->GetRefIdx())) return pos; //total refraction on WIN-GAP boundary
182 Double_t zRad= -0.5*fParam->RadThick()-0.5*fParam->WinThick(); //z position of middle of RAD
183 TVector3 posCkov(fTrkPos.X(),fTrkPos.Y(),zRad); //RAD: photon position is track position @ middle of RAD
184 Propagate(dirCkov,posCkov, -0.5*fParam->WinThick()); //go to RAD-WIN boundary
185 Refract (dirCkov, fParam->GetRefIdx(),fParam->WinIdx()); //RAD-WIN refraction
186 Propagate(dirCkov,posCkov, 0.5*fParam->WinThick()); //go to WIN-GAP boundary
187 Refract (dirCkov, fParam->WinIdx(),fParam->GapIdx()); //WIN-GAP refraction
188 Propagate(dirCkov,posCkov,0.5*fParam->WinThick()+fParam->GapThick()); //go to PC
a591e55f 189 pos.Set(posCkov.X(),posCkov.Y());
190 return pos;
191}//TraceForward()
192//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
193void AliHMPIDRecon::RecPhot(TVector3 dirCkov,Double_t &thetaCer,Double_t &phiCer)
194{
195 //Theta Cerenkov reconstruction
196 // Arguments: (x,y) of initial point in LORS, dirCkov photon vector in LORS
197 // Returns: thetaCer theta cerenkov reconstructed
198// TVector3 dirTrk;
199// dirTrk.SetMagThetaPhi(1,fTrkDir.Theta(),fTrkDir.Phi());
200// Double_t thetaCer = TMath::ACos(dirCkov*dirTrk);
201 TRotation mtheta; mtheta.RotateY(- fTrkDir.Theta());
202 TRotation mphi; mphi.RotateZ(- fTrkDir.Phi());
203 TRotation mrot=mtheta*mphi;
204 TVector3 dirCkovTRS;
205 dirCkovTRS=mrot*dirCkov;
206 phiCer = dirCkovTRS.Phi(); //actual value of the phi of the photon
207 thetaCer= dirCkovTRS.Theta(); //actual value of thetaCerenkov of the photon
d3da6dc4 208}
209//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
4598109f 210Double_t AliHMPIDRecon::FindRingArea(Double_t ckovAng)const
d3da6dc4 211{
4598109f 212// Find area covered in the PC acceptance
213// Arguments: ckovAng - cerenkov angle
d3da6dc4 214// Returns: area of the ring in cm^2 for given theta ckov
215
ffb1ac19 216 const Int_t kN=50;
4598109f 217 TVector2 pos1;
afe12692 218 Double_t area=0;
4598109f 219 Bool_t first=kFALSE;
afe12692 220 for(Int_t i=0;i<kN;i++){
4598109f 221 if(!first) {
222 pos1=TracePhot(ckovAng,Double_t(TMath::TwoPi()*(i+1)/kN)); //find a good trace for the first photon
223 if(pos1.X()==-999) continue; //no area: open ring
ffb1ac19 224 if(!fParam->IsInside(pos1.X(),pos1.Y(),0)) pos1 = IntWithEdge(fMipPos,pos1); // ffind the very first intersection...
4598109f 225 first=kTRUE;
226 continue;
227 }
228 TVector2 pos2=TracePhot(ckovAng,Double_t(TMath::TwoPi()*(i+1)/kN)); //trace the next photon
229 if(pos2.X()==-999) continue; //no area: open ring
ffb1ac19 230 if(!fParam->IsInside(pos2.X(),pos2.Y(),0)) {
4598109f 231 pos2 = IntWithEdge(fMipPos,pos2);
232 }
233 area+=TMath::Abs((pos1-fMipPos).X()*(pos2-fMipPos).Y()-(pos1-fMipPos).Y()*(pos2-fMipPos).X()); //add area of the triangle...
234 pos1 = pos2;
d3da6dc4 235 }
4598109f 236//--- find points from ring
7fc88c5e 237 area*=0.5;
d3da6dc4 238 return area;
239}//FindRingArea()
240//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
4598109f 241TVector2 AliHMPIDRecon::IntWithEdge(TVector2 p1,TVector2 p2)const
242{
243// It finds the intersection of the line for 2 points traced as photons
244// and the edge of a given PC
245// Arguments: 2 points obtained tracing the photons
246// Returns: intersection point with detector (PC) edges
247
4598109f 248 Double_t xmin = (p1.X()<p2.X())? p1.X():p2.X();
249 Double_t xmax = (p1.X()<p2.X())? p2.X():p1.X();
250 Double_t ymin = (p1.Y()<p2.Y())? p1.Y():p2.Y();
251 Double_t ymax = (p1.Y()<p2.Y())? p2.Y():p1.Y();
252
253 Double_t m = TMath::Tan((p2-p1).Phi());
254 TVector2 pint;
255 //intersection with low X
256 pint.Set((Double_t)(p1.X() + (0-p1.Y())/m),0.);
ffb1ac19 257 if(pint.X()>=0 && pint.X()<=fParam->SizeAllX() &&
4598109f 258 pint.X()>=xmin && pint.X()<=xmax &&
259 pint.Y()>=ymin && pint.Y()<=ymax) return pint;
260 //intersection with high X
ffb1ac19 261 pint.Set((Double_t)(p1.X() + (fParam->SizeAllY()-p1.Y())/m),(Double_t)(fParam->SizeAllY()));
262 if(pint.X()>=0 && pint.X()<=fParam->SizeAllX() &&
4598109f 263 pint.X()>=xmin && pint.X()<=xmax &&
264 pint.Y()>=ymin && pint.Y()<=ymax) return pint;
265 //intersection with left Y
266 pint.Set(0.,(Double_t)(p1.Y() + m*(0-p1.X())));
ffb1ac19 267 if(pint.Y()>=0 && pint.Y()<=fParam->SizeAllY() &&
4598109f 268 pint.Y()>=ymin && pint.Y()<=ymax &&
269 pint.X()>=xmin && pint.X()<=xmax) return pint;
270 //intersection with righ Y
ffb1ac19 271 pint.Set((Double_t)(fParam->SizeAllX()),(Double_t)(p1.Y() + m*(fParam->SizeAllX()-p1.X())));
272 if(pint.Y()>=0 && pint.Y()<=fParam->SizeAllY() &&
4598109f 273 pint.Y()>=ymin && pint.Y()<=ymax &&
274 pint.X()>=xmin && pint.X()<=xmax) return pint;
275 return p1;
276}//IntWithEdge()
277//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
d3da6dc4 278Double_t AliHMPIDRecon::FindRingCkov(Int_t)
279{
280// Loops on all Ckov candidates and estimates the best Theta Ckov for a ring formed by those candidates. Also estimates an error for that Theat Ckov
281// collecting errors for all single Ckov candidates thetas. (Assuming they are independent)
282// Arguments: iNclus- total number of clusters in chamber for background estimation
283// Return: best estimation of track Theta ckov
284
285 Double_t wei = 0.;
286 Double_t weightThetaCerenkov = 0.;
287
288 Double_t ckovMin=9999.,ckovMax=0.;
289 Double_t sigma2 = 0; //to collect error squared for this ring
290
291 for(Int_t i=0;i<fPhotCnt;i++){//candidates loop
292 if(fPhotFlag[i] == 2){
a591e55f 293 if(fPhotCkov[i]<ckovMin) ckovMin=fPhotCkov[i]; //find max and min Theta ckov from all candidates within probable window
d3da6dc4 294 if(fPhotCkov[i]>ckovMax) ckovMax=fPhotCkov[i];
a591e55f 295 weightThetaCerenkov += fPhotCkov[i]*fPhotWei[i];
296 wei += fPhotWei[i]; //collect weight as sum of all candidate weghts
d3da6dc4 297
d3da6dc4 298 sigma2 += 1./Sigma2(fPhotCkov[i],fPhotPhi[i]);
299 }
300 }//candidates loop
301
302 if(sigma2>0) fCkovSigma2=1./sigma2;
303 else fCkovSigma2=1e10;
304
b4ad85e9 305 if(wei != 0.) weightThetaCerenkov /= wei; else weightThetaCerenkov = 0.;
d3da6dc4 306 return weightThetaCerenkov;
307}//FindCkovRing()
308//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
309Int_t AliHMPIDRecon::FlagPhot(Double_t ckov)
310{
311// Flag photon candidates if their individual ckov angle is inside the window around ckov angle returned by HoughResponse()
312// Arguments: ckov- value of most probable ckov angle for track as returned by HoughResponse()
313// Returns: number of photon candidates happened to be inside the window
314
a591e55f 315// Photon Flag: Flag = 0 initial set;
316// Flag = 1 good candidate (charge compatible with photon);
317// Flag = 2 photon used for the ring;
d3da6dc4 318
319 Int_t steps = (Int_t)((ckov )/ fDTheta); //how many times we need to have fDTheta to fill the distance between 0 and thetaCkovHough
320
321 Double_t tmin = (Double_t)(steps - 1)*fDTheta;
322 Double_t tmax = (Double_t)(steps)*fDTheta;
323 Double_t tavg = 0.5*(tmin+tmax);
324
325 tmin = tavg - 0.5*fWindowWidth; tmax = tavg + 0.5*fWindowWidth;
326
327 Int_t iInsideCnt = 0; //count photons which Theta ckov inside the window
328 for(Int_t i=0;i<fPhotCnt;i++){//photon candidates loop
afe12692 329 fPhotFlag[i] = 0;
d3da6dc4 330 if(fPhotCkov[i] >= tmin && fPhotCkov[i] <= tmax) {
331 fPhotFlag[i]=2;
332 iInsideCnt++;
333 }
334 }
335 return iInsideCnt;
336}//FlagPhot()
337//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
a591e55f 338TVector2 AliHMPIDRecon::TracePhot(Double_t ckovThe,Double_t ckovPhi)const
d3da6dc4 339{
340// Trace a single Ckov photon from emission point somewhere in radiator up to photocathode taking into account ref indexes of materials it travereses
ffb1ac19 341// Arguments: ckovThe,ckovPhi- photon ckov angles in TRS, [rad]
d3da6dc4 342// Returns: distance between photon point on PC and track projection
343 TRotation mtheta; mtheta.RotateY(fTrkDir.Theta());
344 TRotation mphi; mphi.RotateZ(fTrkDir.Phi());
345 TRotation mrot=mphi*mtheta;
a591e55f 346 TVector3 dirCkov,dirCkovTors;
347
348 dirCkovTors.SetMagThetaPhi(1,ckovThe,ckovPhi); //initially photon is directed according to requested ckov angle
349 dirCkov=mrot*dirCkovTors; //now we know photon direction in LORS
350 return TraceForward(dirCkov);
351}//TracePhot()
d3da6dc4 352//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
a591e55f 353void AliHMPIDRecon::Propagate(const TVector3 dir,TVector3 &pos,Double_t z)const
d3da6dc4 354{
355// Finds an intersection point between a line and XY plane shifted along Z.
356// Arguments: dir,pos - vector along the line and any point of the line
357// z - z coordinate of plain
358// Returns: none
359// On exit: pos is the position if this intesection if any
360 static TVector3 nrm(0,0,1);
361 TVector3 pnt(0,0,z);
362
363 TVector3 diff=pnt-pos;
364 Double_t sint=(nrm*diff)/(nrm*dir);
365 pos+=sint*dir;
366}//Propagate()
367//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
368void AliHMPIDRecon::Refract(TVector3 &dir,Double_t n1,Double_t n2)const
369{
370// Refract direction vector according to Snell law
371// Arguments:
372// n1 - ref idx of first substance
373// n2 - ref idx of second substance
374// Returns: none
375// On exit: dir is new direction
67a1c24c 376 Double_t sinref=(n1/n2)*TMath::Sin(dir.Theta());
76fd1a96 377 if(TMath::Abs(sinref)>1.) dir.SetXYZ(-999,-999,-999);
67a1c24c 378 else dir.SetTheta(TMath::ASin(sinref));
d3da6dc4 379}//Refract()
380//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
381Double_t AliHMPIDRecon::HoughResponse()
382{
383//
611e810d 384// fIdxMip = mipId;
385
d3da6dc4 386//
387 Double_t kThetaMax=0.75;
388 Int_t nChannels = (Int_t)(kThetaMax/fDTheta+0.5);
389 TH1D *phots = new TH1D("Rphot" ,"phots" ,nChannels,0,kThetaMax);
390 TH1D *photsw = new TH1D("RphotWeighted" ,"photsw" ,nChannels,0,kThetaMax);
391 TH1D *resultw = new TH1D("resultw","resultw" ,nChannels,0,kThetaMax);
392 Int_t nBin = (Int_t)(kThetaMax/fDTheta);
393 Int_t nCorrBand = (Int_t)(fWindowWidth/(2*fDTheta));
394
395 for (Int_t i=0; i< fPhotCnt; i++){//photon cadidates loop
396 Double_t angle = fPhotCkov[i]; if(angle<0||angle>kThetaMax) continue;
397 phots->Fill(angle);
398 Int_t bin = (Int_t)(0.5+angle/(fDTheta));
399 Double_t weight=1.;
400 if(fIsWEIGHT){
afe12692 401 Double_t lowerlimit = ((Double_t)bin)*fDTheta - 0.5*fDTheta; Double_t upperlimit = ((Double_t)bin)*fDTheta + 0.5*fDTheta;
4598109f 402 Double_t diffArea = FindRingArea(upperlimit)-FindRingArea(lowerlimit);
d3da6dc4 403 if(diffArea>0) weight = 1./diffArea;
404 }
405 photsw->Fill(angle,weight);
406 fPhotWei[i]=weight;
407 }//photon candidates loop
408
409 for (Int_t i=1; i<=nBin;i++){
410 Int_t bin1= i-nCorrBand;
411 Int_t bin2= i+nCorrBand;
412 if(bin1<1) bin1=1;
413 if(bin2>nBin)bin2=nBin;
414 Double_t sumPhots=phots->Integral(bin1,bin2);
415 if(sumPhots<3) continue; // if less then 3 photons don't trust to this ring
416 Double_t sumPhotsw=photsw->Integral(bin1,bin2);
417 resultw->Fill((Double_t)((i+0.5)*fDTheta),sumPhotsw);
418 }
419// evaluate the "BEST" theta ckov as the maximum value of histogramm
420 Double_t *pVec = resultw->GetArray();
421 Int_t locMax = TMath::LocMax(nBin,pVec);
3ebd8038 422 delete phots;delete photsw;delete resultw; // Reset and delete objects
d3da6dc4 423
424 return (Double_t)(locMax*fDTheta+0.5*fDTheta); //final most probable track theta ckov
425}//HoughResponse()
426//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
427Double_t AliHMPIDRecon::Sigma2(Double_t ckovTh, Double_t ckovPh)const
428{
429// Analithical calculation of total error (as a sum of localization, geometrical and chromatic errors) on Cerenkov angle for a given Cerenkov photon
430// created by a given MIP. Fromulae according to CERN-EP-2000-058
431// Arguments: Cerenkov and azimuthal angles for Cerenkov photon, [radians]
432// dip and azimuthal angles for MIP taken at the entrance to radiator, [radians]
433// MIP beta
434// Returns: absolute error on Cerenkov angle, [radians]
435
436 TVector3 v(-999,-999,-999);
ffb1ac19 437 Double_t trkBeta = 1./(TMath::Cos(ckovTh)*fParam->GetRefIdx());
d3498bf5 438
439 if(trkBeta > 1) trkBeta = 1; //protection against bad measured thetaCer
440 if(trkBeta < 0) trkBeta = 0.0001; //
d3da6dc4 441
442 v.SetX(SigLoc (ckovTh,ckovPh,trkBeta));
443 v.SetY(SigGeom(ckovTh,ckovPh,trkBeta));
444 v.SetZ(SigCrom(ckovTh,ckovPh,trkBeta));
445
446 return v.Mag2();
447}
448//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
449Double_t AliHMPIDRecon::SigLoc(Double_t thetaC, Double_t phiC,Double_t betaM)const
450{
451// Analithical calculation of localization error (due to finite segmentation of PC) on Cerenkov angle for a given Cerenkov photon
452// created by a given MIP. Fromulae according to CERN-EP-2000-058
453// Arguments: Cerenkov and azimuthal angles for Cerenkov photon, [radians]
454// dip and azimuthal angles for MIP taken at the entrance to radiator, [radians]
455// MIP beta
456// Returns: absolute error on Cerenkov angle, [radians]
d3498bf5 457
d3da6dc4 458 Double_t phiDelta = phiC - fTrkDir.Phi();
459
d3498bf5 460 Double_t sint = TMath::Sin(fTrkDir.Theta());
461 Double_t cost = TMath::Cos(fTrkDir.Theta());
462 Double_t sinf = TMath::Sin(fTrkDir.Phi());
463 Double_t cosf = TMath::Cos(fTrkDir.Phi());
464 Double_t sinfd = TMath::Sin(phiDelta);
465 Double_t cosfd = TMath::Cos(phiDelta);
466 Double_t tantheta = TMath::Tan(thetaC);
467
468 Double_t alpha =cost-tantheta*cosfd*sint; // formula (11)
ffb1ac19 469 Double_t k = 1.-fParam->GetRefIdx()*fParam->GetRefIdx()+alpha*alpha/(betaM*betaM); // formula (after 8 in the text)
d3da6dc4 470 if (k<0) return 1e10;
d3498bf5 471 Double_t mu =sint*sinf+tantheta*(cost*cosfd*sinf+sinfd*cosf); // formula (10)
472 Double_t e =sint*cosf+tantheta*(cost*cosfd*cosf-sinfd*sinf); // formula (9)
d3da6dc4 473
ffb1ac19 474 Double_t kk = betaM*TMath::Sqrt(k)/(fParam->GapThick()*alpha); // formula (6) and (7)
d3498bf5 475 Double_t dtdxc = kk*(k*(cosfd*cosf-cost*sinfd*sinf)-(alpha*mu/(betaM*betaM))*sint*sinfd); // formula (6)
476 Double_t dtdyc = kk*(k*(cosfd*sinf+cost*sinfd*cosf)+(alpha* e/(betaM*betaM))*sint*sinfd); // formula (7) pag.4
d3da6dc4 477
d3498bf5 478 Double_t errX = 0.2,errY=0.25; //end of page 7
479 return TMath::Sqrt(errX*errX*dtdxc*dtdxc + errY*errY*dtdyc*dtdyc);
d3da6dc4 480}
481//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
482Double_t AliHMPIDRecon::SigCrom(Double_t thetaC, Double_t phiC,Double_t betaM)const
483{
484// Analithical calculation of chromatic error (due to lack of knowledge of Cerenkov photon energy) on Cerenkov angle for a given Cerenkov photon
485// created by a given MIP. Fromulae according to CERN-EP-2000-058
486// Arguments: Cerenkov and azimuthal angles for Cerenkov photon, [radians]
487// dip and azimuthal angles for MIP taken at the entrance to radiator, [radians]
488// MIP beta
489// Returns: absolute error on Cerenkov angle, [radians]
d3498bf5 490
d3da6dc4 491 Double_t phiDelta = phiC - fTrkDir.Phi();
d3da6dc4 492
d3498bf5 493 Double_t sint = TMath::Sin(fTrkDir.Theta());
494 Double_t cost = TMath::Cos(fTrkDir.Theta());
495 Double_t cosfd = TMath::Cos(phiDelta);
496 Double_t tantheta = TMath::Tan(thetaC);
497
498 Double_t alpha =cost-tantheta*cosfd*sint; // formula (11)
ffb1ac19 499 Double_t dtdn = cost*fParam->GetRefIdx()*betaM*betaM/(alpha*tantheta); // formula (12)
d3da6dc4 500
d3498bf5 501// Double_t f = 0.00928*(7.75-5.635)/TMath::Sqrt(12.);
502 Double_t f = 0.0172*(7.75-5.635)/TMath::Sqrt(24.);
d3da6dc4 503
504 return f*dtdn;
505}//SigCrom()
506//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
507Double_t AliHMPIDRecon::SigGeom(Double_t thetaC, Double_t phiC,Double_t betaM)const
508{
509// Analithical calculation of geometric error (due to lack of knowledge of creation point in radiator) on Cerenkov angle for a given Cerenkov photon
510// created by a given MIP. Formulae according to CERN-EP-2000-058
511// Arguments: Cerenkov and azimuthal angles for Cerenkov photon, [radians]
512// dip and azimuthal angles for MIP taken at the entrance to radiator, [radians]
513// MIP beta
514// Returns: absolute error on Cerenkov angle, [radians]
515
516 Double_t phiDelta = phiC - fTrkDir.Phi();
d3da6dc4 517
d3498bf5 518 Double_t sint = TMath::Sin(fTrkDir.Theta());
519 Double_t cost = TMath::Cos(fTrkDir.Theta());
520 Double_t sinf = TMath::Sin(fTrkDir.Phi());
521 Double_t cosfd = TMath::Cos(phiDelta);
522 Double_t costheta = TMath::Cos(thetaC);
523 Double_t tantheta = TMath::Tan(thetaC);
524
525 Double_t alpha =cost-tantheta*cosfd*sint; // formula (11)
526
ffb1ac19 527 Double_t k = 1.-fParam->GetRefIdx()*fParam->GetRefIdx()+alpha*alpha/(betaM*betaM); // formula (after 8 in the text)
d3da6dc4 528 if (k<0) return 1e10;
529
ffb1ac19 530 Double_t eTr = 0.5*fParam->RadThick()*betaM*TMath::Sqrt(k)/(fParam->GapThick()*alpha); // formula (14)
d3498bf5 531 Double_t lambda = 1.-sint*sint*sinf*sinf; // formula (15)
d3da6dc4 532
d3498bf5 533 Double_t c1 = 1./(1.+ eTr*k/(alpha*alpha*costheta*costheta)); // formula (13.a)
ffb1ac19 534 Double_t c2 = betaM*TMath::Power(k,1.5)*tantheta*lambda/(fParam->GapThick()*alpha*alpha); // formula (13.b)
d3498bf5 535 Double_t c3 = (1.+eTr*k*betaM*betaM)/((1+eTr)*alpha*alpha); // formula (13.c)
ffb1ac19 536 Double_t c4 = TMath::Sqrt(k)*tantheta*(1-lambda)/(fParam->GapThick()*betaM); // formula (13.d)
d3498bf5 537 Double_t dtdT = c1 * (c2+c3*c4);
ffb1ac19 538 Double_t trErr = fParam->RadThick()/(TMath::Sqrt(12.)*cost);
d3da6dc4 539
d3498bf5 540 return trErr*dtdT;
d3da6dc4 541}//SigGeom()