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