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