b3378cb69656eab0ab55f26e7436d3f42248c03d
[u/mrichter/AliRoot.git] / HMPID / AliHMPIDReconHTA.cxx
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   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),
50   fPhiPhot(0),
51   fThetaPhot(0),
52   fClCk(0),
53   fThTrkIn(-999),
54   fPhTrkIn(-999),
55   fThTrkFit(-999),
56   fPhTrkFit(-999),
57   fCkovFit(-999),
58   fNCluFit(0),
59   fCkovSig2(0),
60   fFitStatus(0),
61   fParam(AliHMPIDParam::Instance())
62 {
63 //..
64 //hidden algorithm
65 //..
66   fParam->SetRefIdx(fParam->MeanIdxRad()); // initialization of ref index to a default one
67   InitDatabase();
68 }
69 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
70 AliHMPIDReconHTA::~AliHMPIDReconHTA()
71 {
72   DeleteVars();
73 }
74 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
75 void AliHMPIDReconHTA::InitVars(Int_t n)
76 {
77 //..
78 //Init some variables
79 //..
80   fXClu = new Double_t[n];
81   fYClu = new Double_t[n];
82   fPhiPhot = new Double_t[n];
83   fThetaPhot = new Double_t[n];
84   fClCk = new Bool_t[n];
85   for(Int_t i=0;i<n;i++) fClCk[i] = kTRUE;
86 //
87 }
88 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
89 void AliHMPIDReconHTA::DeleteVars()const
90 {
91 //..
92 //Delete variables
93 //..
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;
99 }
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   AliHMPIDParam *pParam = AliHMPIDParam::Instance(); 
113
114   if(!CluPreFilter(pCluLst)) return kFALSE;
115
116   Int_t nCh=0;
117   Int_t sizeClu=0;
118   
119   fNClu = pCluLst->GetEntriesFast();
120     
121   for (Int_t iClu=0;iClu<fNClu;iClu++){                                                       //clusters loop
122     AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(iClu);                       //get pointer to current cluster    
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
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;
135       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   pParam->SetRefIdx(nmean);
140   
141   //
142   Float_t xra,yra,th,ph; pTrk->GetHMPIDtrk(xra,yra,th,ph);
143   AliDebug(1,Form(" simulated phiTRK %6.2f thetaTRK %6.2f",ph*TMath::RadToDeg(),th*TMath::RadToDeg()));
144   //
145   
146   if(!DoRecHiddenTrk()) {
147     pTrk->SetHMPIDsignal(pParam->kNoPhotAccept);
148     return kFALSE;
149   }                                                                           //Do track and ring reconstruction,if problems returns 1
150   AliDebug(1,Form("    fitted phi %6.2f ",fPhTrkFit*TMath::RadToDeg()));
151   
152   pTrk->SetHMPIDtrk(fRadX,fRadY,fThTrkFit,fPhTrkFit);                                        //store track intersection info
153   pTrk->SetHMPIDmip(fMipX,fMipY,(Int_t)fMipQ,NCluFit());                                     //store mip info + n. phots of the ring 
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
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()));
162   
163   return kTRUE;
164   
165 }//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   Double_t thTrkRec,phiTrkRec,thetaCRec;
178   
179   if(!FindShape(thTrkRec,phiTrkRec,thetaCRec)) {
180     AliDebug(1,Form(" FindShape failed...!"));
181     return kFALSE;
182   }
183   AliDebug(1,Form(" FindShape accepted...!"));
184
185   if(!FitRing(thTrkRec,phiTrkRec)) {
186     AliDebug(1,Form(" FitRing failed...!"));
187     return kFALSE;
188   }
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
197   return kTRUE;
198 }//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   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 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
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   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     
231 // Sort in phi angle...
232   for(Int_t i=0;i<fNClu;i++) {
233     if(!fClCk[i]) {
234       AliDebug(1,Form(" n.%3i  xMIP    %8.3f yMIP %8.3f check %i",i,fMipX,fMipY,fClCk[i]));
235       phiphot[i] = 999.;
236       dist[i]    = 999.;
237       continue;
238     }
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]));
241     AliDebug(1,Form(" n.%3i  phiphot %8.3f dist %8.3f check %i",i,phiphot[i],dist[i],fClCk[i]));
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
263     if(TMath::Abs(dMean-dist[indphi[i]]) > 1.5*rms) {
264       fClCk[indphi[i]] = kFALSE;
265       continue;
266     }
267   }
268
269   AliDebug(1,"Purification of photons...");
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++) {
277     AliDebug(1,Form(" n. %3i phiphot %8.3f dist %8.3f check %i",i,phiphot[indphi[i]],dist[indphi[i]],fClCk[indphi[i]]));
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) {
289     AliDebug(1,Form("FindShape failed: no enough photons = %i...",npeff));
290     delete [] phiphotP;
291     delete [] distP;
292     return kFALSE;
293   }
294
295   Double_t xA,xB;
296   status = kFALSE;
297   
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;}
301
302   Int_t binxDB,binyDB;
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;}
324
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();
338   thTrkRec  *= TMath::DegToRad(); 
339   thetaCRec *= TMath::DegToRad();
340
341   status = kTRUE;
342
343   delete [] phiphotP;
344   delete [] distP;
345   
346   return status;
347 }
348 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
349 Bool_t AliHMPIDReconHTA::ShapeModel(Int_t np,Double_t *phiphot,Double_t *dist,Double_t &xA,Double_t &xB,Double_t &phiStart)
350 {
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
360   TGraph *phigr = new TGraph(np,phiphot,dist);
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));
372
373   phiStart += 180;  
374   if(phiStart>360) phiStart-=360;
375   
376   return kTRUE;
377 }
378 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
379 Double_t AliHMPIDReconHTA::VertParab(Double_t x1,Double_t y1,Double_t x2, Double_t y2, Double_t x3, Double_t y3)const
380 {
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   
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 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
391 Bool_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 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
405 Bool_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   
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   
427   AliDebug(1,Form(" Minimization STARTED with phiTRK %6.2f thetaTRK %6.2f",phiTrkRec,thTrkRec));
428
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
446   SetTrkFit(th,ph);
447   return kTRUE;
448 }
449 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
450 void 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();
477
478   for(Int_t i=0;i<nClTot;i++) {
479     
480     if(pRecHTA->IdxMip() == i) {
481       pRecHTA->SetPhotAngles(i,999.,999.);
482       continue;
483     }
484     
485     if(!(pRecHTA->ClCk(i))) continue;
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);
493     meanCkov  += thetaCer;
494     meanCkov2 += thetaCer*thetaCer;
495     nClAcc++;
496   }
497   
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) {
505     Int_t nClAccStep1 = nClAcc;
506     nClAcc = 0;
507     Double_t meanCkov1=0;
508     Double_t meanCkov3=0;
509     for(Int_t i=0;i<nClTot;i++) {
510       
511       if(!(pRecHTA->ClCk(i))) continue;
512       thetaCer = pRecHTA->PhotTheta(i);
513       if(TMath::Abs(thetaCer-meanCkov)<2*rms) {
514         meanCkov1 += thetaCer;
515         meanCkov3 += thetaCer*thetaCer;
516         nClAcc++;
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;
526     }
527       
528     meanCkov1/=nClAcc;
529     Double_t rms2 = (meanCkov3 - meanCkov*meanCkov*nClAcc)/nClAcc;
530     
531     // get a logger instance
532     // what for??
533     AliLog::GetRootLogger();
534
535     if(nClAcc!=nClAccStep1) pRecHTA->SetFitStatus(kTRUE); else pRecHTA->SetFitStatus(kFALSE);
536     
537     pRecHTA->SetCkovFit(meanCkov1);
538     pRecHTA->SetCkovSig2(rms2);
539     pRecHTA->SetNCluFit(nClAcc);
540   }
541 }//FunMinPhot()
542 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
543 void AliHMPIDReconHTA::InitDatabase()
544 {
545 // Construction a database of ring shapes on fly
546 //   Arguments: none
547 //   Returns  : none
548 //  N.B. fgDB is the distance with x-min from MIP
549 //                                 y-dist from the ring of the MIP perpendicular to major axis
550 //        The content is the packed info of track theta and thetaC in degrees
551 //                        thetaC+1000*thTrk
552 //
553 //  TFile *pout = new TFile("./database.root","recreate");
554
555   static Bool_t isDone = kFALSE;
556   
557   TStopwatch timer;
558   
559   if(isDone) {
560    return;
561   }
562   
563   if(!isDone) {
564     timer.Start();
565   }
566  
567   AliInfo(Form("database HTA is being built.Please, wait..."));
568 //  
569   Double_t x[3]={0,0,0},y[3];
570
571   AliHMPIDRecon rec;
572
573   if(!fParam) fParam=AliHMPIDParam::Instance();
574   Double_t thetaMax = TMath::ACos(1./fParam->MeanIdxRad());
575   Double_t thTrkMax = 1./TMath::ASin(fParam->MeanIdxRad());    
576
577   Int_t nstepx = 1000;
578   Int_t nstepy = 1000;
579
580   //
581   Double_t xrad = 0;
582   Double_t yrad = 0;
583   Double_t phTrk = 0;
584
585   for(Int_t i=0;i<nstepx;i++) {     //loop on thetaC
586     for(Int_t j=0;j<nstepy;j++) {   //loop on theta particle
587       Double_t thetaC = thetaMax/nstepx*((Double_t)i+0.5);
588       Double_t thTrk  = thTrkMax/nstepy*((Double_t)j+0.5);
589       //
590       //mip position
591       //
592       Double_t sizeCh = 0.5*fParam->RadThick()+fParam->WinThick()+fParam->GapThick();
593       Double_t xmip = xrad + sizeCh*TMath::Tan(thTrk)*TMath::Cos(phTrk);
594       Double_t ymip = yrad + sizeCh*TMath::Tan(thTrk)*TMath::Sin(phTrk);
595
596       Double_t dist1,dist2;
597       //
598       //first point at phi=0
599       //
600       rec.SetTrack(xrad,yrad,thTrk,phTrk);
601       TVector2 pos;
602       pos=rec.TracePhot(thetaC,0);
603
604       if(pos.X()==-999) {
605         dist1 = 0;            //open ring...only the distance btw mip and point at 180 will be considered
606       } else {
607         x[0] = pos.X(); y[0] = pos.Y();
608         dist1   = TMath::Sqrt((x[0]-xmip)*(x[0]-xmip)+(y[0]-ymip)*(y[0]-ymip));
609       }
610       //
611       //second point at phi=180
612       //
613       rec.SetTrack(xrad,yrad,thTrk,phTrk);
614       pos=rec.TracePhot(thetaC,TMath::Pi());
615
616       if(pos.X()==-999) {AliDebug(1,Form("it should not happens!Bye"));return;}
617       x[1] = pos.X(); y[1] = pos.Y();
618       if((x[1]-xmip)*(x[0]-xmip)>0) continue; // to avoid circles out mips (for very low ThetaC)
619       dist2   = TMath::Sqrt((x[1]-xmip)*(x[1]-xmip)+(y[1]-ymip)*(y[1]-ymip));
620
621 //      Double_t distA = dist1+dist2;
622       Double_t distA = dist2;     // only the minimum: problem of acceptance
623       //
624       //second point at phi=90
625       //
626       rec.SetTrack(xrad,yrad,thTrk,phTrk);
627       pos=rec.TracePhot(thetaC,TMath::PiOver2());
628
629       if(pos.X()==-999) continue;
630       x[2] = pos.X(); y[2] = pos.Y();
631       Double_t distB   = TMath::Sqrt((x[2]-xmip)*(x[2]-xmip)+(y[2]-ymip)*(y[2]-ymip));
632 // compact the infos...      
633       Int_t compact = (Int_t)(thetaC*TMath::RadToDeg())+1000*(Int_t)(thTrk*TMath::RadToDeg());
634       Int_t binxDB,binyDB;
635       FindBinDB(distA,distB,binxDB,binyDB);
636       if(fgDB[binxDB][binyDB]==0) fgDB[binxDB][binyDB] = compact;
637     }
638   }
639
640   FillZeroChan();
641
642   if(!isDone) {
643     timer.Stop();
644     Double_t nSecs = timer.CpuTime();  
645     AliInfo(Form("database HTA successfully open in %3.1f sec.(CPU). Reconstruction is started.",nSecs));
646     isDone = kTRUE;
647   }
648   
649 //  pout->Write();
650 //  pout->Close();
651   
652 }//InitDatabase()
653 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
654 void AliHMPIDReconHTA::FillZeroChan()const
655 {
656   //If fills eventually channel without entries
657   //inthe histo "database" jyst interpolating the neighboring cells
658   // Arguments: histogram pointer of the database
659   //   Returns: none
660   //
661   const Int_t nxDB = 500;
662   const Int_t nyDB = 150;
663
664   for(Int_t i = 0;i<nxDB;i++) {
665     for(Int_t j = 0;j<nyDB;j++) {
666       if(fgDB[i][j] == 0) {
667         fgDB[i][j] = -1;
668         Int_t nXmin = i-1; Int_t nXmax=i+1;
669         Int_t nYmin = j-1; Int_t nYmax=j+1;
670         Int_t nc = 0;
671         Double_t meanC =0;
672         Double_t meanTrk =0;
673         for(Int_t ix=nXmin;ix<=nXmax;ix++) {
674           if(ix<0||ix>=nxDB) continue;
675           for(Int_t iy=nYmin;iy<=nYmax;iy++) {
676             if(iy<0||iy>=nyDB) continue;
677             meanC  +=  (Int_t)(fgDB[ix][iy]%1000);
678             meanTrk+=  (Int_t)(fgDB[ix][iy]/1000);
679             nc++;
680           }
681           meanC/=nc; meanTrk/=nc;
682           Int_t compact = (Int_t)meanC+1000*(Int_t)meanTrk;
683           if(compact>0)fgDB[i][j] = compact;
684         }
685       }
686     }
687   }
688 }//FillZeroChan()
689 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
690 Int_t AliHMPIDReconHTA::r2(Double_t *coef, Double_t &x1, Double_t &x2)
691 {
692   //2nd deg. equation
693   //solution
694   // Arguments: coef 2 1 0: ax^2+bx+c=0
695   //   Returns: n. of solutions
696   //            x1= 1st sol
697   //            x2= 2nd sol
698   Double_t a,b,c;
699   a = coef[2];
700   b = coef[1];
701   c = coef[0];
702   Double_t delta = b*b-4*a*c;
703   if(delta<0) {return 0;}
704   if(delta==0) {
705     x1=x2=-b/(2*a);
706     return 1;
707   }
708   // delta>0
709   x1 = (-b+TMath::Sqrt(delta))/(2*a);
710   x2 = (-b-TMath::Sqrt(delta))/(2*a);
711   return 2;
712 }//r2()
713 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
714
715 Double_t AliHMPIDReconHTA::FindSimmPhi() 
716 {     
717 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
718 // Reconstruction of phiTRK angle with two methods (in switching)
719 // 
720 // - least square method (for closed rings)
721 // - by minimum distance mip-photon (for open rings)
722   
723   Float_t coeff1ord=0;     Float_t coeff2ord=0;     Float_t coeff0ord=0;  
724   Float_t xrotsumm =0;     Float_t yrotsumm =0;     Float_t xx =0;           
725   Float_t yy =0;           Float_t xy =0;
726   Double_t xmin=0;
727   Double_t ymin=0;
728
729   Int_t np=0;    
730     
731   Double_t distMin = 999.;
732   
733   for(Int_t i=0;i<fNClu;i++) {
734     if(!fClCk[i]) continue;
735     np++;
736     xrotsumm+=fXClu[i];         // summ xi
737     yrotsumm+=fYClu[i];         // summ yi
738     xx+=fXClu[i]*fXClu[i];      // summ xixi     
739     yy+=fYClu[i]*fYClu[i];      // summ yiyi
740     xy+=fXClu[i]*fYClu[i];      // summ yixi     
741     Double_t dist2= (fXClu[i]-fMipX)*(fXClu[i]-fMipX)+(fYClu[i]-fMipY)*(fYClu[i]-fMipY);
742     if(dist2<distMin) {
743       distMin = dist2;
744       xmin = fXClu[i];
745       ymin = fYClu[i];
746     }
747   }
748
749   Double_t AngM = TMath::ATan2(ymin-fMipY,(xmin-fMipX))*TMath::RadToDeg();
750   if (AngM<0) AngM+=360;
751   
752   AliDebug(1,Form(" Simple angle prediction with MIN phi = %f",AngM));
753   
754   //_____calc. met min quadr using effective distance _________________________________________________
755   
756   coeff2ord = xy-xrotsumm*yrotsumm/np;    
757   coeff1ord = yrotsumm*yrotsumm/np - xrotsumm*xrotsumm/np - yy + xx;
758   coeff0ord = -coeff2ord;
759   
760   AliDebug(1,Form(" a = %f b = %f c = %f",coeff2ord,coeff1ord,coeff0ord));
761   
762   Double_t m1=0, m2=0;  Double_t n1=0, n2=0;
763                             // c           // b         // a
764   Double_t coeff[3]={coeff0ord,coeff1ord,coeff2ord};    
765   
766   r2(coeff,m1,m2);
767   
768   n1=(yrotsumm-m1*xrotsumm)/np;                         
769   n2=(yrotsumm-m2*xrotsumm)/np;
770   // 2 solutions.................
771   
772   // negative angles solved...
773   
774   Double_t d1 =(1/(m1*m1+1))*(yy+m1*m1*xx+np*n1*n1-2*m1*xy-2*n1*yrotsumm+2*m1*n1*xrotsumm);
775   Double_t d2 =(1/(m2*m2+1))*(yy+m2*m2*xx+np*n2*n2-2*m2*xy-2*n2*yrotsumm+2*m2*n2*xrotsumm);
776
777   AliDebug(1,Form(" predicted distance d1 %f for angle %6.2f d2 %f for angle %6.2f",d1,TMath::ATan(m1)*TMath::RadToDeg(),
778                                                                                     d2,TMath::ATan(m2)*TMath::RadToDeg()));
779   Double_t mMin;
780   if(d1 > d2) mMin = m2; else mMin = m1;
781   
782   Double_t phiTrk=0;
783   // 
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()+360;}
786   if(ymin >  fMipY && xmin >  fMipX)  {phiTrk  =  TMath::ATan(mMin)*TMath::RadToDeg()+180;}
787   if(ymin <  fMipY && xmin <  fMipX)  {phiTrk  =  TMath::ATan(mMin)*TMath::RadToDeg();}
788   if(ymin == fMipY && xmin >  fMipX)  {phiTrk  =  TMath::ATan(mMin)*TMath::RadToDeg()+180;}
789   if(ymin == fMipY && xmin <  fMipX)  {phiTrk  =  TMath::ATan(mMin)*TMath::RadToDeg();}
790   if(ymin <  fMipY && xmin == fMipX)  {phiTrk  =   90;}
791   if(ymin >  fMipY && xmin == fMipX)  {phiTrk  =  270;}
792   
793   //  ------------------------- choose the best-----------------------
794   
795   
796   if( AngM-40 <=  phiTrk  &&  AngM+40 >= phiTrk)   return phiTrk;  else return AngM;
797 }
798 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
799 void AliHMPIDReconHTA::FindBinDB(Double_t x,Double_t y,Int_t &binX,Int_t &binY)
800 {
801   const Int_t nxDB = 500;
802   const Int_t nyDB = 150;
803   const Double_t xlowDB =  0;   
804   const Double_t xhigDB = 50;   
805   const Double_t ylowDB =  0;
806   const Double_t yhigDB = 30;
807
808   binX = -1;
809   binY = -1;
810   if(x<xlowDB && x>xhigDB &&
811      y<ylowDB && y>yhigDB)    return;
812   binX = Int_t((x-xlowDB)/(xhigDB-xlowDB)*nxDB);
813   binY = Int_t((y-ylowDB)/(yhigDB-ylowDB)*nyDB);
814   if(binX>=nxDB || binY>=nyDB) {
815     binX = -1;
816     binY = -1;
817   }
818   
819 }
820 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
821 Bool_t AliHMPIDReconHTA::UniformDistrib() 
822 {
823   AliHMPIDParam *pParam=AliHMPIDParam::Instance();
824   AliHMPIDRecon pRec;
825
826   Double_t sizeCh = 0.5*pParam->RadThick()+pParam->WinThick()+pParam->GapThick();
827   Double_t xrad = MipX() - sizeCh*TMath::Tan(fThTrkFit)*TMath::Cos(fPhTrkFit);
828   Double_t yrad = MipY() - sizeCh*TMath::Tan(fThTrkFit)*TMath::Sin(fPhTrkFit);
829   
830   pRec.SetTrack(xrad,yrad,fThTrkFit,fPhTrkFit);
831     
832   Int_t npeff=0;
833   Int_t nPhotPerBin = 4;
834   for(Int_t i=0;i<fNClu;i++) {  
835     if(!ClCk(i)) continue;
836     npeff++;
837   }
838   
839   Int_t nPhiBins = npeff/nPhotPerBin+1;
840   if(nPhiBins<=1) nPhiBins = 2;
841
842   Double_t *iPhiBin;
843   iPhiBin = new Double_t[nPhiBins];
844   
845   for(Int_t i=0;i<nPhiBins;i++) iPhiBin[i] =0;
846
847   for(Int_t i=0;i<fNClu;i++) { 
848     if(!ClCk(i)) continue;
849     Double_t phiCer = PhotPhi(i);
850     
851     Double_t PhiProva = phiCer*TMath::RadToDeg();
852     if(PhiProva<0) PhiProva+= 360;
853     Int_t index = (Int_t)((Float_t)nPhiBins*PhiProva/360.);
854     iPhiBin[index]++;
855    }
856
857    Double_t chi2 = 0;
858    for(Int_t i=0;i<nPhiBins;i++) {
859      Double_t theo = (Double_t)npeff/(Double_t)nPhiBins;
860      chi2+= (iPhiBin[i] - theo)*(iPhiBin[i] - theo)/theo;
861    }
862    
863     delete iPhiBin;
864     
865     Double_t prob = TMath::Prob(chi2, nPhiBins-1);
866     AliDebug(1,Form(" Probability for uniform distrib: %6f.3 %s",prob,(prob<0.05) ? "rejected" : "accepted"));
867     if(prob<0.05) return kFALSE;
868     return kTRUE;
869     
870  }
871 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++