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