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