]> 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
Comment...uncomment...
[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():
40 TTask("RichRec","RichPat"),
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//..
94 if(fXClu) delete fXClu;
95 if(fYClu) delete fYClu;
2e81abb0 96 if(fPhiPhot) delete fPhiPhot;
97 if(fThetaPhot) delete fThetaPhot;
5a3482a0 98 if(fClCk) delete fClCk;
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
5a3482a0 225 Double_t *phiphot = new Double_t[fNClu];
226 Double_t *dist = new Double_t[fNClu];
227 Int_t *indphi = new Int_t[fNClu];
228
229 Bool_t status;
230
9785d5fb 231// Sort in phi angle...
5a3482a0 232 for(Int_t i=0;i<fNClu;i++) {
9785d5fb 233 if(!fClCk[i]) {
2e81abb0 234 AliDebug(1,Form(" n.%3i xMIP %8.3f yMIP %8.3f check %i",i,fMipX,fMipY,fClCk[i]));
9785d5fb 235 phiphot[i] = 999.;
2e81abb0 236 dist[i] = 999.;
9785d5fb 237 continue;
238 }
5a3482a0 239 phiphot[i] = (TMath::ATan2(fMipY-fYClu[i],fMipX-fXClu[i])+TMath::Pi())*TMath::RadToDeg();
240 dist[i]=TMath::Sqrt((fMipX-fXClu[i])*(fMipX-fXClu[i])+(fMipY-fYClu[i])*(fMipY-fYClu[i]));
2e81abb0 241 AliDebug(1,Form(" n.%3i phiphot %8.3f dist %8.3f check %i",i,phiphot[i],dist[i],fClCk[i]));
5a3482a0 242 }
243
244 TMath::Sort(fNClu,phiphot,indphi,kFALSE);
245
246// Purify with a truncated mean;
247 Int_t np=0;
248 Double_t dMean = 0;
249 Double_t dMean2 = 0;
250 for(Int_t i=0;i<fNClu;i++) {
251 if(!fClCk[indphi[i]]) continue; // Check if a good photon candidate or not
252 dMean +=dist[indphi[i]];
253 dMean2+=dist[indphi[i]]*dist[indphi[i]];
254 np++;
255 }
256
257 dMean /=(Double_t)np;
258 dMean2 /=(Double_t)np;
259 Double_t rms = TMath::Sqrt(dMean2 - dMean*dMean);
260
261 for(Int_t i=0;i<fNClu;i++) {
262 if(!fClCk[indphi[i]]) continue; // Check if a good photon candidate or not
2e81abb0 263 if(TMath::Abs(dMean-dist[indphi[i]]) > 1.5*rms) {
5a3482a0 264 fClCk[indphi[i]] = kFALSE;
265 continue;
266 }
267 }
268
2e81abb0 269 AliDebug(1,"Purification of photons...");
5a3482a0 270//
271// purify vectors for good photon candidates
272//
273 Int_t npeff=0;
274 Double_t *phiphotP = new Double_t[fNClu+1];
275 Double_t *distP = new Double_t[fNClu+1];
276 for(Int_t i=0;i<fNClu;i++) {
2e81abb0 277 AliDebug(1,Form(" n. %3i phiphot %8.3f dist %8.3f check %i",i,phiphot[indphi[i]],dist[indphi[i]],fClCk[indphi[i]]));
5a3482a0 278 if(!fClCk[indphi[i]]) continue; // Check if a good photon candidate or not
279 phiphotP[npeff] = phiphot[indphi[i]];
280 distP[npeff] = dist[indphi[i]];
281 npeff++;
282 }
283
284 delete [] phiphot;
285 delete [] dist;
286 delete [] indphi;
287
288 if(npeff<3) {
2e81abb0 289 AliDebug(1,Form("FindShape failed: no enough photons = %i...",npeff));
5a3482a0 290 delete [] phiphotP;
291 delete [] distP;
292 return kFALSE;
293 }
294
5a3482a0 295 Double_t xA,xB;
9785d5fb 296 status = kFALSE;
ee2f3539 297
2e81abb0 298 if(!ShapeModel(npeff,phiphotP,distP,xA,xB,phiTrkRec)) {AliDebug(1,Form("ShapeModel failed ")); return kFALSE;}
299
300// if(xA > 50 || xB > 15) {AliDebug(1,Form("xA and xB failed out of range")); return kFALSE;}
5a3482a0 301
170a4fc5 302 Int_t binxDB,binyDB;
2e81abb0 303 Int_t compactDB=-1;
304
305 if(xA > xB) { //geometrically not possible, try to recover on a mean values...
306
307 FindBinDB(xA,xA,binxDB,binyDB);
308 if(binxDB<0 || binyDB<0) {AliDebug(1,Form("bin < 0 ! failed ")); return kFALSE;}
309 Int_t compactDB1 = CompactDB(binxDB,binyDB);
310 FindBinDB(xB,xB,binxDB,binyDB);
311 if(binxDB<0 || binyDB<0) {AliDebug(1,Form("bin < 0 ! failed ")); return kFALSE;}
312 Int_t compactDB2 = CompactDB(binxDB,binyDB);
313 Double_t thetaCRec1 = (Double_t)(compactDB1%1000);
314 Double_t thetaCRec2 = (Double_t)(compactDB2%1000);
315 Double_t thTrkRec1 = (Double_t)(compactDB1/1000);
316 Double_t thTrkRec2 = (Double_t)(compactDB2/1000);
317 thetaCRec = 0.5*(thetaCRec1+thetaCRec2);
318 thTrkRec = 0.5*( thTrkRec1+ thTrkRec2);
319
320 } else {
321
322 FindBinDB(xA,xB,binxDB,binyDB);
323 if(binxDB<0 || binyDB<0) {AliDebug(1,Form("bin < 0 ! failed ")); return kFALSE;}
5a3482a0 324
2e81abb0 325 compactDB = CompactDB(binxDB,binyDB);
326
327 if(compactDB<0) {AliDebug(1,Form("compact< 0! failed ")); return kFALSE;}
328 //
329 //
330 thetaCRec = (Double_t)(compactDB%1000);
331 thTrkRec = (Double_t)(compactDB/1000);
332
333 }
334
335 AliDebug(1,Form(" CompactDB %i thTrkRec %8.3f thetaCRec %8.3f ",compactDB,thTrkRec,thetaCRec));
336
337 phiTrkRec *= TMath::DegToRad();
ee2f3539 338 thTrkRec *= TMath::DegToRad();
339 thetaCRec *= TMath::DegToRad();
9785d5fb 340
ee2f3539 341 status = kTRUE;
5a3482a0 342
343 delete [] phiphotP;
344 delete [] distP;
345
346 return status;
347}
348//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
349Bool_t AliHMPIDReconHTA::ShapeModel(Int_t np,Double_t *phiphot,Double_t *dist,Double_t &xA,Double_t &xB,Double_t &phiStart)
350{
aa00f952 351// Find a Spline curve to define dist. vs. phi angle
352// in order to estimate the phi of the track
353// Arguments: np - # points corresponding to # photon candidates
354// dist - distance of each photon from MIP
355// phiphot - phi of the photon in the DRS
356// Returns: xA - min. distance from MIP
357// xB - dist. from mip perpedicular to the major axis
358// phiStart- estimate of the track phi
359
5a3482a0 360 TGraph *phigr = new TGraph(np,phiphot,dist);
2e81abb0 361 phiStart = FindSimmPhi();
362
363 Double_t phiStart1 = phiStart;
364 if(phiStart1 > 360) phiStart1 -= 360;
365 Double_t phiStart2 = phiStart+90;
366 if(phiStart2 > 360) phiStart2 -= 360;
367 xA = phigr->Eval(phiStart1);
368 xB = phigr->Eval(phiStart2);
369 //---
370 AliDebug(1,Form("phiStart %f phiStart1 %f phiStart2 %f ",phiStart,phiStart1,phiStart2));
371 AliDebug(1,Form("xA %f xB %f",xA,xB));
5a3482a0 372
2e81abb0 373 phiStart += 180;
374 if(phiStart>360) phiStart-=360;
5a3482a0 375
5a3482a0 376 return kTRUE;
377}
378//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
aa00f952 379Double_t AliHMPIDReconHTA::VertParab(Double_t x1,Double_t y1,Double_t x2, Double_t y2, Double_t x3, Double_t y3)const
5a3482a0 380{
aa00f952 381// It uses parabola from 3 points to evaluate the x-coord of the parab
382// Arguments: xi,yi - points
383// Returns: x-coord of the vertex
384
5a3482a0 385 Double_t a = ((x1-x3)*(y1-y2)-(x1-x2)*(y1-y3))/((x1*x1-x2*x2)*(x1-x3)-(x1*x1-x3*x3)*(x1-x2));
386 Double_t b = (y1-y2 - a*(x1*x1-x2*x2))/(x1-x2);
387// Double_t c = y1 - a*x1*x1-b*x1;
388 return -0.5*b/a;
389}
390//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2e81abb0 391Bool_t AliHMPIDReconHTA::FitRing(Double_t thTrkRec,Double_t phiTrkRec)
392{
393 Double_t th = thTrkRec;
394 Double_t ph = phiTrkRec;
395
396 FitFree(th,ph);
397 while(FitStatus()) {
398 th = fThTrkFit;
399 ph = fPhTrkFit;
400 FitFree(th,ph);
401 }
402 return kTRUE;
403}
404//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
5a3482a0 405Bool_t AliHMPIDReconHTA::FitFree(Double_t thTrkRec,Double_t phiTrkRec)
406{
407// Fit performed by minimizing RMS/sqrt(n) of the
408// photons reconstructed. First phi is fixed and theta
409// is fouond, then (th,ph) of the track
410// as free parameters
411// Arguments: PhiRec phi of the track
412// Returns: none
413
5a3482a0 414 TMinuit *pMinuit = new TMinuit(2);
415 pMinuit->mncler();
416 gMinuit->SetObjectFit((TObject*)this); gMinuit->SetFCN(AliHMPIDReconHTA::FunMinPhot); //set fit function
417 Double_t aArg=-1,parStep,parLow,parHigh; Int_t iErrFlg; //tmp vars for TMinuit
418 Double_t d1,d2,d3;
419 TString sName;
420 Double_t th,ph;
421
422 pMinuit->mnexcm("SET PRI",&aArg,1,iErrFlg); //suspend all printout from TMinuit
423 pMinuit->mnexcm("SET NOW",&aArg,0,iErrFlg);
424
425 if(thTrkRec==0) thTrkRec = 3.*TMath::DegToRad(); // not to start from the edge...
426
2e81abb0 427 AliDebug(1,Form(" Minimization STARTED with phiTRK %6.2f thetaTRK %6.2f",phiTrkRec,thTrkRec));
428
5a3482a0 429 pMinuit->mnparm(0," thTrk ",thTrkRec ,parStep=0.01,parLow=0,parHigh=TMath::PiOver4(),iErrFlg);
430 pMinuit->mnparm(1," phiTrk ",phiTrkRec,parStep=0.01,parLow=0,parHigh=TMath::TwoPi(),iErrFlg);
431
432 pMinuit->FixParameter(1);
433 pMinuit->mnexcm("SIMPLEX" ,&aArg,0,iErrFlg);
434 pMinuit->mnexcm("MIGRAD" ,&aArg,0,iErrFlg);
435 pMinuit->Release(1);
436 pMinuit->mnexcm("MIGRAD" ,&aArg,0,iErrFlg);
437
438 pMinuit->mnpout(0,sName,th,d1,d2,d3,iErrFlg);
439 pMinuit->mnpout(1,sName,ph,d1,d2,d3,iErrFlg);
440
441 Double_t f,par[2];
442 Double_t *grad=0x0;
443 par[0] = th;par[1] = ph;
444 pMinuit->Eval(2,grad,f,par,3);
445
5a3482a0 446 SetTrkFit(th,ph);
447 return kTRUE;
448}
449//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
450void AliHMPIDReconHTA::FunMinPhot(Int_t &/* */,Double_t* /* */,Double_t &f,Double_t *par,Int_t iflag)
451{
452// Minimization function to find best track and thetaC parameters
453// Arguments: f = function value to minimize
454// par = list of parameter to find
455// iflag = flag status. See Minuit instructions
456// Returns: none
457//
458// Note: it is necessary to call an instance of AlihMPIDParam. Not possible to use fParam
459// because of the static instantiation of the function in Minuit
460
461 AliHMPIDParam *pParam=AliHMPIDParam::Instance();
462 AliHMPIDReconHTA *pRecHTA=(AliHMPIDReconHTA*)gMinuit->GetObjectFit();
463 AliHMPIDRecon pRec;
464 Double_t sizeCh = 0.5*pParam->RadThick()+pParam->WinThick()+pParam->GapThick();
465 Double_t thTrk = par[0];
466 Double_t phTrk = par[1];
467 Double_t xrad = pRecHTA->MipX() - sizeCh*TMath::Tan(thTrk)*TMath::Cos(phTrk);
468 Double_t yrad = pRecHTA->MipY() - sizeCh*TMath::Tan(thTrk)*TMath::Sin(phTrk);
469 pRecHTA->SetRadXY(xrad,yrad);
470 pRec.SetTrack(xrad,yrad,thTrk,phTrk);
471
472 Double_t meanCkov =0;
473 Double_t meanCkov2=0;
474 Double_t thetaCer,phiCer;
475 Int_t nClAcc = 0;
476 Int_t nClTot=pRecHTA->NClu();
2e81abb0 477
5a3482a0 478 for(Int_t i=0;i<nClTot;i++) {
2e81abb0 479
480 if(pRecHTA->IdxMip() == i) {
481 pRecHTA->SetPhotAngles(i,999.,999.);
482 continue;
483 }
484
5a3482a0 485 if(!(pRecHTA->ClCk(i))) continue;
2e81abb0 486
487 Bool_t status = pRec.FindPhotCkov(pRecHTA->XClu(i),pRecHTA->YClu(i),thetaCer,phiCer);
488 if(!status) {
489 pRecHTA->SetPhotAngles(i,999.,999.);
490 continue;
491 }
492 pRecHTA->SetPhotAngles(i,thetaCer,phiCer);
5a3482a0 493 meanCkov += thetaCer;
494 meanCkov2 += thetaCer*thetaCer;
495 nClAcc++;
496 }
2e81abb0 497
5a3482a0 498 if(nClAcc==0) {f=999;return;}
499 meanCkov /=(Double_t)nClAcc;
500 meanCkov2 /=(Double_t)nClAcc;
501 Double_t rms = TMath::Sqrt(TMath::Abs(meanCkov2 - meanCkov*meanCkov));
502 f = rms/TMath::Sqrt((Double_t)nClAcc);
503
504 if(iflag==3) {
2e81abb0 505 Int_t nClAccStep1 = nClAcc;
5a3482a0 506 nClAcc = 0;
507 Double_t meanCkov1=0;
876a9d77 508 Double_t meanCkov3=0;
5a3482a0 509 for(Int_t i=0;i<nClTot;i++) {
2e81abb0 510
511 if(!(pRecHTA->ClCk(i))) continue;
512 thetaCer = pRecHTA->PhotTheta(i);
5a3482a0 513 if(TMath::Abs(thetaCer-meanCkov)<2*rms) {
514 meanCkov1 += thetaCer;
876a9d77 515 meanCkov3 += thetaCer*thetaCer;
5a3482a0 516 nClAcc++;
2e81abb0 517 } else pRecHTA->SetClCk(i,kFALSE);
518 }
519
520 if(nClAcc<3) {
521 pRecHTA->SetFitStatus(kFALSE);
522 pRecHTA->SetCkovFit(meanCkov);
523 pRecHTA->SetCkovSig2(rms*rms);
524 pRecHTA->SetNCluFit(nClAccStep1);
525 return;
5a3482a0 526 }
2e81abb0 527
5a3482a0 528 meanCkov1/=nClAcc;
876a9d77 529 Double_t rms2 = (meanCkov3 - meanCkov*meanCkov*nClAcc)/nClAcc;
2e81abb0 530 AliLog::Instance();
531 if(nClAcc!=nClAccStep1) pRecHTA->SetFitStatus(kTRUE); else pRecHTA->SetFitStatus(kFALSE);
532
533 pRecHTA->SetCkovFit(meanCkov1);
534 pRecHTA->SetCkovSig2(rms2);
535 pRecHTA->SetNCluFit(nClAcc);
5a3482a0 536 }
537}//FunMinPhot()
538//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
539void AliHMPIDReconHTA::InitDatabase()
540{
aa00f952 541// Construction a database of ring shapes on fly
542// Arguments: none
543// Returns : none
170a4fc5 544// N.B. fgDB is the distance with x-min from MIP
545// y-dist from the ring of the MIP perpendicular to major axis
aa00f952 546// The content is the packed info of track theta and thetaC in degrees
547// thetaC+1000*thTrk
548//
9785d5fb 549// TFile *pout = new TFile("./database.root","recreate");
ee2f3539 550
170a4fc5 551 static Bool_t isDone = kFALSE;
552
5a3482a0 553 TStopwatch timer;
170a4fc5 554
555 if(isDone) {
5565f017 556 return;
557 }
170a4fc5 558
559 if(!isDone) {
560 timer.Start();
561 }
562
5a3482a0 563 AliInfo(Form("database HTA is being built.Please, wait..."));
564//
c61a7285 565 Double_t x[3]={0,0,0},y[3];
5a3482a0 566
567 AliHMPIDRecon rec;
568
569 if(!fParam) fParam=AliHMPIDParam::Instance();
570 Double_t thetaMax = TMath::ACos(1./fParam->MeanIdxRad());
571 Double_t thTrkMax = 1./TMath::ASin(fParam->MeanIdxRad());
572
573 Int_t nstepx = 1000;
574 Int_t nstepy = 1000;
575
5a3482a0 576 //
577 Double_t xrad = 0;
578 Double_t yrad = 0;
579 Double_t phTrk = 0;
580
581 for(Int_t i=0;i<nstepx;i++) { //loop on thetaC
582 for(Int_t j=0;j<nstepy;j++) { //loop on theta particle
583 Double_t thetaC = thetaMax/nstepx*((Double_t)i+0.5);
584 Double_t thTrk = thTrkMax/nstepy*((Double_t)j+0.5);
585 //
586 //mip position
587 //
fcaff63d 588 Double_t sizeCh = 0.5*fParam->RadThick()+fParam->WinThick()+fParam->GapThick();
5a3482a0 589 Double_t xmip = xrad + sizeCh*TMath::Tan(thTrk)*TMath::Cos(phTrk);
590 Double_t ymip = yrad + sizeCh*TMath::Tan(thTrk)*TMath::Sin(phTrk);
591
592 Double_t dist1,dist2;
593 //
594 //first point at phi=0
595 //
596 rec.SetTrack(xrad,yrad,thTrk,phTrk);
597 TVector2 pos;
598 pos=rec.TracePhot(thetaC,0);
599
600 if(pos.X()==-999) {
ee2f3539 601 dist1 = 0; //open ring...only the distance btw mip and point at 180 will be considered
5a3482a0 602 } else {
603 x[0] = pos.X(); y[0] = pos.Y();
604 dist1 = TMath::Sqrt((x[0]-xmip)*(x[0]-xmip)+(y[0]-ymip)*(y[0]-ymip));
605 }
606 //
607 //second point at phi=180
608 //
609 rec.SetTrack(xrad,yrad,thTrk,phTrk);
610 pos=rec.TracePhot(thetaC,TMath::Pi());
611
2e81abb0 612 if(pos.X()==-999) {AliDebug(1,Form("it should not happens!Bye"));return;}
5a3482a0 613 x[1] = pos.X(); y[1] = pos.Y();
614 if((x[1]-xmip)*(x[0]-xmip)>0) continue; // to avoid circles out mips (for very low ThetaC)
615 dist2 = TMath::Sqrt((x[1]-xmip)*(x[1]-xmip)+(y[1]-ymip)*(y[1]-ymip));
616
617// Double_t distA = dist1+dist2;
618 Double_t distA = dist2; // only the minimum: problem of acceptance
619 //
620 //second point at phi=90
621 //
622 rec.SetTrack(xrad,yrad,thTrk,phTrk);
623 pos=rec.TracePhot(thetaC,TMath::PiOver2());
624
625 if(pos.X()==-999) continue;
626 x[2] = pos.X(); y[2] = pos.Y();
627 Double_t distB = TMath::Sqrt((x[2]-xmip)*(x[2]-xmip)+(y[2]-ymip)*(y[2]-ymip));
628// compact the infos...
629 Int_t compact = (Int_t)(thetaC*TMath::RadToDeg())+1000*(Int_t)(thTrk*TMath::RadToDeg());
ee2f3539 630 Int_t binxDB,binyDB;
631 FindBinDB(distA,distB,binxDB,binyDB);
170a4fc5 632 if(fgDB[binxDB][binyDB]==0) fgDB[binxDB][binyDB] = compact;
5a3482a0 633 }
634 }
635
9785d5fb 636 FillZeroChan();
5a3482a0 637
170a4fc5 638 if(!isDone) {
639 timer.Stop();
640 Double_t nSecs = timer.CpuTime();
641 AliInfo(Form("database HTA successfully open in %3.1f sec.(CPU). Reconstruction is started.",nSecs));
642 isDone = kTRUE;
643 }
9785d5fb 644
645// pout->Write();
646// pout->Close();
647
9a573d52 648}//InitDatabase()
5a3482a0 649//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
9785d5fb 650void AliHMPIDReconHTA::FillZeroChan()const
5a3482a0 651{
606697a8 652 //If fills eventually channel without entries
653 //inthe histo "database" jyst interpolating the neighboring cells
654 // Arguments: histogram pointer of the database
655 // Returns: none
656 //
170a4fc5 657 const Int_t nxDB = 500;
658 const Int_t nyDB = 150;
659
660 for(Int_t i = 0;i<nxDB;i++) {
661 for(Int_t j = 0;j<nyDB;j++) {
662 if(fgDB[i][j] == 0) {
663 fgDB[i][j] = -1;
5a3482a0 664 Int_t nXmin = i-1; Int_t nXmax=i+1;
665 Int_t nYmin = j-1; Int_t nYmax=j+1;
666 Int_t nc = 0;
667 Double_t meanC =0;
668 Double_t meanTrk =0;
669 for(Int_t ix=nXmin;ix<=nXmax;ix++) {
170a4fc5 670 if(ix<0||ix>=nxDB) continue;
5a3482a0 671 for(Int_t iy=nYmin;iy<=nYmax;iy++) {
170a4fc5 672 if(iy<0||iy>=nyDB) continue;
673 meanC += (Int_t)(fgDB[ix][iy]%1000);
674 meanTrk+= (Int_t)(fgDB[ix][iy]/1000);
5a3482a0 675 nc++;
676 }
677 meanC/=nc; meanTrk/=nc;
678 Int_t compact = (Int_t)meanC+1000*(Int_t)meanTrk;
170a4fc5 679 if(compact>0)fgDB[i][j] = compact;
5a3482a0 680 }
681 }
682 }
683 }
9a573d52 684}//FillZeroChan()
5a3482a0 685//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
ee2f3539 686Int_t AliHMPIDReconHTA::r2(Double_t *coef, Double_t &x1, Double_t &x2)
687{
688 //2nd deg. equation
689 //solution
690 // Arguments: coef 2 1 0: ax^2+bx+c=0
691 // Returns: n. of solutions
692 // x1= 1st sol
693 // x2= 2nd sol
694 Double_t a,b,c;
695 a = coef[2];
696 b = coef[1];
697 c = coef[0];
698 Double_t delta = b*b-4*a*c;
699 if(delta<0) {return 0;}
700 if(delta==0) {
701 x1=x2=-b/(2*a);
702 return 1;
703 }
704 // delta>0
705 x1 = (-b+TMath::Sqrt(delta))/(2*a);
706 x2 = (-b-TMath::Sqrt(delta))/(2*a);
707 return 2;
708}//r2()
709//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
9a573d52 710
711Double_t AliHMPIDReconHTA::FindSimmPhi()
ee2f3539 712{
713//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2e81abb0 714// Reconstruction of phiTRK angle with two methods (in switching)
715//
716// - least square method (for closed rings)
717// - by minimum distance mip-photon (for open rings)
ee2f3539 718
719 Float_t coeff1ord=0; Float_t coeff2ord=0; Float_t coeff0ord=0;
9a573d52 720 Float_t xrotsumm =0; Float_t yrotsumm =0; Float_t xx =0;
2e81abb0 721 Float_t yy =0; Float_t xy =0;
ee2f3539 722 Double_t xmin=0;
723 Double_t ymin=0;
9a573d52 724
725 Int_t np=0;
ee2f3539 726
727 Double_t distMin = 999.;
728
9a573d52 729 for(Int_t i=0;i<fNClu;i++) {
730 if(!fClCk[i]) continue;
731 np++;
732 xrotsumm+=fXClu[i]; // summ xi
733 yrotsumm+=fYClu[i]; // summ yi
734 xx+=fXClu[i]*fXClu[i]; // summ xixi
735 yy+=fYClu[i]*fYClu[i]; // summ yiyi
736 xy+=fXClu[i]*fYClu[i]; // summ yixi
ee2f3539 737 Double_t dist2= (fXClu[i]-fMipX)*(fXClu[i]-fMipX)+(fYClu[i]-fMipY)*(fYClu[i]-fMipY);
738 if(dist2<distMin) {
739 distMin = dist2;
740 xmin = fXClu[i];
741 ymin = fYClu[i];
742 }
9a573d52 743 }
ee2f3539 744
2e81abb0 745 Double_t AngM = TMath::ATan2(ymin-fMipY,(xmin-fMipX))*TMath::RadToDeg();
746 if (AngM<0) AngM+=360;
747
748 AliDebug(1,Form(" Simple angle prediction with MIN phi = %f",AngM));
9a573d52 749
ee2f3539 750 //_____calc. met min quadr using effective distance _________________________________________________
9a573d52 751
2e81abb0 752 coeff2ord = xy-xrotsumm*yrotsumm/np;
753 coeff1ord = yrotsumm*yrotsumm/np - xrotsumm*xrotsumm/np - yy + xx;
754 coeff0ord = -coeff2ord;
755
756 AliDebug(1,Form(" a = %f b = %f c = %f",coeff2ord,coeff1ord,coeff0ord));
757
ee2f3539 758 Double_t m1=0, m2=0; Double_t n1=0, n2=0;
759 // c // b // a
760 Double_t coeff[3]={coeff0ord,coeff1ord,coeff2ord};
9a573d52 761
762 r2(coeff,m1,m2);
763
764 n1=(yrotsumm-m1*xrotsumm)/np;
765 n2=(yrotsumm-m2*xrotsumm)/np;
ee2f3539 766 // 2 solutions.................
767
ee2f3539 768 // negative angles solved...
769
9a573d52 770 Double_t d1 =(1/(m1*m1+1))*(yy+m1*m1*xx+np*n1*n1-2*m1*xy-2*n1*yrotsumm+2*m1*n1*xrotsumm);
771 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 772
773 AliDebug(1,Form(" predicted distance d1 %f for angle %6.2f d2 %f for angle %6.2f",d1,TMath::ATan(m1)*TMath::RadToDeg(),
774 d2,TMath::ATan(m2)*TMath::RadToDeg()));
9a573d52 775 Double_t mMin;
776 if(d1 > d2) mMin = m2; else mMin = m1;
777
2e81abb0 778 Double_t phiTrk=0;
ee2f3539 779 //
2e81abb0 780 if(ymin < fMipY && xmin > fMipX) {phiTrk = TMath::ATan(mMin)*TMath::RadToDeg()+180;}
781 if(ymin > fMipY && xmin < fMipX) {phiTrk = TMath::ATan(mMin)*TMath::RadToDeg()+360;}
782 if(ymin > fMipY && xmin > fMipX) {phiTrk = TMath::ATan(mMin)*TMath::RadToDeg()+180;}
783 if(ymin < fMipY && xmin < fMipX) {phiTrk = TMath::ATan(mMin)*TMath::RadToDeg();}
784 if(ymin == fMipY && xmin > fMipX) {phiTrk = TMath::ATan(mMin)*TMath::RadToDeg()+180;}
785 if(ymin == fMipY && xmin < fMipX) {phiTrk = TMath::ATan(mMin)*TMath::RadToDeg();}
786 if(ymin < fMipY && xmin == fMipX) {phiTrk = 90;}
787 if(ymin > fMipY && xmin == fMipX) {phiTrk = 270;}
9a573d52 788
ee2f3539 789 // ------------------------- choose the best-----------------------
9a573d52 790
ee2f3539 791
2e81abb0 792 if( AngM-40 <= phiTrk && AngM+40 >= phiTrk) return phiTrk; else return AngM;
9a573d52 793}
9a573d52 794//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
ee2f3539 795void AliHMPIDReconHTA::FindBinDB(Double_t x,Double_t y,Int_t &binX,Int_t &binY)
9a573d52 796{
170a4fc5 797 const Int_t nxDB = 500;
798 const Int_t nyDB = 150;
ee2f3539 799 const Double_t xlowDB = 0;
800 const Double_t xhigDB = 50;
801 const Double_t ylowDB = 0;
2e81abb0 802 const Double_t yhigDB = 30;
ee2f3539 803
804 binX = -1;
805 binY = -1;
170a4fc5 806 if(x<xlowDB && x>xhigDB &&
807 y<ylowDB && y>yhigDB) return;
808 binX = Int_t((x-xlowDB)/(xhigDB-xlowDB)*nxDB);
809 binY = Int_t((y-ylowDB)/(yhigDB-ylowDB)*nyDB);
810 if(binX>=nxDB || binY>=nyDB) {
811 binX = -1;
812 binY = -1;
813 }
814
ee2f3539 815}
2e81abb0 816//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
817Bool_t AliHMPIDReconHTA::UniformDistrib()
818{
819 AliHMPIDParam *pParam=AliHMPIDParam::Instance();
820 AliHMPIDRecon pRec;
821
822 Double_t sizeCh = 0.5*pParam->RadThick()+pParam->WinThick()+pParam->GapThick();
823 Double_t xrad = MipX() - sizeCh*TMath::Tan(fThTrkFit)*TMath::Cos(fPhTrkFit);
824 Double_t yrad = MipY() - sizeCh*TMath::Tan(fThTrkFit)*TMath::Sin(fPhTrkFit);
825
826 pRec.SetTrack(xrad,yrad,fThTrkFit,fPhTrkFit);
827
828 Int_t npeff=0;
829 Int_t nPhotPerBin = 4;
830 for(Int_t i=0;i<fNClu;i++) {
831 if(!ClCk(i)) continue;
832 npeff++;
833 }
834
835 Int_t nPhiBins = npeff/nPhotPerBin+1;
836 if(nPhiBins<=1) nPhiBins = 2;
837
838 Double_t *iPhiBin;
839 iPhiBin = new Double_t[nPhiBins];
840
841 for(Int_t i=0;i<nPhiBins;i++) iPhiBin[i] =0;
842
843 for(Int_t i=0;i<fNClu;i++) {
844 if(!ClCk(i)) continue;
845 Double_t phiCer = PhotPhi(i);
846
847 Double_t PhiProva = phiCer*TMath::RadToDeg();
848 if(PhiProva<0) PhiProva+= 360;
849 Int_t index = (Int_t)((Float_t)nPhiBins*PhiProva/360.);
850 iPhiBin[index]++;
851 }
852
853 Double_t chi2 = 0;
854 for(Int_t i=0;i<nPhiBins;i++) {
855 Double_t theo = (Double_t)npeff/(Double_t)nPhiBins;
856 chi2+= (iPhiBin[i] - theo)*(iPhiBin[i] - theo)/theo;
857 }
858
859 delete iPhiBin;
860
861 Double_t prob = TMath::Prob(chi2, nPhiBins-1);
862 AliDebug(1,Form(" Probability for uniform distrib: %6f.3 %s",prob,(prob<0.05) ? "rejected" : "accepted"));
863 if(prob<0.05) return kFALSE;
864 return kTRUE;
865
866 }
867//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++