Coding rule violations
[u/mrichter/AliRoot.git] / HMPID / AliHMPIDReconHTA.cxx
CommitLineData
5a3482a0 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// AliHMPIDReconHTA //
19// //
20// HMPID class to perfom pattern recognition based on Hough transfrom //
21// for single chamber //
22//////////////////////////////////////////////////////////////////////////
23
24#include "AliHMPIDReconHTA.h"//class header
25#include "AliHMPIDCluster.h" //CkovHiddenTrk()
26#include "AliHMPIDRecon.h" //FunMinPhot()
27#include <TMinuit.h> //FitFree()
28#include <TClonesArray.h> //CkovHiddenTrk()
29#include <AliESDtrack.h> //CkovHiddenTrk()
30#include <TH2I.h> //InitDatabase()
31#include <TGraph.h> //ShapeModel()
32#include <TFile.h> //ShapeModel()
33#include <TSpline.h> //ShapeModel()
34#include "TStopwatch.h" //
35
36TH2I* AliHMPIDReconHTA::fDatabase = 0x0;
37
38//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
39AliHMPIDReconHTA::AliHMPIDReconHTA():TTask("RichRec","RichPat")
40{
41//..
42//hidden algorithm
43//..
44 fMipX=fMipY=fThTrkFit=fPhTrkFit=fCkovFit=fMipQ=fRadX=fRadY=-999;
45 fIdxMip=fNClu=0;
46 fCkovSig2=0;
47 fXClu = 0x0;
48 fYClu = 0x0;
49 fClCk = 0x0;
50
51 fParam=AliHMPIDParam::Instance();
52
53 fParam->SetRefIdx(fParam->MeanIdxRad()); // initialization of ref index to a default one
54 if(!fDatabase) InitDatabase();
55}
56//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
57AliHMPIDReconHTA::~AliHMPIDReconHTA()
58{
59 DeleteVars();
60}
61//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
62void AliHMPIDReconHTA::InitVars(Int_t n)
63{
64//..
65//Init some variables
66//..
67 fXClu = new Double_t[n];
68 fYClu = new Double_t[n];
69 fClCk = new Bool_t[n];
70 for(Int_t i=0;i<n;i++) fClCk[i] = kTRUE;
71//
72}
73//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
74void AliHMPIDReconHTA::DeleteVars()
75{
76//..
77//Delete variables
78//..
79 if(fXClu) delete fXClu;
80 if(fYClu) delete fYClu;
81 if(fClCk) delete fClCk;
82}
83//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
84Bool_t AliHMPIDReconHTA::CkovHiddenTrk(AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t nmean, Double_t qthre)
85{
86// Pattern recognition method without any infos from tracking:HTA (Hidden Track Algorithm)...
87// The method finds in the chmuber the cluster with the highest charge
88// compatibile with a MIP, then the strategy is applied
89// Arguments: pTrk - pointer to ESD track
90// pCluLs - list of clusters for a given chamber
91// nmean - mean freon ref. index
92// Returns: - 0=ok,1=not fitted
93
94 AliHMPIDParam *pParam = AliHMPIDParam::Instance();
95 pParam->SetRefIdx(nmean);
96
97 if(!CluPreFilter(pCluLst)) {return kFALSE;}
98
99 Float_t mipX=-1,mipY=-1;Int_t mipId=-1,mipQ=-1;
100 Double_t qRef = 0;
101 Int_t nCh=0;
102 for (Int_t iClu=0;iClu<pCluLst->GetEntriesFast();iClu++){ //clusters loop
103 AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(iClu); //get pointer to current cluster
104 nCh = pClu->Ch();
105 fXClu[iClu] = pClu->X();fYClu[iClu] = pClu->Y(); //store x,y for fitting procedure
106 fClCk[iClu] = kTRUE; //all cluster are accepted at this stage to be reconstructed
107 if(pClu->Q()>qthre) fClCk[iClu] = kFALSE; // not a good photon in any case (multiple MIPs)
108 if(pClu->Q()>qRef){ //searching the highest charge to select a MIP
109 qRef = pClu->Q();
110 mipId=iClu; mipX=pClu->X();mipY=pClu->Y();mipQ=(Int_t)pClu->Q();
111 }
112// Printf(" n. %d x %f y %f Q %f",iClu,pClu->X(),pClu->Y(),pClu->Q());
113 }//clusters loop
114
115 fNClu = pCluLst->GetEntriesFast();
116 if(qRef>qthre){ //charge compartible with MIP clusters
117 fIdxMip = mipId;
118 fClCk[mipId] = kFALSE;
119 fMipX = mipX; fMipY=mipY; fMipQ = qRef;
120// Printf(" mipId %d x %f y %f Q %f",fIdxMip,fMipX,fMipY,fMipQ);
121 if(!DoRecHiddenTrk()) {
122 pTrk->SetHMPIDsignal(pParam->kNoPhotAccept);
123 return kFALSE;
124 } //Do track and ring reconstruction,if problems returns 1
125 pTrk->SetHMPIDtrk(fRadX,fRadY,fThTrkFit,fPhTrkFit); //store track intersection info
126 pTrk->SetHMPIDmip(fMipX,fMipY,(Int_t)fMipQ,fNClu); //store mip info
127 pTrk->SetHMPIDcluIdx(nCh,fIdxMip); //set cham number and index of cluster
128 pTrk->SetHMPIDsignal(fCkovFit); //find best Theta ckov for ring i.e. track
129 pTrk->SetHMPIDchi2(fCkovSig2); //errors squared
130// Printf(" n clusters tot %i accepted %i",pCluLst->GetEntriesFast(),fNClu);
131// Printf("CkovHiddenTrk: thetaC %f th %f ph %f",fCkovFit,fThTrkFit,fPhTrkFit);
132 return kTRUE;
133 }
134
135 return kFALSE;
136}//CkovHiddenTrk()
137//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
138Bool_t AliHMPIDReconHTA::DoRecHiddenTrk()
139{
140// Pattern recognition method without any infos from tracking...
141// First a preclustering filter to avoid part of the noise
142// Then only ellipsed-rings are fitted (no possibility,
143// for the moment, to reconstruct very inclined tracks)
144// Finally a fitting with (th,ph) free, starting by very close values
145// previously evaluated.
146// Arguments: none
147// Returns: none
148 Double_t thTrkRec,phiTrkRec,thetaCRec;
149
150 if(!FindShape(thTrkRec,phiTrkRec,thetaCRec)) return kFALSE;
151
152 if(!FitFree(thTrkRec,phiTrkRec)) {return kFALSE;}
153
154 return kTRUE;
155}//DoRecHiddenTrk()
156/*
157//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
158Bool_t AliHMPIDReconHTA::DoRecHiddenTrk(TClonesArray *pCluLst)
159{
160// Pattern recognition method without any infos from tracking...
161// First a preclustering filter to avoid part of the noise
162// Then only ellipsed-rings are fitted (no possibility,
163// for the moment, to reconstruct very inclined tracks)
164// Finally a fitting with (th,ph) free, starting by very close values
165// previously evaluated.
166// Arguments: none
167// Returns: none
168 Double_t thTrkRec,phiTrkRec,thetaCRec;
169
170 if(!FindShape(thTrkRec,phiTrkRec,thetaCRec)) return kFALSE;
171
172// Printf("thTrkRec %f phiTrkRec %f ThetaCRec %f",thTrkRec*TMath::RadToDeg(),phiTrkRec*TMath::RadToDeg(),thetaCRec*TMath::RadToDeg());
173 Int_t nClTmp1 = pCluLst->GetEntriesFast()-1; //minus MIP...
174 Int_t nClTmp2 = 0;
175
176 while(nClTmp1 != nClTmp2){
177 SetNClu(pCluLst->GetEntriesFast());
178 if(!FitFree(thTrkRec,phiTrkRec)) {return kFALSE;}
179 nClTmp2 = NClu();
180 if(nClTmp2!=nClTmp1) {nClTmp1=nClTmp2;nClTmp2=0;}
181 }
182
183 fNClu = nClTmp2;
184 return kTRUE;
185}//DoRecHiddenTrk()
186*/
187//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
188Bool_t AliHMPIDReconHTA::CluPreFilter(TClonesArray *pCluLst)
189{
190// Filter of bkg clusters
191// based on elliptical-shapes...
192//
193 Int_t nClusTot = pCluLst->GetEntriesFast();
194 if(nClusTot<4||nClusTot>100) {
195 return kFALSE;
196 } else {
197 InitVars(nClusTot);
198 return kTRUE;
199 }
200}
201//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
202Bool_t AliHMPIDReconHTA::FindShape(Double_t &thTrkRec,Double_t &phiTrkRec,Double_t &thetaCRec)
203{
204 Double_t *phiphot = new Double_t[fNClu];
205 Double_t *dist = new Double_t[fNClu];
206 Int_t *indphi = new Int_t[fNClu];
207
208 Bool_t status;
209
210// Sort in phi angle...
211 for(Int_t i=0;i<fNClu;i++) {
212 phiphot[i] = (TMath::ATan2(fMipY-fYClu[i],fMipX-fXClu[i])+TMath::Pi())*TMath::RadToDeg();
213 dist[i]=TMath::Sqrt((fMipX-fXClu[i])*(fMipX-fXClu[i])+(fMipY-fYClu[i])*(fMipY-fYClu[i]));
214 }
215
216 TMath::Sort(fNClu,phiphot,indphi,kFALSE);
217
218// Purify with a truncated mean;
219 Int_t np=0;
220 Double_t dMean = 0;
221 Double_t dMean2 = 0;
222 for(Int_t i=0;i<fNClu;i++) {
223 if(!fClCk[indphi[i]]) continue; // Check if a good photon candidate or not
224 dMean +=dist[indphi[i]];
225 dMean2+=dist[indphi[i]]*dist[indphi[i]];
226 np++;
227 }
228
229 dMean /=(Double_t)np;
230 dMean2 /=(Double_t)np;
231 Double_t rms = TMath::Sqrt(dMean2 - dMean*dMean);
232
233 for(Int_t i=0;i<fNClu;i++) {
234 if(!fClCk[indphi[i]]) continue; // Check if a good photon candidate or not
235 if(TMath::Abs(dMean-dist[indphi[i]]) > 2*rms) {
236 fClCk[indphi[i]] = kFALSE;
237 continue;
238 }
239 }
240
241//
242// purify vectors for good photon candidates
243//
244 Int_t npeff=0;
245 Double_t *phiphotP = new Double_t[fNClu+1];
246 Double_t *distP = new Double_t[fNClu+1];
247 for(Int_t i=0;i<fNClu;i++) {
248 if(!fClCk[indphi[i]]) continue; // Check if a good photon candidate or not
249 phiphotP[npeff] = phiphot[indphi[i]];
250 distP[npeff] = dist[indphi[i]];
251 npeff++;
252 }
253
254 delete [] phiphot;
255 delete [] dist;
256 delete [] indphi;
257
258 if(npeff<3) {
259 delete [] phiphotP;
260 delete [] distP;
261 return kFALSE;
262 }
263
264// for(Int_t i=0;i<npeff;i++) {Printf(" n. %d phiphot %f dist %f",i,phiphotP[i],distP[i]);}
265
266 Double_t xA,xB;
267 if(ShapeModel(npeff,phiphotP,distP,xA,xB,phiTrkRec)) {
268
269// Printf("FindShape: phi start %f",phiTrkRec*TMath::RadToDeg());
270
271 Int_t bin = fDatabase->FindBin(xA,xB);
272 Int_t compact = (Int_t)fDatabase->GetBinContent(bin);
273 thetaCRec = (Double_t)(compact%1000);
274 thTrkRec = (Double_t)(compact/1000);
275
276 thTrkRec *= TMath::DegToRad();
277 thetaCRec *= TMath::DegToRad();
278
279// Printf("FindShape: xA %f xB %f compact %d thTrk %f thC %f",xA,xB,compact,thTrkRec*TMath::RadToDeg(),thetaCRec*TMath::RadToDeg());
280
281 status = kTRUE;
282
283 } else {
284
285 status = kFALSE;
286
287 }
288
289 delete [] phiphotP;
290 delete [] distP;
291
292 return status;
293}
294//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
295Bool_t AliHMPIDReconHTA::ShapeModel(Int_t np,Double_t *phiphot,Double_t *dist,Double_t &xA,Double_t &xB,Double_t &phiStart)
296{
297 TGraph *phigr = new TGraph(np,phiphot,dist);
298 TSpline3 *sphi = new TSpline3("sphi",phigr);
299 if(!sphi) {Printf("Spline not created!Bye.");return kFALSE;}
300
301 Int_t locMin = TMath::LocMin(np,dist);
302 Int_t locMax = TMath::LocMax(np,dist);
303
304 Double_t minX = phiphot[locMin];
305// Double_t minY = dist[locMin];
306 Double_t maxX = phiphot[locMax];
307// Double_t maxY = dist[locMax];
308
309 Int_t ip[3] = {-1,0,1};
310 if(locMin==0 ) {ip[0]= 0;ip[1]= 1;ip[2]= 2;}
311 if(locMin==np-1) {ip[0]=-2;ip[1]=-1;ip[2]= 0;}
312
313 Double_t minXf = VertParab(phiphot[locMin+ip[0]],dist[locMin+ip[0]],
314 phiphot[locMin+ip[1]],dist[locMin+ip[1]],
315 phiphot[locMin+ip[2]],dist[locMin+ip[2]]);
316 if(minXf< phiphot[locMin+ip[0]] || minXf > phiphot[locMin+ip[2]]) minXf = minX;
317
318 ip[0]=-1;ip[1]=0;ip[2]=1;
319 if(locMax==0 ) {ip[0]= 0;ip[1]= 1;ip[2]= 2;}
320 if(locMax==np-1) {ip[0]=-2;ip[1]=-1;ip[2]= 0;}
321
322 Double_t maxXf = VertParab(phiphot[locMax+ip[0]],dist[locMax+ip[0]],
323 phiphot[locMax+ip[1]],dist[locMax+ip[1]],
324 phiphot[locMax+ip[2]],dist[locMax+ip[2]]);
325 if(maxXf< phiphot[locMax+ip[0]] || maxXf > phiphot[locMax+ip[2]]) maxXf = maxX;
326
327// Printf(" phi at mindist %f and found %f",minX,minXf);
328// Printf(" phi at maxdist %f and found %f",maxX,maxXf);
329//
330 if(TMath::Abs(maxXf-minXf)>30) {
331 xA = sphi->Eval(minXf);
332 if(xA < 0) xA = dist[sphi->FindX(xA)];
333 xB = sphi->Eval(minXf-90);
334 if(xB < 0) xB = dist[sphi->FindX(xB)];
335 phiStart = minXf-180; //open ring or acceptance effect...so believe to min phi angle!
336 } else {
337 phiStart = 0.5*(maxXf-180+minXf);
338 xA = sphi->Eval(phiStart+180);
339 if(xA < 0) xA = dist[sphi->FindX(xA)];
340 xB = sphi->Eval(phiStart+90);
341 if(xB < 0) xB = dist[sphi->FindX(xB)];
342 }
343 //
344// Printf("ShapeModel: phiStart %f xA %f xB %f",phiStart,xA,xB);
345
346 phiStart*=TMath::DegToRad();
347 return kTRUE;
348}
349//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
350Double_t AliHMPIDReconHTA::VertParab(Double_t x1,Double_t y1,Double_t x2, Double_t y2, Double_t x3, Double_t y3)
351{
352 Double_t a = ((x1-x3)*(y1-y2)-(x1-x2)*(y1-y3))/((x1*x1-x2*x2)*(x1-x3)-(x1*x1-x3*x3)*(x1-x2));
353 Double_t b = (y1-y2 - a*(x1*x1-x2*x2))/(x1-x2);
354// Double_t c = y1 - a*x1*x1-b*x1;
355 return -0.5*b/a;
356}
357//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
358Bool_t AliHMPIDReconHTA::FitFree(Double_t thTrkRec,Double_t phiTrkRec)
359{
360// Fit performed by minimizing RMS/sqrt(n) of the
361// photons reconstructed. First phi is fixed and theta
362// is fouond, then (th,ph) of the track
363// as free parameters
364// Arguments: PhiRec phi of the track
365// Returns: none
366
367
368 TMinuit *pMinuit = new TMinuit(2);
369 pMinuit->mncler();
370 gMinuit->SetObjectFit((TObject*)this); gMinuit->SetFCN(AliHMPIDReconHTA::FunMinPhot); //set fit function
371 Double_t aArg=-1,parStep,parLow,parHigh; Int_t iErrFlg; //tmp vars for TMinuit
372 Double_t d1,d2,d3;
373 TString sName;
374 Double_t th,ph;
375
376 pMinuit->mnexcm("SET PRI",&aArg,1,iErrFlg); //suspend all printout from TMinuit
377 pMinuit->mnexcm("SET NOW",&aArg,0,iErrFlg);
378
379 if(thTrkRec==0) thTrkRec = 3.*TMath::DegToRad(); // not to start from the edge...
380
381 pMinuit->mnparm(0," thTrk ",thTrkRec ,parStep=0.01,parLow=0,parHigh=TMath::PiOver4(),iErrFlg);
382 pMinuit->mnparm(1," phiTrk ",phiTrkRec,parStep=0.01,parLow=0,parHigh=TMath::TwoPi(),iErrFlg);
383
384 pMinuit->FixParameter(1);
385 pMinuit->mnexcm("SIMPLEX" ,&aArg,0,iErrFlg);
386 pMinuit->mnexcm("MIGRAD" ,&aArg,0,iErrFlg);
387 pMinuit->Release(1);
388 pMinuit->mnexcm("MIGRAD" ,&aArg,0,iErrFlg);
389
390 pMinuit->mnpout(0,sName,th,d1,d2,d3,iErrFlg);
391 pMinuit->mnpout(1,sName,ph,d1,d2,d3,iErrFlg);
392
393 Double_t f,par[2];
394 Double_t *grad=0x0;
395 par[0] = th;par[1] = ph;
396 pMinuit->Eval(2,grad,f,par,3);
397
398// Printf("FitFree: theta %f phi %f",th,ph);
399
400 SetTrkFit(th,ph);
401 return kTRUE;
402}
403//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
404void AliHMPIDReconHTA::FunMinPhot(Int_t &/* */,Double_t* /* */,Double_t &f,Double_t *par,Int_t iflag)
405{
406// Minimization function to find best track and thetaC parameters
407// Arguments: f = function value to minimize
408// par = list of parameter to find
409// iflag = flag status. See Minuit instructions
410// Returns: none
411//
412// Note: it is necessary to call an instance of AlihMPIDParam. Not possible to use fParam
413// because of the static instantiation of the function in Minuit
414
415 AliHMPIDParam *pParam=AliHMPIDParam::Instance();
416 AliHMPIDReconHTA *pRecHTA=(AliHMPIDReconHTA*)gMinuit->GetObjectFit();
417 AliHMPIDRecon pRec;
418 Double_t sizeCh = 0.5*pParam->RadThick()+pParam->WinThick()+pParam->GapThick();
419 Double_t thTrk = par[0];
420 Double_t phTrk = par[1];
421 Double_t xrad = pRecHTA->MipX() - sizeCh*TMath::Tan(thTrk)*TMath::Cos(phTrk);
422 Double_t yrad = pRecHTA->MipY() - sizeCh*TMath::Tan(thTrk)*TMath::Sin(phTrk);
423 pRecHTA->SetRadXY(xrad,yrad);
424 pRec.SetTrack(xrad,yrad,thTrk,phTrk);
425
426 Double_t meanCkov =0;
427 Double_t meanCkov2=0;
428 Double_t thetaCer,phiCer;
429 Int_t nClAcc = 0;
430 Int_t nClTot=pRecHTA->NClu();
431
432 for(Int_t i=0;i<nClTot;i++) {
433 if(!(pRecHTA->ClCk(i))) continue;
434 pRec.FindPhotCkov(pRecHTA->XClu(i),pRecHTA->YClu(i),thetaCer,phiCer);
435 meanCkov += thetaCer;
436 meanCkov2 += thetaCer*thetaCer;
437 nClAcc++;
438 }
439 if(nClAcc==0) {f=999;return;}
440 meanCkov /=(Double_t)nClAcc;
441 meanCkov2 /=(Double_t)nClAcc;
442 Double_t rms = TMath::Sqrt(TMath::Abs(meanCkov2 - meanCkov*meanCkov));
443 f = rms/TMath::Sqrt((Double_t)nClAcc);
444
445 if(iflag==3) {
446/*
447 Printf("FunMinPhot before: photons candidates %i used %i",nClTot,nClAcc);
448 nClAcc = 0;
449 Double_t meanCkov1=0;
450 Double_t meanCkov2=0;
451 for(Int_t i=0;i<nClTot;i++) {
452 if(!(pRec->ClCk(i))) continue;
453 pRec->FindPhotCkov(pRec->XClu(i),pRec->YClu(i),thetaCer,phiCer);
454 if(TMath::Abs(thetaCer-meanCkov)<2*rms) {
455 meanCkov1 += thetaCer;
456 meanCkov2 += thetaCer*thetaCer;
457 nClAcc++;
458 } else pRec->SetClCk(i,kFALSE);
459 }
460 meanCkov1/=nClAcc;
461 Double_t rms2 = (meanCkov2 - meanCkov*meanCkov*nClAcc)/nClAcc;
462 Printf("FunMinPhot after: photons candidates %i used %i thetaCer %f",nClTot,nClAcc,meanCkov1);
463 pRec->SetCkovFit(meanCkov1);
464 pRec->SetCkovSig2(rms2);
465 pRec->SetNClu(nClAcc);
466*/
467// Printf("FunMinPhot: reconstructed theta Cerenkov %f with %d photons",meanCkov,nClAcc);
468 pRecHTA->SetCkovFit(meanCkov);
469 pRecHTA->SetCkovSig2(rms*rms);
470 pRecHTA->SetNClu(nClAcc);
471 }
472}//FunMinPhot()
473//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
474void AliHMPIDReconHTA::InitDatabase()
475{
476 TStopwatch timer;
477 timer.Start();
478
479 AliInfo(Form("database HTA is being built.Please, wait..."));
480//
481 Double_t x[3],y[3];
482
483 AliHMPIDRecon rec;
484
485 if(!fParam) fParam=AliHMPIDParam::Instance();
486 Double_t thetaMax = TMath::ACos(1./fParam->MeanIdxRad());
487 Double_t thTrkMax = 1./TMath::ASin(fParam->MeanIdxRad());
488
489 Int_t nstepx = 1000;
490 Int_t nstepy = 1000;
491
492 TH2I *deconv = new TH2I("deconv","database;d1;d2;thC+1000*thTrk",500,0,50,150,0,15);
493 //
494 Double_t xrad = 0;
495 Double_t yrad = 0;
496 Double_t phTrk = 0;
497
498 for(Int_t i=0;i<nstepx;i++) { //loop on thetaC
499 for(Int_t j=0;j<nstepy;j++) { //loop on theta particle
500 Double_t thetaC = thetaMax/nstepx*((Double_t)i+0.5);
501 Double_t thTrk = thTrkMax/nstepy*((Double_t)j+0.5);
502 //
503 //mip position
504 //
505 Double_t sizeCh = fParam->RadThick()+fParam->WinThick()+fParam->GapThick();
506 Double_t xmip = xrad + sizeCh*TMath::Tan(thTrk)*TMath::Cos(phTrk);
507 Double_t ymip = yrad + sizeCh*TMath::Tan(thTrk)*TMath::Sin(phTrk);
508
509 Double_t dist1,dist2;
510 //
511 //first point at phi=0
512 //
513 rec.SetTrack(xrad,yrad,thTrk,phTrk);
514 TVector2 pos;
515 pos=rec.TracePhot(thetaC,0);
516
517 if(pos.X()==-999) {
518 dist1 = 0; //open ring...anly the distance btw mip and point at 180 will be considered
519 } else {
520 x[0] = pos.X(); y[0] = pos.Y();
521 dist1 = TMath::Sqrt((x[0]-xmip)*(x[0]-xmip)+(y[0]-ymip)*(y[0]-ymip));
522 }
523 //
524 //second point at phi=180
525 //
526 rec.SetTrack(xrad,yrad,thTrk,phTrk);
527 pos=rec.TracePhot(thetaC,TMath::Pi());
528
529 if(pos.X()==-999) {Printf("it should not happens!Bye");return;}
530 x[1] = pos.X(); y[1] = pos.Y();
531 if((x[1]-xmip)*(x[0]-xmip)>0) continue; // to avoid circles out mips (for very low ThetaC)
532 dist2 = TMath::Sqrt((x[1]-xmip)*(x[1]-xmip)+(y[1]-ymip)*(y[1]-ymip));
533
534// Double_t distA = dist1+dist2;
535 Double_t distA = dist2; // only the minimum: problem of acceptance
536 //
537 //second point at phi=90
538 //
539 rec.SetTrack(xrad,yrad,thTrk,phTrk);
540 pos=rec.TracePhot(thetaC,TMath::PiOver2());
541
542 if(pos.X()==-999) continue;
543 x[2] = pos.X(); y[2] = pos.Y();
544 Double_t distB = TMath::Sqrt((x[2]-xmip)*(x[2]-xmip)+(y[2]-ymip)*(y[2]-ymip));
545// compact the infos...
546 Int_t compact = (Int_t)(thetaC*TMath::RadToDeg())+1000*(Int_t)(thTrk*TMath::RadToDeg());
547 Int_t bin = deconv->FindBin(distA,distB);
548 if(deconv->GetBinContent(bin)==0) deconv->Fill(distA,distB,compact);
549 }
550 }
551
552 FillZeroChan(deconv);
553 fDatabase = deconv;
554
555 timer.Stop();
556 Double_t nSecs = timer.CpuTime();
557 AliInfo(Form("database HTA successfully open in %3.1f sec.(CPU). Reconstruction is started.",nSecs));
558}
559//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
560void AliHMPIDReconHTA::FillZeroChan(TH2I *deconv)
561{
562 Int_t nbinx = deconv->GetNbinsX();
563 Int_t nbiny = deconv->GetNbinsY();
564 for(Int_t i = 0;i<nbinx;i++) {
565 for(Int_t j = 0;j<nbiny;j++) {
566 if(deconv->GetBinContent(i,j) == 0) {
567 Int_t nXmin = i-1; Int_t nXmax=i+1;
568 Int_t nYmin = j-1; Int_t nYmax=j+1;
569 Int_t nc = 0;
570 Double_t meanC =0;
571 Double_t meanTrk =0;
572 for(Int_t ix=nXmin;ix<=nXmax;ix++) {
573 if(ix<0||ix>nbinx) continue;
574 for(Int_t iy=nYmin;iy<=nYmax;iy++) {
575 if(iy<0||iy>nbiny) continue;
576 meanC += (Int_t)deconv->GetBinContent(ix,iy)%1000;
577 meanTrk+= (Int_t)deconv->GetBinContent(ix,iy)/1000;
578 nc++;
579 }
580 meanC/=nc; meanTrk/=nc;
581 Int_t compact = (Int_t)meanC+1000*(Int_t)meanTrk;
582 if(compact>0)deconv->SetCellContent(i,j,compact);
583 }
584 }
585 }
586 }
587}
588//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++