]> 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
Fixing linking problems for the dHLTdumpraw binary.
[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
9785d5fb 36TH2F* AliHMPIDReconHTA::fgDatabase = new TH2F("deconv","database;d1;d2;thC+1000*thTrk",500,0,50,150,0,15);
37
5a3482a0 38
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),
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
9785d5fb 62 if(fgDatabase->GetEntries()<1) 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
132 if(!DoRecHiddenTrk()) {
133 pTrk->SetHMPIDsignal(pParam->kNoPhotAccept);
134 return kFALSE;
135 } //Do track and ring reconstruction,if problems returns 1
136
137 pTrk->SetHMPIDtrk(fRadX,fRadY,fThTrkFit,fPhTrkFit); //store track intersection info
138 pTrk->SetHMPIDmip(fMipX,fMipY,(Int_t)fMipQ,fNClu); //store mip info
139 pTrk->SetHMPIDcluIdx(nCh,fIdxMip+1000*sizeClu); //set cham number, index of cluster + cluster size
140 pTrk->SetHMPIDsignal(fCkovFit); //find best Theta ckov for ring i.e. track
141 pTrk->SetHMPIDchi2(fCkovSig2); //errors squared
5a3482a0 142// Printf(" n clusters tot %i accepted %i",pCluLst->GetEntriesFast(),fNClu);
143// Printf("CkovHiddenTrk: thetaC %f th %f ph %f",fCkovFit,fThTrkFit,fPhTrkFit);
5a3482a0 144
9785d5fb 145 return kTRUE;
146
5a3482a0 147}//CkovHiddenTrk()
148//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
149Bool_t AliHMPIDReconHTA::DoRecHiddenTrk()
150{
151// Pattern recognition method without any infos from tracking...
152// First a preclustering filter to avoid part of the noise
153// Then only ellipsed-rings are fitted (no possibility,
154// for the moment, to reconstruct very inclined tracks)
155// Finally a fitting with (th,ph) free, starting by very close values
156// previously evaluated.
157// Arguments: none
158// Returns: none
159 Double_t thTrkRec,phiTrkRec,thetaCRec;
160
161 if(!FindShape(thTrkRec,phiTrkRec,thetaCRec)) return kFALSE;
162
163 if(!FitFree(thTrkRec,phiTrkRec)) {return kFALSE;}
164
165 return kTRUE;
166}//DoRecHiddenTrk()
5a3482a0 167//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
168Bool_t AliHMPIDReconHTA::CluPreFilter(TClonesArray *pCluLst)
169{
aa00f952 170// Pre-filter of bkg clusters
171// Arguments: pSluLst - List of the clusters for a given chamber
172// Returns: status - TRUE if filtering leaves enough photons, FALSE if not
5a3482a0 173//
174 Int_t nClusTot = pCluLst->GetEntriesFast();
175 if(nClusTot<4||nClusTot>100) {
176 return kFALSE;
177 } else {
178 InitVars(nClusTot);
179 return kTRUE;
180 }
181}
182//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
183Bool_t AliHMPIDReconHTA::FindShape(Double_t &thTrkRec,Double_t &phiTrkRec,Double_t &thetaCRec)
184{
aa00f952 185// Finds the estimates for phi and theta of the track and the ThetaCerenkov
186// by using a database of the shapes of the rings
187// Arguments: none
188// Returns: thTrkRec - estimate of theta track
189// phiTrkRec - estimate of phi track
190// thetaCRec - estimate of ThetaCerenkov
191// status - TRUE if a good solution is found, FALSE if not
192
5a3482a0 193 Double_t *phiphot = new Double_t[fNClu];
194 Double_t *dist = new Double_t[fNClu];
195 Int_t *indphi = new Int_t[fNClu];
196
197 Bool_t status;
198
9785d5fb 199// Sort in phi angle...
200// Printf(" mipX %f mipy %f",fMipX,fMipY);
5a3482a0 201 for(Int_t i=0;i<fNClu;i++) {
9785d5fb 202 if(!fClCk[i]) {
203 phiphot[i] = 999.;
204 dist[i] = 999.;
205 continue;
206 }
5a3482a0 207 phiphot[i] = (TMath::ATan2(fMipY-fYClu[i],fMipX-fXClu[i])+TMath::Pi())*TMath::RadToDeg();
208 dist[i]=TMath::Sqrt((fMipX-fXClu[i])*(fMipX-fXClu[i])+(fMipY-fYClu[i])*(fMipY-fYClu[i]));
9785d5fb 209// Printf(" n.%3i phiphot %f dist %f check %i",i,phiphot[i],dist[i],fClCk[i]);
5a3482a0 210 }
211
212 TMath::Sort(fNClu,phiphot,indphi,kFALSE);
213
214// Purify with a truncated mean;
215 Int_t np=0;
216 Double_t dMean = 0;
217 Double_t dMean2 = 0;
218 for(Int_t i=0;i<fNClu;i++) {
219 if(!fClCk[indphi[i]]) continue; // Check if a good photon candidate or not
220 dMean +=dist[indphi[i]];
221 dMean2+=dist[indphi[i]]*dist[indphi[i]];
222 np++;
223 }
224
225 dMean /=(Double_t)np;
226 dMean2 /=(Double_t)np;
227 Double_t rms = TMath::Sqrt(dMean2 - dMean*dMean);
228
229 for(Int_t i=0;i<fNClu;i++) {
230 if(!fClCk[indphi[i]]) continue; // Check if a good photon candidate or not
231 if(TMath::Abs(dMean-dist[indphi[i]]) > 2*rms) {
232 fClCk[indphi[i]] = kFALSE;
233 continue;
234 }
235 }
236
237//
238// purify vectors for good photon candidates
239//
240 Int_t npeff=0;
241 Double_t *phiphotP = new Double_t[fNClu+1];
242 Double_t *distP = new Double_t[fNClu+1];
243 for(Int_t i=0;i<fNClu;i++) {
244 if(!fClCk[indphi[i]]) continue; // Check if a good photon candidate or not
245 phiphotP[npeff] = phiphot[indphi[i]];
246 distP[npeff] = dist[indphi[i]];
9785d5fb 247// Printf("n. %2i phi %f dist %f",npeff,phiphotP[npeff],distP[npeff]);
5a3482a0 248 npeff++;
249 }
250
251 delete [] phiphot;
252 delete [] dist;
253 delete [] indphi;
254
255 if(npeff<3) {
256 delete [] phiphotP;
257 delete [] distP;
258 return kFALSE;
259 }
260
261// for(Int_t i=0;i<npeff;i++) {Printf(" n. %d phiphot %f dist %f",i,phiphotP[i],distP[i]);}
262
263 Double_t xA,xB;
9785d5fb 264 status = kFALSE;
5a3482a0 265 if(ShapeModel(npeff,phiphotP,distP,xA,xB,phiTrkRec)) {
266
9785d5fb 267// Printf("FindShape: phi start %f xA %f yA %f",phiTrkRec*TMath::RadToDeg(),xA,xB);
268 if(xA < 50 && xB < 15) { // limits of the Database. See TH2F in InitDatabase()
5a3482a0 269
9785d5fb 270 Int_t bin = fgDatabase->FindBin(xA,xB);
271 if(bin>0) {
272 Int_t compact = (Int_t)fgDatabase->GetBinContent(bin);
273 thetaCRec = (Double_t)(compact%1000);
274 thTrkRec = (Double_t)(compact/1000);
5a3482a0 275
9785d5fb 276 thTrkRec *= TMath::DegToRad();
277 thetaCRec *= TMath::DegToRad();
5a3482a0 278
9785d5fb 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 }
5a3482a0 284 }
285
286 delete [] phiphotP;
287 delete [] distP;
288
289 return status;
290}
291//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
292Bool_t AliHMPIDReconHTA::ShapeModel(Int_t np,Double_t *phiphot,Double_t *dist,Double_t &xA,Double_t &xB,Double_t &phiStart)
293{
aa00f952 294// Find a Spline curve to define dist. vs. phi angle
295// in order to estimate the phi of the track
296// Arguments: np - # points corresponding to # photon candidates
297// dist - distance of each photon from MIP
298// phiphot - phi of the photon in the DRS
299// Returns: xA - min. distance from MIP
300// xB - dist. from mip perpedicular to the major axis
301// phiStart- estimate of the track phi
302
5a3482a0 303 TGraph *phigr = new TGraph(np,phiphot,dist);
304 TSpline3 *sphi = new TSpline3("sphi",phigr);
305 if(!sphi) {Printf("Spline not created!Bye.");return kFALSE;}
306
307 Int_t locMin = TMath::LocMin(np,dist);
308 Int_t locMax = TMath::LocMax(np,dist);
309
310 Double_t minX = phiphot[locMin];
311// Double_t minY = dist[locMin];
312 Double_t maxX = phiphot[locMax];
313// Double_t maxY = dist[locMax];
314
315 Int_t ip[3] = {-1,0,1};
316 if(locMin==0 ) {ip[0]= 0;ip[1]= 1;ip[2]= 2;}
317 if(locMin==np-1) {ip[0]=-2;ip[1]=-1;ip[2]= 0;}
318
319 Double_t minXf = VertParab(phiphot[locMin+ip[0]],dist[locMin+ip[0]],
320 phiphot[locMin+ip[1]],dist[locMin+ip[1]],
321 phiphot[locMin+ip[2]],dist[locMin+ip[2]]);
322 if(minXf< phiphot[locMin+ip[0]] || minXf > phiphot[locMin+ip[2]]) minXf = minX;
323
324 ip[0]=-1;ip[1]=0;ip[2]=1;
325 if(locMax==0 ) {ip[0]= 0;ip[1]= 1;ip[2]= 2;}
326 if(locMax==np-1) {ip[0]=-2;ip[1]=-1;ip[2]= 0;}
327
328 Double_t maxXf = VertParab(phiphot[locMax+ip[0]],dist[locMax+ip[0]],
329 phiphot[locMax+ip[1]],dist[locMax+ip[1]],
330 phiphot[locMax+ip[2]],dist[locMax+ip[2]]);
331 if(maxXf< phiphot[locMax+ip[0]] || maxXf > phiphot[locMax+ip[2]]) maxXf = maxX;
332
5a3482a0 333//
334 if(TMath::Abs(maxXf-minXf)>30) {
335 xA = sphi->Eval(minXf);
336 if(xA < 0) xA = dist[sphi->FindX(xA)];
337 xB = sphi->Eval(minXf-90);
338 if(xB < 0) xB = dist[sphi->FindX(xB)];
339 phiStart = minXf-180; //open ring or acceptance effect...so believe to min phi angle!
340 } else {
341 phiStart = 0.5*(maxXf-180+minXf);
342 xA = sphi->Eval(phiStart+180);
343 if(xA < 0) xA = dist[sphi->FindX(xA)];
344 xB = sphi->Eval(phiStart+90);
345 if(xB < 0) xB = dist[sphi->FindX(xB)];
346 }
347 //
348// Printf("ShapeModel: phiStart %f xA %f xB %f",phiStart,xA,xB);
349
350 phiStart*=TMath::DegToRad();
9a573d52 351
352 Double_t phitest = FindSimmPhi();
353
354 Printf("phiStart %f phiTest %f",phiStart*TMath::RadToDeg(),phitest*TMath::RadToDeg());
9785d5fb 355
5a3482a0 356 return kTRUE;
357}
358//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
aa00f952 359Double_t AliHMPIDReconHTA::VertParab(Double_t x1,Double_t y1,Double_t x2, Double_t y2, Double_t x3, Double_t y3)const
5a3482a0 360{
aa00f952 361// It uses parabola from 3 points to evaluate the x-coord of the parab
362// Arguments: xi,yi - points
363// Returns: x-coord of the vertex
364
5a3482a0 365 Double_t a = ((x1-x3)*(y1-y2)-(x1-x2)*(y1-y3))/((x1*x1-x2*x2)*(x1-x3)-(x1*x1-x3*x3)*(x1-x2));
366 Double_t b = (y1-y2 - a*(x1*x1-x2*x2))/(x1-x2);
367// Double_t c = y1 - a*x1*x1-b*x1;
368 return -0.5*b/a;
369}
370//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
371Bool_t AliHMPIDReconHTA::FitFree(Double_t thTrkRec,Double_t phiTrkRec)
372{
373// Fit performed by minimizing RMS/sqrt(n) of the
374// photons reconstructed. First phi is fixed and theta
375// is fouond, then (th,ph) of the track
376// as free parameters
377// Arguments: PhiRec phi of the track
378// Returns: none
379
5a3482a0 380 TMinuit *pMinuit = new TMinuit(2);
381 pMinuit->mncler();
382 gMinuit->SetObjectFit((TObject*)this); gMinuit->SetFCN(AliHMPIDReconHTA::FunMinPhot); //set fit function
383 Double_t aArg=-1,parStep,parLow,parHigh; Int_t iErrFlg; //tmp vars for TMinuit
384 Double_t d1,d2,d3;
385 TString sName;
386 Double_t th,ph;
387
388 pMinuit->mnexcm("SET PRI",&aArg,1,iErrFlg); //suspend all printout from TMinuit
389 pMinuit->mnexcm("SET NOW",&aArg,0,iErrFlg);
390
391 if(thTrkRec==0) thTrkRec = 3.*TMath::DegToRad(); // not to start from the edge...
392
393 pMinuit->mnparm(0," thTrk ",thTrkRec ,parStep=0.01,parLow=0,parHigh=TMath::PiOver4(),iErrFlg);
394 pMinuit->mnparm(1," phiTrk ",phiTrkRec,parStep=0.01,parLow=0,parHigh=TMath::TwoPi(),iErrFlg);
395
396 pMinuit->FixParameter(1);
397 pMinuit->mnexcm("SIMPLEX" ,&aArg,0,iErrFlg);
398 pMinuit->mnexcm("MIGRAD" ,&aArg,0,iErrFlg);
399 pMinuit->Release(1);
400 pMinuit->mnexcm("MIGRAD" ,&aArg,0,iErrFlg);
401
402 pMinuit->mnpout(0,sName,th,d1,d2,d3,iErrFlg);
403 pMinuit->mnpout(1,sName,ph,d1,d2,d3,iErrFlg);
404
405 Double_t f,par[2];
406 Double_t *grad=0x0;
407 par[0] = th;par[1] = ph;
408 pMinuit->Eval(2,grad,f,par,3);
409
410// Printf("FitFree: theta %f phi %f",th,ph);
411
412 SetTrkFit(th,ph);
413 return kTRUE;
414}
415//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
416void AliHMPIDReconHTA::FunMinPhot(Int_t &/* */,Double_t* /* */,Double_t &f,Double_t *par,Int_t iflag)
417{
418// Minimization function to find best track and thetaC parameters
419// Arguments: f = function value to minimize
420// par = list of parameter to find
421// iflag = flag status. See Minuit instructions
422// Returns: none
423//
424// Note: it is necessary to call an instance of AlihMPIDParam. Not possible to use fParam
425// because of the static instantiation of the function in Minuit
426
427 AliHMPIDParam *pParam=AliHMPIDParam::Instance();
428 AliHMPIDReconHTA *pRecHTA=(AliHMPIDReconHTA*)gMinuit->GetObjectFit();
429 AliHMPIDRecon pRec;
430 Double_t sizeCh = 0.5*pParam->RadThick()+pParam->WinThick()+pParam->GapThick();
431 Double_t thTrk = par[0];
432 Double_t phTrk = par[1];
433 Double_t xrad = pRecHTA->MipX() - sizeCh*TMath::Tan(thTrk)*TMath::Cos(phTrk);
434 Double_t yrad = pRecHTA->MipY() - sizeCh*TMath::Tan(thTrk)*TMath::Sin(phTrk);
435 pRecHTA->SetRadXY(xrad,yrad);
436 pRec.SetTrack(xrad,yrad,thTrk,phTrk);
437
438 Double_t meanCkov =0;
439 Double_t meanCkov2=0;
440 Double_t thetaCer,phiCer;
441 Int_t nClAcc = 0;
442 Int_t nClTot=pRecHTA->NClu();
443
444 for(Int_t i=0;i<nClTot;i++) {
445 if(!(pRecHTA->ClCk(i))) continue;
446 pRec.FindPhotCkov(pRecHTA->XClu(i),pRecHTA->YClu(i),thetaCer,phiCer);
447 meanCkov += thetaCer;
448 meanCkov2 += thetaCer*thetaCer;
449 nClAcc++;
450 }
451 if(nClAcc==0) {f=999;return;}
452 meanCkov /=(Double_t)nClAcc;
453 meanCkov2 /=(Double_t)nClAcc;
454 Double_t rms = TMath::Sqrt(TMath::Abs(meanCkov2 - meanCkov*meanCkov));
455 f = rms/TMath::Sqrt((Double_t)nClAcc);
456
457 if(iflag==3) {
458/*
459 Printf("FunMinPhot before: photons candidates %i used %i",nClTot,nClAcc);
460 nClAcc = 0;
461 Double_t meanCkov1=0;
462 Double_t meanCkov2=0;
463 for(Int_t i=0;i<nClTot;i++) {
464 if(!(pRec->ClCk(i))) continue;
465 pRec->FindPhotCkov(pRec->XClu(i),pRec->YClu(i),thetaCer,phiCer);
466 if(TMath::Abs(thetaCer-meanCkov)<2*rms) {
467 meanCkov1 += thetaCer;
468 meanCkov2 += thetaCer*thetaCer;
469 nClAcc++;
470 } else pRec->SetClCk(i,kFALSE);
471 }
472 meanCkov1/=nClAcc;
473 Double_t rms2 = (meanCkov2 - meanCkov*meanCkov*nClAcc)/nClAcc;
474 Printf("FunMinPhot after: photons candidates %i used %i thetaCer %f",nClTot,nClAcc,meanCkov1);
475 pRec->SetCkovFit(meanCkov1);
476 pRec->SetCkovSig2(rms2);
477 pRec->SetNClu(nClAcc);
478*/
479// Printf("FunMinPhot: reconstructed theta Cerenkov %f with %d photons",meanCkov,nClAcc);
480 pRecHTA->SetCkovFit(meanCkov);
481 pRecHTA->SetCkovSig2(rms*rms);
482 pRecHTA->SetNClu(nClAcc);
483 }
484}//FunMinPhot()
485//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
486void AliHMPIDReconHTA::InitDatabase()
487{
aa00f952 488// Construction a database of ring shapes on fly
489// Arguments: none
490// Returns : none
491// N.B. fgDatabase points to a TH2I with x-min dist from MIP
492// y-dist from the ring of the MIP perpendicular to major axis
493// The content is the packed info of track theta and thetaC in degrees
494// thetaC+1000*thTrk
495//
9785d5fb 496// TFile *pout = new TFile("./database.root","recreate");
497
5a3482a0 498 TStopwatch timer;
499 timer.Start();
500
501 AliInfo(Form("database HTA is being built.Please, wait..."));
502//
c61a7285 503 Double_t x[3]={0,0,0},y[3];
5a3482a0 504
505 AliHMPIDRecon rec;
506
507 if(!fParam) fParam=AliHMPIDParam::Instance();
508 Double_t thetaMax = TMath::ACos(1./fParam->MeanIdxRad());
509 Double_t thTrkMax = 1./TMath::ASin(fParam->MeanIdxRad());
510
511 Int_t nstepx = 1000;
512 Int_t nstepy = 1000;
513
9785d5fb 514// TH2F *fgDatabase = new TH2F("deconv","database;d1;d2;thC+1000*thTrk",500,0,50,150,0,15);
5a3482a0 515 //
516 Double_t xrad = 0;
517 Double_t yrad = 0;
518 Double_t phTrk = 0;
519
520 for(Int_t i=0;i<nstepx;i++) { //loop on thetaC
521 for(Int_t j=0;j<nstepy;j++) { //loop on theta particle
522 Double_t thetaC = thetaMax/nstepx*((Double_t)i+0.5);
523 Double_t thTrk = thTrkMax/nstepy*((Double_t)j+0.5);
524 //
525 //mip position
526 //
fcaff63d 527 Double_t sizeCh = 0.5*fParam->RadThick()+fParam->WinThick()+fParam->GapThick();
5a3482a0 528 Double_t xmip = xrad + sizeCh*TMath::Tan(thTrk)*TMath::Cos(phTrk);
529 Double_t ymip = yrad + sizeCh*TMath::Tan(thTrk)*TMath::Sin(phTrk);
530
531 Double_t dist1,dist2;
532 //
533 //first point at phi=0
534 //
535 rec.SetTrack(xrad,yrad,thTrk,phTrk);
536 TVector2 pos;
537 pos=rec.TracePhot(thetaC,0);
538
539 if(pos.X()==-999) {
540 dist1 = 0; //open ring...anly the distance btw mip and point at 180 will be considered
541 } else {
542 x[0] = pos.X(); y[0] = pos.Y();
543 dist1 = TMath::Sqrt((x[0]-xmip)*(x[0]-xmip)+(y[0]-ymip)*(y[0]-ymip));
544 }
545 //
546 //second point at phi=180
547 //
548 rec.SetTrack(xrad,yrad,thTrk,phTrk);
549 pos=rec.TracePhot(thetaC,TMath::Pi());
550
551 if(pos.X()==-999) {Printf("it should not happens!Bye");return;}
552 x[1] = pos.X(); y[1] = pos.Y();
553 if((x[1]-xmip)*(x[0]-xmip)>0) continue; // to avoid circles out mips (for very low ThetaC)
554 dist2 = TMath::Sqrt((x[1]-xmip)*(x[1]-xmip)+(y[1]-ymip)*(y[1]-ymip));
555
556// Double_t distA = dist1+dist2;
557 Double_t distA = dist2; // only the minimum: problem of acceptance
558 //
559 //second point at phi=90
560 //
561 rec.SetTrack(xrad,yrad,thTrk,phTrk);
562 pos=rec.TracePhot(thetaC,TMath::PiOver2());
563
564 if(pos.X()==-999) continue;
565 x[2] = pos.X(); y[2] = pos.Y();
566 Double_t distB = TMath::Sqrt((x[2]-xmip)*(x[2]-xmip)+(y[2]-ymip)*(y[2]-ymip));
567// compact the infos...
568 Int_t compact = (Int_t)(thetaC*TMath::RadToDeg())+1000*(Int_t)(thTrk*TMath::RadToDeg());
9785d5fb 569 Int_t bin = fgDatabase->FindBin(distA,distB);
570 if(fgDatabase->GetBinContent(bin)==0) fgDatabase->Fill(distA,distB,compact);
5a3482a0 571 }
572 }
573
9785d5fb 574 FillZeroChan();
575// fgDatabase = deconv;
5a3482a0 576
577 timer.Stop();
578 Double_t nSecs = timer.CpuTime();
579 AliInfo(Form("database HTA successfully open in %3.1f sec.(CPU). Reconstruction is started.",nSecs));
9785d5fb 580
581// pout->Write();
582// pout->Close();
583
9a573d52 584}//InitDatabase()
5a3482a0 585//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
9785d5fb 586void AliHMPIDReconHTA::FillZeroChan()const
5a3482a0 587{
606697a8 588 //If fills eventually channel without entries
589 //inthe histo "database" jyst interpolating the neighboring cells
590 // Arguments: histogram pointer of the database
591 // Returns: none
592 //
9785d5fb 593 Int_t nbinx = fgDatabase->GetNbinsX();
594 Int_t nbiny = fgDatabase->GetNbinsY();
5a3482a0 595 for(Int_t i = 0;i<nbinx;i++) {
596 for(Int_t j = 0;j<nbiny;j++) {
9785d5fb 597 if(fgDatabase->GetBinContent(i,j) == 0) {
5a3482a0 598 Int_t nXmin = i-1; Int_t nXmax=i+1;
599 Int_t nYmin = j-1; Int_t nYmax=j+1;
600 Int_t nc = 0;
601 Double_t meanC =0;
602 Double_t meanTrk =0;
603 for(Int_t ix=nXmin;ix<=nXmax;ix++) {
604 if(ix<0||ix>nbinx) continue;
605 for(Int_t iy=nYmin;iy<=nYmax;iy++) {
606 if(iy<0||iy>nbiny) continue;
9785d5fb 607 meanC += (Int_t)fgDatabase->GetBinContent(ix,iy)%1000;
608 meanTrk+= (Int_t)fgDatabase->GetBinContent(ix,iy)/1000;
5a3482a0 609 nc++;
610 }
611 meanC/=nc; meanTrk/=nc;
612 Int_t compact = (Int_t)meanC+1000*(Int_t)meanTrk;
9785d5fb 613 if(compact>0)fgDatabase->SetCellContent(i,j,compact);
5a3482a0 614 }
615 }
616 }
617 }
9a573d52 618}//FillZeroChan()
5a3482a0 619//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
9a573d52 620
621// stima gli angoli con il metodo dei minimi quadrati che sfrutta le distanze...
622
623Double_t AliHMPIDReconHTA::FindSimmPhi()
624{
625 //It finds the phi of the ring
626 //by using the min. dist. algorithm
627 // Arguments: none
628 // Returns: phi
629 //
630
631 Float_t xrotsumm =0; Float_t yrotsumm =0; Float_t xx =0;
632 Float_t yy =0; Float_t xy =0;
633
634 Int_t np=0;
635
636 for(Int_t i=0;i<fNClu;i++) {
637 if(!fClCk[i]) continue;
638 np++;
639 xrotsumm+=fXClu[i]; // summ xi
640 yrotsumm+=fYClu[i]; // summ yi
641 xx+=fXClu[i]*fXClu[i]; // summ xixi
642 yy+=fYClu[i]*fYClu[i]; // summ yiyi
643 xy+=fXClu[i]*fYClu[i]; // summ yixi
644 }
645
646 //_____calc. met min quadr using effective distance _________________________________________________
647
648 Double_t coeff[3];
649 coeff[0]= xy-xrotsumm*yrotsumm/np;
650 coeff[1]= yrotsumm*yrotsumm/np - xrotsumm*xrotsumm/np - yy + xx;
651 coeff[2]= xrotsumm*yrotsumm/np- xy;
652
653 //____________________________________________________________________________________________________
654
655 Double_t m1=0, m2=0;
656 Double_t n1=0, n2=0;
657
658 r2(coeff,m1,m2);
659
660 n1=(yrotsumm-m1*xrotsumm)/np;
661 n2=(yrotsumm-m2*xrotsumm)/np;
662
663 // le due soluzioni.................
664
665 Double_t d1 =(1/(m1*m1+1))*(yy+m1*m1*xx+np*n1*n1-2*m1*xy-2*n1*yrotsumm+2*m1*n1*xrotsumm);
666 Double_t d2 =(1/(m2*m2+1))*(yy+m2*m2*xx+np*n2*n2-2*m2*xy-2*n2*yrotsumm+2*m2*n2*xrotsumm);
667 Double_t mMin;
668 if(d1 > d2) mMin = m2; else mMin = m1;
669
670 Double_t PhiTrk= TMath::ATan(mMin);
671
672 // positive angles...
673 if(PhiTrk<0) PhiTrk+=180*TMath::DegToRad();
674
675 return PhiTrk;
676}
677
678//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
679Int_t AliHMPIDReconHTA::r2(Double_t *coef, Double_t &x1, Double_t &x2)
680{
681 //2nd deg. equation
682 //solution
683 // Arguments: coef 2 1 0: ax^2+bx+c=0
684 // Returns: n. of solutions
685 // x1= 1st sol
686 // x2= 2nd sol
687 Double_t a,b,c;
688 a = coef[2];
689 b = coef[1];
690 c = coef[0];
691 Double_t delta = b*b-4*a*c;
692 if(delta<0) {return 0;}
693 if(delta==0) {
694 x1=x2=-b/(2*a);
695 return 1;
696 }
697 // delta>0
698 x1 = (-b+TMath::Sqrt(delta))/(2*a);
699 x2 = (-b-TMath::Sqrt(delta))/(2*a);
700 return 2;
701}//r2()