]> 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
Database as TH2F removed + minors.
[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()
9785d5fb 27#include <TFile.h> //Database()
5a3482a0 28#include <TMinuit.h> //FitFree()
29#include <TClonesArray.h> //CkovHiddenTrk()
30#include <AliESDtrack.h> //CkovHiddenTrk()
9785d5fb 31#include <TH2F.h> //InitDatabase()
5a3482a0 32#include <TGraph.h> //ShapeModel()
5a3482a0 33#include <TSpline.h> //ShapeModel()
34#include "TStopwatch.h" //
35
170a4fc5 36Int_t AliHMPIDReconHTA::fgDB[500][150]={75000*0};
5a3482a0 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),
ee2f3539 50 fThTrkIn(-999),
51 fPhTrkIn(-999),
c61a7285 52 fThTrkFit(-999),
53 fPhTrkFit(-999),
54 fCkovFit(-999),
55 fCkovSig2(0),
56 fParam(AliHMPIDParam::Instance())
5a3482a0 57{
58//..
59//hidden algorithm
60//..
5a3482a0 61 fParam->SetRefIdx(fParam->MeanIdxRad()); // initialization of ref index to a default one
5565f017 62 InitDatabase();
5a3482a0 63}
64//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
65AliHMPIDReconHTA::~AliHMPIDReconHTA()
66{
67 DeleteVars();
68}
69//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
70void AliHMPIDReconHTA::InitVars(Int_t n)
71{
72//..
73//Init some variables
74//..
75 fXClu = new Double_t[n];
76 fYClu = new Double_t[n];
77 fClCk = new Bool_t[n];
78 for(Int_t i=0;i<n;i++) fClCk[i] = kTRUE;
79//
80}
81//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
aa00f952 82void AliHMPIDReconHTA::DeleteVars()const
5a3482a0 83{
84//..
85//Delete variables
86//..
87 if(fXClu) delete fXClu;
88 if(fYClu) delete fYClu;
89 if(fClCk) delete fClCk;
90}
91//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
9785d5fb 92Bool_t AliHMPIDReconHTA::CkovHiddenTrk(AliESDtrack *pTrk,TClonesArray *pCluLst,Int_t index, Double_t nmean)
5a3482a0 93{
94// Pattern recognition method without any infos from tracking:HTA (Hidden Track Algorithm)...
95// The method finds in the chmuber the cluster with the highest charge
96// compatibile with a MIP, then the strategy is applied
97// Arguments: pTrk - pointer to ESD track
98// pCluLs - list of clusters for a given chamber
9785d5fb 99// pNmean - pointer to ref. index
100// pQthre - pointer to qthre
5a3482a0 101// Returns: - 0=ok,1=not fitted
102
103 AliHMPIDParam *pParam = AliHMPIDParam::Instance();
5a3482a0 104
9785d5fb 105 if(!CluPreFilter(pCluLst)) return kFALSE;
106
5a3482a0 107 Int_t nCh=0;
43d3333b 108 Int_t sizeClu=0;
9785d5fb 109
110 fNClu = pCluLst->GetEntriesFast();
111
112 for (Int_t iClu=0;iClu<fNClu;iClu++){ //clusters loop
5a3482a0 113 AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(iClu); //get pointer to current cluster
5a3482a0 114 fXClu[iClu] = pClu->X();fYClu[iClu] = pClu->Y(); //store x,y for fitting procedure
115 fClCk[iClu] = kTRUE; //all cluster are accepted at this stage to be reconstructed
9785d5fb 116
117 if(iClu == index) {
118
119 fMipX = pClu->X();
120 fMipY = pClu->Y();
121 fMipQ = pClu->Q();
122 sizeClu = pClu->Size();
123 nCh = pClu->Ch();
124 fClCk[index] = kFALSE;
125 fIdxMip = index;
126 // Printf(" n. %d x %f y %f Q %f",iClu,pClu->X(),pClu->Y(),pClu->Q());
127 }
5a3482a0 128 }//clusters loop
9785d5fb 129
130 pParam->SetRefIdx(nmean);
131
e56b695f 132 //
133 Float_t xra,yra,th,ph; pTrk->GetHMPIDtrk(xra,yra,th,ph);
ee2f3539 134// Printf(" simulated phi %6.2f ",ph*TMath::RadToDeg());
e56b695f 135 //
136
9785d5fb 137 if(!DoRecHiddenTrk()) {
138 pTrk->SetHMPIDsignal(pParam->kNoPhotAccept);
139 return kFALSE;
140 } //Do track and ring reconstruction,if problems returns 1
ee2f3539 141// Printf(" fitted phi %6.2f ",fPhTrkFit*TMath::RadToDeg());
9785d5fb 142
143 pTrk->SetHMPIDtrk(fRadX,fRadY,fThTrkFit,fPhTrkFit); //store track intersection info
144 pTrk->SetHMPIDmip(fMipX,fMipY,(Int_t)fMipQ,fNClu); //store mip info
145 pTrk->SetHMPIDcluIdx(nCh,fIdxMip+1000*sizeClu); //set cham number, index of cluster + cluster size
146 pTrk->SetHMPIDsignal(fCkovFit); //find best Theta ckov for ring i.e. track
147 pTrk->SetHMPIDchi2(fCkovSig2); //errors squared
5a3482a0 148// Printf(" n clusters tot %i accepted %i",pCluLst->GetEntriesFast(),fNClu);
149// Printf("CkovHiddenTrk: thetaC %f th %f ph %f",fCkovFit,fThTrkFit,fPhTrkFit);
5a3482a0 150
9785d5fb 151 return kTRUE;
152
5a3482a0 153}//CkovHiddenTrk()
154//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
155Bool_t AliHMPIDReconHTA::DoRecHiddenTrk()
156{
157// Pattern recognition method without any infos from tracking...
158// First a preclustering filter to avoid part of the noise
159// Then only ellipsed-rings are fitted (no possibility,
160// for the moment, to reconstruct very inclined tracks)
161// Finally a fitting with (th,ph) free, starting by very close values
162// previously evaluated.
163// Arguments: none
164// Returns: none
165 Double_t thTrkRec,phiTrkRec,thetaCRec;
166
ee2f3539 167 if(!FindShape(thTrkRec,phiTrkRec,thetaCRec)) {
168// Printf("FindShape failed...!");
169 return kFALSE;
170 }
5a3482a0 171
ee2f3539 172 if(!FitFree(thTrkRec,phiTrkRec)) {
173// Printf("FitFree failed...!");
174 return kFALSE;
175 }
176
5a3482a0 177 return kTRUE;
178}//DoRecHiddenTrk()
5a3482a0 179//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
180Bool_t AliHMPIDReconHTA::CluPreFilter(TClonesArray *pCluLst)
181{
aa00f952 182// Pre-filter of bkg clusters
183// Arguments: pSluLst - List of the clusters for a given chamber
184// Returns: status - TRUE if filtering leaves enough photons, FALSE if not
5a3482a0 185//
186 Int_t nClusTot = pCluLst->GetEntriesFast();
187 if(nClusTot<4||nClusTot>100) {
188 return kFALSE;
189 } else {
190 InitVars(nClusTot);
191 return kTRUE;
192 }
193}
194//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
195Bool_t AliHMPIDReconHTA::FindShape(Double_t &thTrkRec,Double_t &phiTrkRec,Double_t &thetaCRec)
196{
aa00f952 197// Finds the estimates for phi and theta of the track and the ThetaCerenkov
198// by using a database of the shapes of the rings
199// Arguments: none
200// Returns: thTrkRec - estimate of theta track
201// phiTrkRec - estimate of phi track
202// thetaCRec - estimate of ThetaCerenkov
203// status - TRUE if a good solution is found, FALSE if not
204
5a3482a0 205 Double_t *phiphot = new Double_t[fNClu];
206 Double_t *dist = new Double_t[fNClu];
207 Int_t *indphi = new Int_t[fNClu];
208
209 Bool_t status;
210
9785d5fb 211// Sort in phi angle...
212// Printf(" mipX %f mipy %f",fMipX,fMipY);
5a3482a0 213 for(Int_t i=0;i<fNClu;i++) {
9785d5fb 214 if(!fClCk[i]) {
215 phiphot[i] = 999.;
216 dist[i] = 999.;
217 continue;
218 }
5a3482a0 219 phiphot[i] = (TMath::ATan2(fMipY-fYClu[i],fMipX-fXClu[i])+TMath::Pi())*TMath::RadToDeg();
220 dist[i]=TMath::Sqrt((fMipX-fXClu[i])*(fMipX-fXClu[i])+(fMipY-fYClu[i])*(fMipY-fYClu[i]));
9785d5fb 221// Printf(" n.%3i phiphot %f dist %f check %i",i,phiphot[i],dist[i],fClCk[i]);
5a3482a0 222 }
223
224 TMath::Sort(fNClu,phiphot,indphi,kFALSE);
225
226// Purify with a truncated mean;
227 Int_t np=0;
228 Double_t dMean = 0;
229 Double_t dMean2 = 0;
230 for(Int_t i=0;i<fNClu;i++) {
231 if(!fClCk[indphi[i]]) continue; // Check if a good photon candidate or not
232 dMean +=dist[indphi[i]];
233 dMean2+=dist[indphi[i]]*dist[indphi[i]];
234 np++;
235 }
236
237 dMean /=(Double_t)np;
238 dMean2 /=(Double_t)np;
239 Double_t rms = TMath::Sqrt(dMean2 - dMean*dMean);
240
241 for(Int_t i=0;i<fNClu;i++) {
242 if(!fClCk[indphi[i]]) continue; // Check if a good photon candidate or not
243 if(TMath::Abs(dMean-dist[indphi[i]]) > 2*rms) {
244 fClCk[indphi[i]] = kFALSE;
245 continue;
246 }
247 }
248
249//
250// purify vectors for good photon candidates
251//
252 Int_t npeff=0;
253 Double_t *phiphotP = new Double_t[fNClu+1];
254 Double_t *distP = new Double_t[fNClu+1];
255 for(Int_t i=0;i<fNClu;i++) {
256 if(!fClCk[indphi[i]]) continue; // Check if a good photon candidate or not
257 phiphotP[npeff] = phiphot[indphi[i]];
258 distP[npeff] = dist[indphi[i]];
9785d5fb 259// Printf("n. %2i phi %f dist %f",npeff,phiphotP[npeff],distP[npeff]);
5a3482a0 260 npeff++;
261 }
262
263 delete [] phiphot;
264 delete [] dist;
265 delete [] indphi;
266
267 if(npeff<3) {
ee2f3539 268 Printf("FindShape failed: no enough photons = %i...",npeff);
5a3482a0 269 delete [] phiphotP;
270 delete [] distP;
271 return kFALSE;
272 }
273
274// for(Int_t i=0;i<npeff;i++) {Printf(" n. %d phiphot %f dist %f",i,phiphotP[i],distP[i]);}
275
276 Double_t xA,xB;
9785d5fb 277 status = kFALSE;
ee2f3539 278
170a4fc5 279 if(!ShapeModel(npeff,phiphotP,distP,xA,xB,phiTrkRec)) {/*Printf("ShapeModel failed ");*/ return kFALSE;}
5a3482a0 280
ee2f3539 281// if(xA > 50 || xB > 15) {Printf("xA and xB failed out of range"); return kFALSE;}
5a3482a0 282
170a4fc5 283 Int_t binxDB,binyDB;
284 FindBinDB(xA,xB,binxDB,binyDB);
285 if(binxDB<0 || binyDB<0) {/*Printf("bin < 0 ! failed ");*/ return kFALSE;}
ee2f3539 286
170a4fc5 287 Int_t compactDB = fgDB[binxDB][binyDB];
ee2f3539 288
170a4fc5 289 if(compactDB<0) {/*Printf("compact< 0! failed ");*/ return kFALSE;}
ee2f3539 290 //
ee2f3539 291 //
170a4fc5 292 thetaCRec = (Double_t)(compactDB%1000);
293 thTrkRec = (Double_t)(compactDB/1000);
5a3482a0 294
ee2f3539 295 thTrkRec *= TMath::DegToRad();
296 thetaCRec *= TMath::DegToRad();
9785d5fb 297
ee2f3539 298 status = kTRUE;
5a3482a0 299
300 delete [] phiphotP;
301 delete [] distP;
302
303 return status;
304}
305//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
306Bool_t AliHMPIDReconHTA::ShapeModel(Int_t np,Double_t *phiphot,Double_t *dist,Double_t &xA,Double_t &xB,Double_t &phiStart)
307{
aa00f952 308// Find a Spline curve to define dist. vs. phi angle
309// in order to estimate the phi of the track
310// Arguments: np - # points corresponding to # photon candidates
311// dist - distance of each photon from MIP
312// phiphot - phi of the photon in the DRS
313// Returns: xA - min. distance from MIP
314// xB - dist. from mip perpedicular to the major axis
315// phiStart- estimate of the track phi
316
5a3482a0 317 TGraph *phigr = new TGraph(np,phiphot,dist);
318 TSpline3 *sphi = new TSpline3("sphi",phigr);
319 if(!sphi) {Printf("Spline not created!Bye.");return kFALSE;}
320
321 Int_t locMin = TMath::LocMin(np,dist);
322 Int_t locMax = TMath::LocMax(np,dist);
323
324 Double_t minX = phiphot[locMin];
325// Double_t minY = dist[locMin];
326 Double_t maxX = phiphot[locMax];
327// Double_t maxY = dist[locMax];
328
329 Int_t ip[3] = {-1,0,1};
330 if(locMin==0 ) {ip[0]= 0;ip[1]= 1;ip[2]= 2;}
331 if(locMin==np-1) {ip[0]=-2;ip[1]=-1;ip[2]= 0;}
332
333 Double_t minXf = VertParab(phiphot[locMin+ip[0]],dist[locMin+ip[0]],
334 phiphot[locMin+ip[1]],dist[locMin+ip[1]],
335 phiphot[locMin+ip[2]],dist[locMin+ip[2]]);
336 if(minXf< phiphot[locMin+ip[0]] || minXf > phiphot[locMin+ip[2]]) minXf = minX;
337
338 ip[0]=-1;ip[1]=0;ip[2]=1;
339 if(locMax==0 ) {ip[0]= 0;ip[1]= 1;ip[2]= 2;}
340 if(locMax==np-1) {ip[0]=-2;ip[1]=-1;ip[2]= 0;}
341
342 Double_t maxXf = VertParab(phiphot[locMax+ip[0]],dist[locMax+ip[0]],
343 phiphot[locMax+ip[1]],dist[locMax+ip[1]],
344 phiphot[locMax+ip[2]],dist[locMax+ip[2]]);
345 if(maxXf< phiphot[locMax+ip[0]] || maxXf > phiphot[locMax+ip[2]]) maxXf = maxX;
346
5a3482a0 347//
348 if(TMath::Abs(maxXf-minXf)>30) {
349 xA = sphi->Eval(minXf);
350 if(xA < 0) xA = dist[sphi->FindX(xA)];
351 xB = sphi->Eval(minXf-90);
352 if(xB < 0) xB = dist[sphi->FindX(xB)];
353 phiStart = minXf-180; //open ring or acceptance effect...so believe to min phi angle!
354 } else {
355 phiStart = 0.5*(maxXf-180+minXf);
356 xA = sphi->Eval(phiStart+180);
357 if(xA < 0) xA = dist[sphi->FindX(xA)];
358 xB = sphi->Eval(phiStart+90);
359 if(xB < 0) xB = dist[sphi->FindX(xB)];
360 }
361 //
362// Printf("ShapeModel: phiStart %f xA %f xB %f",phiStart,xA,xB);
363
9a573d52 364
e56b695f 365 phiStart*=TMath::DegToRad();
ee2f3539 366 //----
367 Double_t phitest = FindSimmPhi();
368 phiStart = phitest;
369 //---
e56b695f 370// Printf("phiStart %f phiTest %f",phiStart*TMath::RadToDeg(),phitest*TMath::RadToDeg());
9785d5fb 371
5a3482a0 372 return kTRUE;
373}
374//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
aa00f952 375Double_t AliHMPIDReconHTA::VertParab(Double_t x1,Double_t y1,Double_t x2, Double_t y2, Double_t x3, Double_t y3)const
5a3482a0 376{
aa00f952 377// It uses parabola from 3 points to evaluate the x-coord of the parab
378// Arguments: xi,yi - points
379// Returns: x-coord of the vertex
380
5a3482a0 381 Double_t a = ((x1-x3)*(y1-y2)-(x1-x2)*(y1-y3))/((x1*x1-x2*x2)*(x1-x3)-(x1*x1-x3*x3)*(x1-x2));
382 Double_t b = (y1-y2 - a*(x1*x1-x2*x2))/(x1-x2);
383// Double_t c = y1 - a*x1*x1-b*x1;
384 return -0.5*b/a;
385}
386//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
387Bool_t AliHMPIDReconHTA::FitFree(Double_t thTrkRec,Double_t phiTrkRec)
388{
389// Fit performed by minimizing RMS/sqrt(n) of the
390// photons reconstructed. First phi is fixed and theta
391// is fouond, then (th,ph) of the track
392// as free parameters
393// Arguments: PhiRec phi of the track
394// Returns: none
395
5a3482a0 396 TMinuit *pMinuit = new TMinuit(2);
397 pMinuit->mncler();
398 gMinuit->SetObjectFit((TObject*)this); gMinuit->SetFCN(AliHMPIDReconHTA::FunMinPhot); //set fit function
399 Double_t aArg=-1,parStep,parLow,parHigh; Int_t iErrFlg; //tmp vars for TMinuit
400 Double_t d1,d2,d3;
401 TString sName;
402 Double_t th,ph;
403
404 pMinuit->mnexcm("SET PRI",&aArg,1,iErrFlg); //suspend all printout from TMinuit
405 pMinuit->mnexcm("SET NOW",&aArg,0,iErrFlg);
406
407 if(thTrkRec==0) thTrkRec = 3.*TMath::DegToRad(); // not to start from the edge...
408
409 pMinuit->mnparm(0," thTrk ",thTrkRec ,parStep=0.01,parLow=0,parHigh=TMath::PiOver4(),iErrFlg);
410 pMinuit->mnparm(1," phiTrk ",phiTrkRec,parStep=0.01,parLow=0,parHigh=TMath::TwoPi(),iErrFlg);
411
412 pMinuit->FixParameter(1);
413 pMinuit->mnexcm("SIMPLEX" ,&aArg,0,iErrFlg);
414 pMinuit->mnexcm("MIGRAD" ,&aArg,0,iErrFlg);
415 pMinuit->Release(1);
416 pMinuit->mnexcm("MIGRAD" ,&aArg,0,iErrFlg);
417
418 pMinuit->mnpout(0,sName,th,d1,d2,d3,iErrFlg);
419 pMinuit->mnpout(1,sName,ph,d1,d2,d3,iErrFlg);
420
421 Double_t f,par[2];
422 Double_t *grad=0x0;
423 par[0] = th;par[1] = ph;
424 pMinuit->Eval(2,grad,f,par,3);
425
5a3482a0 426 SetTrkFit(th,ph);
427 return kTRUE;
428}
429//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
430void AliHMPIDReconHTA::FunMinPhot(Int_t &/* */,Double_t* /* */,Double_t &f,Double_t *par,Int_t iflag)
431{
432// Minimization function to find best track and thetaC parameters
433// Arguments: f = function value to minimize
434// par = list of parameter to find
435// iflag = flag status. See Minuit instructions
436// Returns: none
437//
438// Note: it is necessary to call an instance of AlihMPIDParam. Not possible to use fParam
439// because of the static instantiation of the function in Minuit
440
441 AliHMPIDParam *pParam=AliHMPIDParam::Instance();
442 AliHMPIDReconHTA *pRecHTA=(AliHMPIDReconHTA*)gMinuit->GetObjectFit();
443 AliHMPIDRecon pRec;
444 Double_t sizeCh = 0.5*pParam->RadThick()+pParam->WinThick()+pParam->GapThick();
445 Double_t thTrk = par[0];
446 Double_t phTrk = par[1];
447 Double_t xrad = pRecHTA->MipX() - sizeCh*TMath::Tan(thTrk)*TMath::Cos(phTrk);
448 Double_t yrad = pRecHTA->MipY() - sizeCh*TMath::Tan(thTrk)*TMath::Sin(phTrk);
449 pRecHTA->SetRadXY(xrad,yrad);
450 pRec.SetTrack(xrad,yrad,thTrk,phTrk);
451
452 Double_t meanCkov =0;
453 Double_t meanCkov2=0;
454 Double_t thetaCer,phiCer;
455 Int_t nClAcc = 0;
456 Int_t nClTot=pRecHTA->NClu();
457
458 for(Int_t i=0;i<nClTot;i++) {
459 if(!(pRecHTA->ClCk(i))) continue;
460 pRec.FindPhotCkov(pRecHTA->XClu(i),pRecHTA->YClu(i),thetaCer,phiCer);
461 meanCkov += thetaCer;
462 meanCkov2 += thetaCer*thetaCer;
463 nClAcc++;
464 }
465 if(nClAcc==0) {f=999;return;}
466 meanCkov /=(Double_t)nClAcc;
467 meanCkov2 /=(Double_t)nClAcc;
468 Double_t rms = TMath::Sqrt(TMath::Abs(meanCkov2 - meanCkov*meanCkov));
469 f = rms/TMath::Sqrt((Double_t)nClAcc);
470
471 if(iflag==3) {
472/*
473 Printf("FunMinPhot before: photons candidates %i used %i",nClTot,nClAcc);
474 nClAcc = 0;
475 Double_t meanCkov1=0;
476 Double_t meanCkov2=0;
477 for(Int_t i=0;i<nClTot;i++) {
478 if(!(pRec->ClCk(i))) continue;
479 pRec->FindPhotCkov(pRec->XClu(i),pRec->YClu(i),thetaCer,phiCer);
480 if(TMath::Abs(thetaCer-meanCkov)<2*rms) {
481 meanCkov1 += thetaCer;
482 meanCkov2 += thetaCer*thetaCer;
483 nClAcc++;
484 } else pRec->SetClCk(i,kFALSE);
485 }
486 meanCkov1/=nClAcc;
487 Double_t rms2 = (meanCkov2 - meanCkov*meanCkov*nClAcc)/nClAcc;
488 Printf("FunMinPhot after: photons candidates %i used %i thetaCer %f",nClTot,nClAcc,meanCkov1);
489 pRec->SetCkovFit(meanCkov1);
490 pRec->SetCkovSig2(rms2);
491 pRec->SetNClu(nClAcc);
492*/
493// Printf("FunMinPhot: reconstructed theta Cerenkov %f with %d photons",meanCkov,nClAcc);
494 pRecHTA->SetCkovFit(meanCkov);
495 pRecHTA->SetCkovSig2(rms*rms);
496 pRecHTA->SetNClu(nClAcc);
497 }
498}//FunMinPhot()
499//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
500void AliHMPIDReconHTA::InitDatabase()
501{
aa00f952 502// Construction a database of ring shapes on fly
503// Arguments: none
504// Returns : none
170a4fc5 505// N.B. fgDB is the distance with x-min from MIP
506// y-dist from the ring of the MIP perpendicular to major axis
aa00f952 507// The content is the packed info of track theta and thetaC in degrees
508// thetaC+1000*thTrk
509//
9785d5fb 510// TFile *pout = new TFile("./database.root","recreate");
ee2f3539 511
170a4fc5 512 static Bool_t isDone = kFALSE;
513
5a3482a0 514 TStopwatch timer;
170a4fc5 515
516 if(isDone) {
5565f017 517 return;
518 }
170a4fc5 519
520 if(!isDone) {
521 timer.Start();
522 }
523
5a3482a0 524 AliInfo(Form("database HTA is being built.Please, wait..."));
525//
c61a7285 526 Double_t x[3]={0,0,0},y[3];
5a3482a0 527
528 AliHMPIDRecon rec;
529
530 if(!fParam) fParam=AliHMPIDParam::Instance();
531 Double_t thetaMax = TMath::ACos(1./fParam->MeanIdxRad());
532 Double_t thTrkMax = 1./TMath::ASin(fParam->MeanIdxRad());
533
534 Int_t nstepx = 1000;
535 Int_t nstepy = 1000;
536
5a3482a0 537 //
538 Double_t xrad = 0;
539 Double_t yrad = 0;
540 Double_t phTrk = 0;
541
542 for(Int_t i=0;i<nstepx;i++) { //loop on thetaC
543 for(Int_t j=0;j<nstepy;j++) { //loop on theta particle
544 Double_t thetaC = thetaMax/nstepx*((Double_t)i+0.5);
545 Double_t thTrk = thTrkMax/nstepy*((Double_t)j+0.5);
546 //
547 //mip position
548 //
fcaff63d 549 Double_t sizeCh = 0.5*fParam->RadThick()+fParam->WinThick()+fParam->GapThick();
5a3482a0 550 Double_t xmip = xrad + sizeCh*TMath::Tan(thTrk)*TMath::Cos(phTrk);
551 Double_t ymip = yrad + sizeCh*TMath::Tan(thTrk)*TMath::Sin(phTrk);
552
553 Double_t dist1,dist2;
554 //
555 //first point at phi=0
556 //
557 rec.SetTrack(xrad,yrad,thTrk,phTrk);
558 TVector2 pos;
559 pos=rec.TracePhot(thetaC,0);
560
561 if(pos.X()==-999) {
ee2f3539 562 dist1 = 0; //open ring...only the distance btw mip and point at 180 will be considered
5a3482a0 563 } else {
564 x[0] = pos.X(); y[0] = pos.Y();
565 dist1 = TMath::Sqrt((x[0]-xmip)*(x[0]-xmip)+(y[0]-ymip)*(y[0]-ymip));
566 }
567 //
568 //second point at phi=180
569 //
570 rec.SetTrack(xrad,yrad,thTrk,phTrk);
571 pos=rec.TracePhot(thetaC,TMath::Pi());
572
573 if(pos.X()==-999) {Printf("it should not happens!Bye");return;}
574 x[1] = pos.X(); y[1] = pos.Y();
575 if((x[1]-xmip)*(x[0]-xmip)>0) continue; // to avoid circles out mips (for very low ThetaC)
576 dist2 = TMath::Sqrt((x[1]-xmip)*(x[1]-xmip)+(y[1]-ymip)*(y[1]-ymip));
577
578// Double_t distA = dist1+dist2;
579 Double_t distA = dist2; // only the minimum: problem of acceptance
580 //
581 //second point at phi=90
582 //
583 rec.SetTrack(xrad,yrad,thTrk,phTrk);
584 pos=rec.TracePhot(thetaC,TMath::PiOver2());
585
586 if(pos.X()==-999) continue;
587 x[2] = pos.X(); y[2] = pos.Y();
588 Double_t distB = TMath::Sqrt((x[2]-xmip)*(x[2]-xmip)+(y[2]-ymip)*(y[2]-ymip));
589// compact the infos...
590 Int_t compact = (Int_t)(thetaC*TMath::RadToDeg())+1000*(Int_t)(thTrk*TMath::RadToDeg());
ee2f3539 591 Int_t binxDB,binyDB;
592 FindBinDB(distA,distB,binxDB,binyDB);
170a4fc5 593 if(fgDB[binxDB][binyDB]==0) fgDB[binxDB][binyDB] = compact;
5a3482a0 594 }
595 }
596
9785d5fb 597 FillZeroChan();
5a3482a0 598
170a4fc5 599 if(!isDone) {
600 timer.Stop();
601 Double_t nSecs = timer.CpuTime();
602 AliInfo(Form("database HTA successfully open in %3.1f sec.(CPU). Reconstruction is started.",nSecs));
603 isDone = kTRUE;
604 }
9785d5fb 605
606// pout->Write();
607// pout->Close();
608
9a573d52 609}//InitDatabase()
5a3482a0 610//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
9785d5fb 611void AliHMPIDReconHTA::FillZeroChan()const
5a3482a0 612{
606697a8 613 //If fills eventually channel without entries
614 //inthe histo "database" jyst interpolating the neighboring cells
615 // Arguments: histogram pointer of the database
616 // Returns: none
617 //
170a4fc5 618 const Int_t nxDB = 500;
619 const Int_t nyDB = 150;
620
621 for(Int_t i = 0;i<nxDB;i++) {
622 for(Int_t j = 0;j<nyDB;j++) {
623 if(fgDB[i][j] == 0) {
624 fgDB[i][j] = -1;
5a3482a0 625 Int_t nXmin = i-1; Int_t nXmax=i+1;
626 Int_t nYmin = j-1; Int_t nYmax=j+1;
627 Int_t nc = 0;
628 Double_t meanC =0;
629 Double_t meanTrk =0;
630 for(Int_t ix=nXmin;ix<=nXmax;ix++) {
170a4fc5 631 if(ix<0||ix>=nxDB) continue;
5a3482a0 632 for(Int_t iy=nYmin;iy<=nYmax;iy++) {
170a4fc5 633 if(iy<0||iy>=nyDB) continue;
634 meanC += (Int_t)(fgDB[ix][iy]%1000);
635 meanTrk+= (Int_t)(fgDB[ix][iy]/1000);
5a3482a0 636 nc++;
637 }
638 meanC/=nc; meanTrk/=nc;
639 Int_t compact = (Int_t)meanC+1000*(Int_t)meanTrk;
170a4fc5 640 if(compact>0)fgDB[i][j] = compact;
5a3482a0 641 }
642 }
643 }
644 }
9a573d52 645}//FillZeroChan()
5a3482a0 646//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
ee2f3539 647Int_t AliHMPIDReconHTA::r2(Double_t *coef, Double_t &x1, Double_t &x2)
648{
649 //2nd deg. equation
650 //solution
651 // Arguments: coef 2 1 0: ax^2+bx+c=0
652 // Returns: n. of solutions
653 // x1= 1st sol
654 // x2= 2nd sol
655 Double_t a,b,c;
656 a = coef[2];
657 b = coef[1];
658 c = coef[0];
659 Double_t delta = b*b-4*a*c;
660 if(delta<0) {return 0;}
661 if(delta==0) {
662 x1=x2=-b/(2*a);
663 return 1;
664 }
665 // delta>0
666 x1 = (-b+TMath::Sqrt(delta))/(2*a);
667 x2 = (-b-TMath::Sqrt(delta))/(2*a);
668 return 2;
669}//r2()
670//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
9a573d52 671
672Double_t AliHMPIDReconHTA::FindSimmPhi()
ee2f3539 673{
674//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
675// RESTITUISCE IN OUTPUT IL VALORE FINALE DELL'ANGOLO RICOSTRUITO
676// CON I DUE METODI:
677// - metodo dei minimi quadrati con le distanze effettive;................................(PER RING CHIUSI)
678// - metodo della determin della pendenza individuando la distanza minima mip-fotone;.....(PER RING APERTI)
679
680 Float_t coeff1ord=0; Float_t coeff2ord=0; Float_t coeff0ord=0;
9a573d52 681 Float_t xrotsumm =0; Float_t yrotsumm =0; Float_t xx =0;
ee2f3539 682 Float_t yy =0; Float_t xy =0; Float_t yx =0;
683 Double_t xmin=0;
684 Double_t ymin=0;
9a573d52 685
686 Int_t np=0;
ee2f3539 687
688 Double_t distMin = 999.;
689
9a573d52 690 for(Int_t i=0;i<fNClu;i++) {
691 if(!fClCk[i]) continue;
692 np++;
693 xrotsumm+=fXClu[i]; // summ xi
694 yrotsumm+=fYClu[i]; // summ yi
695 xx+=fXClu[i]*fXClu[i]; // summ xixi
696 yy+=fYClu[i]*fYClu[i]; // summ yiyi
697 xy+=fXClu[i]*fYClu[i]; // summ yixi
ee2f3539 698 Double_t dist2= (fXClu[i]-fMipX)*(fXClu[i]-fMipX)+(fYClu[i]-fMipY)*(fYClu[i]-fMipY);
699 if(dist2<distMin) {
700 distMin = dist2;
701 xmin = fXClu[i];
702 ymin = fYClu[i];
703 }
9a573d52 704 }
ee2f3539 705
706 Double_t AngM=0;
707 if(ymin < fMipY && xmin > fMipX) {AngM = TMath::ATan((ymin-fMipY)/(xmin-fMipX))*TMath::RadToDeg()+180;}
708 if(ymin > fMipY && xmin < fMipX) {AngM = TMath::ATan((ymin-fMipY)/(xmin-fMipX))*TMath::RadToDeg()+360;}
709 if(ymin > fMipY && xmin > fMipX) {AngM = TMath::ATan((ymin-fMipY)/(xmin-fMipX))*TMath::RadToDeg()+180;}
710 if(ymin < fMipY && xmin < fMipX) {AngM = TMath::ATan((ymin-fMipY)/(xmin-fMipX))*TMath::RadToDeg();}
711 if(ymin == fMipY && xmin > fMipX) {AngM = TMath::ATan((ymin-fMipY)/(xmin-fMipX))*TMath::RadToDeg()+180;}
712 if(ymin == fMipY && xmin < fMipX) {AngM = TMath::ATan((ymin-fMipY)/(xmin-fMipX))*TMath::RadToDeg();}
713 if(ymin < fMipY && xmin == fMipX) {AngM = 90;}
714 if(ymin > fMipY && xmin == fMipX) {AngM = 270;}
9a573d52 715
ee2f3539 716 //_____calc. met min quadr using effective distance _________________________________________________
9a573d52 717
ee2f3539 718 coeff2ord= xy-xrotsumm*yrotsumm/np;
719 coeff1ord= yrotsumm*yrotsumm/np - xrotsumm*xrotsumm/np - yy + xx;
720 coeff0ord= xrotsumm*yrotsumm/np - yx;
721 Double_t m1=0, m2=0; Double_t n1=0, n2=0;
722 // c // b // a
723 Double_t coeff[3]={coeff0ord,coeff1ord,coeff2ord};
9a573d52 724
725 r2(coeff,m1,m2);
726
727 n1=(yrotsumm-m1*xrotsumm)/np;
728 n2=(yrotsumm-m2*xrotsumm)/np;
ee2f3539 729 // 2 solutions.................
730
731 Double_t PhiTrk1= TMath::ATan(m1);
4df7ec8c 732 //Double_t PhiTrk2= TMath::ATan(m2);
9a573d52 733
ee2f3539 734 // negative angles solved...
735
736 Double_t PhiTrk1Positive=0;
737
738 if(PhiTrk1<0) PhiTrk1Positive= PhiTrk1 + 180*TMath::DegToRad();
739 if(PhiTrk1>=0) PhiTrk1Positive= PhiTrk1;
9a573d52 740
741 Double_t d1 =(1/(m1*m1+1))*(yy+m1*m1*xx+np*n1*n1-2*m1*xy-2*n1*yrotsumm+2*m1*n1*xrotsumm);
742 Double_t d2 =(1/(m2*m2+1))*(yy+m2*m2*xx+np*n2*n2-2*m2*xy-2*n2*yrotsumm+2*m2*n2*xrotsumm);
ee2f3539 743
9a573d52 744 Double_t mMin;
745 if(d1 > d2) mMin = m2; else mMin = m1;
746
4df7ec8c 747 //Double_t PhiTrk = TMath::ATan(mMin)*TMath::RadToDeg();
ee2f3539 748 Double_t PhiTrkPositive=0;
749 //
750 if(ymin < fMipY && xmin > fMipX) {PhiTrkPositive = TMath::ATan(mMin)*TMath::RadToDeg()+180;}
751 if(ymin > fMipY && xmin < fMipX) {PhiTrkPositive = TMath::ATan(mMin)*TMath::RadToDeg()+360;}
752 if(ymin > fMipY && xmin > fMipX) {PhiTrkPositive = TMath::ATan(mMin)*TMath::RadToDeg()+180;}
753 if(ymin < fMipY && xmin < fMipX) {PhiTrkPositive = TMath::ATan(mMin)*TMath::RadToDeg();}
754 if(ymin == fMipY && xmin > fMipX) {PhiTrkPositive = TMath::ATan(mMin)*TMath::RadToDeg()+180;}
755 if(ymin == fMipY && xmin < fMipX) {PhiTrkPositive = TMath::ATan(mMin)*TMath::RadToDeg();}
756 if(ymin < fMipY && xmin == fMipX) {PhiTrkPositive = 90;}
757 if(ymin > fMipY && xmin == fMipX) {PhiTrkPositive = 270;}
9a573d52 758
ee2f3539 759 // ------------------------- choose the best-----------------------
9a573d52 760
ee2f3539 761
762 Double_t PhiTrkFinal=0;
763 if( AngM-40 <= PhiTrkPositive && AngM+40 >= PhiTrkPositive) PhiTrkFinal = PhiTrkPositive; else PhiTrkFinal = AngM;
764
765 return PhiTrkFinal*TMath::DegToRad();
9a573d52 766}
9a573d52 767//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
ee2f3539 768void AliHMPIDReconHTA::FindBinDB(Double_t x,Double_t y,Int_t &binX,Int_t &binY)
9a573d52 769{
170a4fc5 770 const Int_t nxDB = 500;
771 const Int_t nyDB = 150;
ee2f3539 772 const Double_t xlowDB = 0;
773 const Double_t xhigDB = 50;
774 const Double_t ylowDB = 0;
775 const Double_t yhigDB = 15;
776
777 binX = -1;
778 binY = -1;
170a4fc5 779 if(x<xlowDB && x>xhigDB &&
780 y<ylowDB && y>yhigDB) return;
781 binX = Int_t((x-xlowDB)/(xhigDB-xlowDB)*nxDB);
782 binY = Int_t((y-ylowDB)/(yhigDB-ylowDB)*nyDB);
783 if(binX>=nxDB || binY>=nyDB) {
784 binX = -1;
785 binY = -1;
786 }
787
ee2f3539 788}