MIP cluster size added (compacted word) and infos re-arranganged in HMPID
[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 <TMinuit.h>         //FitFree()
28 #include <TClonesArray.h>    //CkovHiddenTrk()
29 #include <AliESDtrack.h>     //CkovHiddenTrk()
30 #include <TH2I.h>            //InitDatabase()
31 #include <TGraph.h>          //ShapeModel()
32 #include <TSpline.h>         //ShapeModel()
33 #include "TStopwatch.h"      //
34
35 TH2I* AliHMPIDReconHTA::fgDatabase = 0x0;
36
37 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
38 AliHMPIDReconHTA::AliHMPIDReconHTA():
39   TTask("RichRec","RichPat"),
40   fMipX(-999),
41   fMipY(-999),
42   fMipQ(-999),
43   fRadX(-999),
44   fRadY(-999),
45   fIdxMip(0),
46   fNClu(0),
47   fXClu(0),
48   fYClu(0),
49   fClCk(0),
50   fThTrkFit(-999),
51   fPhTrkFit(-999),
52   fCkovFit(-999),
53   fCkovSig2(0),
54   fParam(AliHMPIDParam::Instance())
55 {
56 //..
57 //hidden algorithm
58 //..
59   fParam->SetRefIdx(fParam->MeanIdxRad()); // initialization of ref index to a default one
60   if(!fgDatabase) InitDatabase();
61 }
62 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
63 AliHMPIDReconHTA::~AliHMPIDReconHTA()
64 {
65   DeleteVars();
66 }
67 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
68 void AliHMPIDReconHTA::InitVars(Int_t n)
69 {
70 //..
71 //Init some variables
72 //..
73   fXClu = new Double_t[n];
74   fYClu = new Double_t[n];
75   fClCk = new Bool_t[n];
76   for(Int_t i=0;i<n;i++) fClCk[i] = kTRUE;
77 //
78 }
79 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
80 void AliHMPIDReconHTA::DeleteVars()const
81 {
82 //..
83 //Delete variables
84 //..
85   if(fXClu) delete fXClu;
86   if(fYClu) delete fYClu;
87   if(fClCk) delete fClCk;
88 }
89 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
90 Bool_t AliHMPIDReconHTA::CkovHiddenTrk(AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t nmean, Double_t qthre)
91 {
92 // Pattern recognition method without any infos from tracking:HTA (Hidden Track Algorithm)...
93 // The method finds in the chmuber the cluster with the highest charge
94 // compatibile with a MIP, then the strategy is applied
95 // Arguments:  pTrk     - pointer to ESD track
96 //             pCluLs   - list of clusters for a given chamber 
97 //             nmean    - mean freon ref. index
98 //   Returns:           - 0=ok,1=not fitted 
99   
100   AliHMPIDParam *pParam = AliHMPIDParam::Instance(); 
101   pParam->SetRefIdx(nmean);
102
103   if(!CluPreFilter(pCluLst)) {return kFALSE;}
104   
105   Float_t mipX=-1,mipY=-1;Int_t mipId=-1,mipQ=-1;                                                                           
106   Double_t qRef = 0;
107   Int_t nCh=0;
108   Int_t sizeClu=0;
109   for (Int_t iClu=0;iClu<pCluLst->GetEntriesFast();iClu++){                                   //clusters loop
110     AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(iClu);                       //get pointer to current cluster    
111     nCh = pClu->Ch();
112     fXClu[iClu] = pClu->X();fYClu[iClu] = pClu->Y();                                          //store x,y for fitting procedure
113     fClCk[iClu] = kTRUE;                                                                      //all cluster are accepted at this stage to be reconstructed
114     if(pClu->Q()>qthre) fClCk[iClu] = kFALSE;                                                 // not a good photon in any case (multiple MIPs)
115     if(pClu->Q()>qRef){                                                                       //searching the highest charge to select a MIP      
116       qRef = pClu->Q();
117       mipId=iClu; mipX=pClu->X();mipY=pClu->Y();mipQ=(Int_t)pClu->Q();
118       sizeClu=pClu->Size();
119     }                                                                                    
120 //    Printf(" n. %d x %f y %f Q %f",iClu,pClu->X(),pClu->Y(),pClu->Q());
121   }//clusters loop
122
123   fNClu = pCluLst->GetEntriesFast();
124   if(qRef>qthre){                                                                     //charge compartible with MIP clusters
125     fIdxMip = mipId;
126     fClCk[mipId] = kFALSE;
127     fMipX = mipX; fMipY=mipY; fMipQ = qRef;
128 //    Printf(" mipId %d x %f y %f Q %f",fIdxMip,fMipX,fMipY,fMipQ);
129     if(!DoRecHiddenTrk()) {
130       pTrk->SetHMPIDsignal(pParam->kNoPhotAccept);
131       return kFALSE;
132     }                                                                           //Do track and ring reconstruction,if problems returns 1
133     pTrk->SetHMPIDtrk(fRadX,fRadY,fThTrkFit,fPhTrkFit);                                        //store track intersection info
134     pTrk->SetHMPIDmip(fMipX,fMipY,(Int_t)fMipQ,fNClu);                                         //store mip info 
135     pTrk->SetHMPIDcluIdx(nCh,fIdxMip+1000*sizeClu);                                            //set cham number, index of cluster + cluster size
136     pTrk->SetHMPIDsignal(fCkovFit);                                                            //find best Theta ckov for ring i.e. track
137     pTrk->SetHMPIDchi2(fCkovSig2);                                                             //errors squared
138 //    Printf(" n clusters tot %i accepted %i",pCluLst->GetEntriesFast(),fNClu);
139 //    Printf("CkovHiddenTrk: thetaC %f th %f ph %f",fCkovFit,fThTrkFit,fPhTrkFit);
140     return kTRUE;
141   }
142   
143   return kFALSE;
144 }//CkovHiddenTrk()
145 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
146 Bool_t AliHMPIDReconHTA::DoRecHiddenTrk()
147 {
148 // Pattern recognition method without any infos from tracking...
149 // First a preclustering filter to avoid part of the noise
150 // Then only ellipsed-rings are fitted (no possibility, 
151 // for the moment, to reconstruct very inclined tracks)
152 // Finally a fitting with (th,ph) free, starting by very close values
153 // previously evaluated.
154 // Arguments:   none
155 //   Returns:   none
156   Double_t thTrkRec,phiTrkRec,thetaCRec;
157   
158   if(!FindShape(thTrkRec,phiTrkRec,thetaCRec)) return kFALSE;
159
160   if(!FitFree(thTrkRec,phiTrkRec)) {return kFALSE;}
161   
162   return kTRUE;
163 }//DoRecHiddenTrk()
164 /*
165 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
166 Bool_t AliHMPIDReconHTA::DoRecHiddenTrk(TClonesArray *pCluLst)
167 {
168 // Pattern recognition method without any infos from tracking...
169 // First a preclustering filter to avoid part of the noise
170 // Then only ellipsed-rings are fitted (no possibility, 
171 // for the moment, to reconstruct very inclined tracks)
172 // Finally a fitting with (th,ph) free, starting by very close values
173 // previously evaluated.
174 // Arguments:   none
175 //   Returns:   none
176   Double_t thTrkRec,phiTrkRec,thetaCRec;
177   
178   if(!FindShape(thTrkRec,phiTrkRec,thetaCRec)) return kFALSE;
179
180 //  Printf("thTrkRec %f phiTrkRec %f ThetaCRec %f",thTrkRec*TMath::RadToDeg(),phiTrkRec*TMath::RadToDeg(),thetaCRec*TMath::RadToDeg());  
181   Int_t nClTmp1 = pCluLst->GetEntriesFast()-1;  //minus MIP...
182   Int_t nClTmp2 = 0;
183   
184   while(nClTmp1 != nClTmp2){
185     SetNClu(pCluLst->GetEntriesFast());
186     if(!FitFree(thTrkRec,phiTrkRec)) {return kFALSE;}
187     nClTmp2 = NClu();
188     if(nClTmp2!=nClTmp1) {nClTmp1=nClTmp2;nClTmp2=0;}
189   }
190   
191   fNClu = nClTmp2;
192   return kTRUE;
193 }//DoRecHiddenTrk()
194 */
195 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
196 Bool_t AliHMPIDReconHTA::CluPreFilter(TClonesArray *pCluLst)
197 {
198 // Pre-filter of bkg clusters
199 // Arguments:    pSluLst  -  List of the clusters for a given chamber
200 //   Returns:    status   -  TRUE if filtering leaves enough photons, FALSE if not
201 //
202   Int_t nClusTot = pCluLst->GetEntriesFast();
203   if(nClusTot<4||nClusTot>100) {
204     return kFALSE; 
205   } else { 
206     InitVars(nClusTot);
207     return kTRUE;
208   }
209 }
210 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
211 Bool_t AliHMPIDReconHTA::FindShape(Double_t &thTrkRec,Double_t &phiTrkRec,Double_t &thetaCRec)
212 {
213 // Finds the estimates for phi and theta of the track and the ThetaCerenkov
214 // by using a database of the shapes of the rings
215 // Arguments:   none  
216 //   Returns:   thTrkRec  - estimate of theta track
217 //              phiTrkRec - estimate of phi   track
218 //              thetaCRec - estimate of ThetaCerenkov
219 //              status    - TRUE if a good solution is found, FALSE if not
220
221   Double_t *phiphot = new Double_t[fNClu];  
222   Double_t *dist    = new Double_t[fNClu];  
223   Int_t    *indphi  = new    Int_t[fNClu];  
224
225   Bool_t status;
226     
227 // Sort in phi angle...  
228   for(Int_t i=0;i<fNClu;i++) {
229     phiphot[i] = (TMath::ATan2(fMipY-fYClu[i],fMipX-fXClu[i])+TMath::Pi())*TMath::RadToDeg();
230     dist[i]=TMath::Sqrt((fMipX-fXClu[i])*(fMipX-fXClu[i])+(fMipY-fYClu[i])*(fMipY-fYClu[i]));
231   }
232   
233   TMath::Sort(fNClu,phiphot,indphi,kFALSE);
234   
235 // Purify with a truncated mean;
236   Int_t np=0;
237   Double_t dMean  = 0;
238   Double_t dMean2 = 0;
239   for(Int_t i=0;i<fNClu;i++) {
240     if(!fClCk[indphi[i]]) continue;                                                  // Check if a good photon candidate or not
241     dMean +=dist[indphi[i]];
242     dMean2+=dist[indphi[i]]*dist[indphi[i]];
243     np++;
244   }
245   
246   dMean  /=(Double_t)np;
247   dMean2 /=(Double_t)np;
248   Double_t rms = TMath::Sqrt(dMean2 - dMean*dMean);
249   
250   for(Int_t i=0;i<fNClu;i++) {
251     if(!fClCk[indphi[i]]) continue;                                                  // Check if a good photon candidate or not
252     if(TMath::Abs(dMean-dist[indphi[i]]) > 2*rms) {
253       fClCk[indphi[i]] = kFALSE;
254       continue;
255     }
256   }
257
258 //
259 //  purify vectors for good photon candidates
260 //
261   Int_t npeff=0;
262   Double_t *phiphotP = new Double_t[fNClu+1];  
263   Double_t *distP    = new Double_t[fNClu+1];  
264   for(Int_t i=0;i<fNClu;i++) {
265     if(!fClCk[indphi[i]]) continue;                                                  // Check if a good photon candidate or not
266     phiphotP[npeff] = phiphot[indphi[i]];
267     distP[npeff]    = dist[indphi[i]];
268     npeff++;
269   }
270   
271   delete [] phiphot;
272   delete [] dist;
273   delete [] indphi;
274
275   if(npeff<3) {
276     delete [] phiphotP;
277     delete [] distP;
278     return kFALSE;
279   }
280
281 //  for(Int_t i=0;i<npeff;i++) {Printf(" n. %d phiphot %f dist %f",i,phiphotP[i],distP[i]);}
282   
283   Double_t xA,xB;
284   if(ShapeModel(npeff,phiphotP,distP,xA,xB,phiTrkRec)) {
285     
286 //    Printf("FindShape: phi start %f",phiTrkRec*TMath::RadToDeg());
287
288     Int_t bin = fgDatabase->FindBin(xA,xB);
289     Int_t compact = (Int_t)fgDatabase->GetBinContent(bin);
290     thetaCRec =  (Double_t)(compact%1000);
291     thTrkRec  =  (Double_t)(compact/1000);
292
293     thTrkRec *= TMath::DegToRad(); 
294     thetaCRec *= TMath::DegToRad();
295
296 //    Printf("FindShape: xA %f xB %f compact %d thTrk %f thC %f",xA,xB,compact,thTrkRec*TMath::RadToDeg(),thetaCRec*TMath::RadToDeg());
297         
298     status = kTRUE;
299     
300   } else {
301     
302     status = kFALSE;
303     
304   }
305
306   delete [] phiphotP;
307   delete [] distP;
308   
309   return status;
310 }
311 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
312 Bool_t AliHMPIDReconHTA::ShapeModel(Int_t np,Double_t *phiphot,Double_t *dist,Double_t &xA,Double_t &xB,Double_t &phiStart)
313 {
314 // Find a Spline curve to define dist. vs. phi angle
315 // in order to estimate the phi of the track
316 // Arguments:   np     - # points corresponding to # photon candidates
317 //             dist    - distance of each photon from MIP
318 //             phiphot - phi of the photon in the DRS
319 //   Returns:  xA      - min. distance from MIP
320 //             xB      - dist. from mip perpedicular to the major axis 
321 //             phiStart- estimate of the track phi
322
323   TGraph *phigr = new TGraph(np,phiphot,dist);
324   TSpline3 *sphi = new TSpline3("sphi",phigr);
325   if(!sphi) {Printf("Spline not created!Bye.");return kFALSE;}
326
327   Int_t locMin = TMath::LocMin(np,dist);
328   Int_t locMax = TMath::LocMax(np,dist);
329   
330   Double_t minX = phiphot[locMin];
331 //  Double_t minY =    dist[locMin];
332   Double_t maxX = phiphot[locMax];
333 //  Double_t maxY =    dist[locMax];
334   
335   Int_t ip[3] = {-1,0,1};
336   if(locMin==0   ) {ip[0]= 0;ip[1]= 1;ip[2]= 2;}
337   if(locMin==np-1) {ip[0]=-2;ip[1]=-1;ip[2]= 0;}
338   
339   Double_t minXf = VertParab(phiphot[locMin+ip[0]],dist[locMin+ip[0]], 
340                              phiphot[locMin+ip[1]],dist[locMin+ip[1]],
341                              phiphot[locMin+ip[2]],dist[locMin+ip[2]]);
342   if(minXf< phiphot[locMin+ip[0]] || minXf > phiphot[locMin+ip[2]]) minXf = minX;
343   
344   ip[0]=-1;ip[1]=0;ip[2]=1;
345   if(locMax==0   ) {ip[0]= 0;ip[1]= 1;ip[2]= 2;}
346   if(locMax==np-1) {ip[0]=-2;ip[1]=-1;ip[2]= 0;}
347   
348   Double_t maxXf = VertParab(phiphot[locMax+ip[0]],dist[locMax+ip[0]], 
349                              phiphot[locMax+ip[1]],dist[locMax+ip[1]],
350                              phiphot[locMax+ip[2]],dist[locMax+ip[2]]);
351   if(maxXf< phiphot[locMax+ip[0]] || maxXf > phiphot[locMax+ip[2]]) maxXf = maxX;
352   
353 //  Printf(" phi at mindist %f and found %f",minX,minXf);
354 //  Printf(" phi at maxdist %f and found %f",maxX,maxXf);
355 //  
356   if(TMath::Abs(maxXf-minXf)>30) {
357     xA = sphi->Eval(minXf);
358     if(xA < 0) xA = dist[sphi->FindX(xA)];
359     xB = sphi->Eval(minXf-90);
360     if(xB < 0) xB = dist[sphi->FindX(xB)];
361     phiStart = minXf-180;  //open ring or acceptance effect...so believe to min phi angle!
362   } else {
363     phiStart = 0.5*(maxXf-180+minXf);
364     xA = sphi->Eval(phiStart+180);
365     if(xA < 0) xA = dist[sphi->FindX(xA)];
366     xB = sphi->Eval(phiStart+90);
367     if(xB < 0) xB = dist[sphi->FindX(xB)];
368   }
369   //
370 //  Printf("ShapeModel: phiStart %f xA %f xB %f",phiStart,xA,xB);
371   
372   phiStart*=TMath::DegToRad();
373   return kTRUE;
374 }
375 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
376 Double_t AliHMPIDReconHTA::VertParab(Double_t x1,Double_t y1,Double_t x2, Double_t y2, Double_t x3, Double_t y3)const
377 {
378 // It uses parabola from 3 points to evaluate the x-coord of the parab 
379 // Arguments:    xi,yi - points
380 //   Returns:    x-coord of the vertex 
381   
382   Double_t a = ((x1-x3)*(y1-y2)-(x1-x2)*(y1-y3))/((x1*x1-x2*x2)*(x1-x3)-(x1*x1-x3*x3)*(x1-x2));
383   Double_t b = (y1-y2 - a*(x1*x1-x2*x2))/(x1-x2);
384 //  Double_t c = y1 - a*x1*x1-b*x1;
385   return -0.5*b/a;
386 }
387 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
388 Bool_t AliHMPIDReconHTA::FitFree(Double_t thTrkRec,Double_t phiTrkRec)
389 {
390 // Fit performed by minimizing RMS/sqrt(n) of the
391 // photons reconstructed. First phi is fixed and theta
392 // is fouond, then (th,ph) of the track
393 // as free parameters
394 // Arguments:    PhiRec phi of the track
395 //   Returns:    none
396   
397   TMinuit *pMinuit = new TMinuit(2);
398   pMinuit->mncler();
399   gMinuit->SetObjectFit((TObject*)this);  gMinuit->SetFCN(AliHMPIDReconHTA::FunMinPhot);  //set fit function
400   Double_t aArg=-1,parStep,parLow,parHigh;     Int_t iErrFlg;                 //tmp vars for TMinuit
401   Double_t d1,d2,d3;
402   TString sName;
403   Double_t th,ph;
404   
405   pMinuit->mnexcm("SET PRI",&aArg,1,iErrFlg);                                          //suspend all printout from TMinuit
406   pMinuit->mnexcm("SET NOW",&aArg,0,iErrFlg);
407
408   if(thTrkRec==0) thTrkRec = 3.*TMath::DegToRad();    // not to start from the edge...
409   
410   pMinuit->mnparm(0," thTrk  ",thTrkRec ,parStep=0.01,parLow=0,parHigh=TMath::PiOver4(),iErrFlg);
411   pMinuit->mnparm(1," phiTrk ",phiTrkRec,parStep=0.01,parLow=0,parHigh=TMath::TwoPi(),iErrFlg);
412   
413   pMinuit->FixParameter(1);
414   pMinuit->mnexcm("SIMPLEX" ,&aArg,0,iErrFlg);   
415   pMinuit->mnexcm("MIGRAD"  ,&aArg,0,iErrFlg);
416   pMinuit->Release(1);  
417   pMinuit->mnexcm("MIGRAD"  ,&aArg,0,iErrFlg);
418   
419   pMinuit->mnpout(0,sName,th,d1,d2,d3,iErrFlg);
420   pMinuit->mnpout(1,sName,ph,d1,d2,d3,iErrFlg);   
421   
422   Double_t f,par[2];
423   Double_t *grad=0x0;
424   par[0] = th;par[1] = ph;
425   pMinuit->Eval(2,grad,f,par,3);
426
427 //  Printf("FitFree: theta %f phi %f",th,ph);
428   
429   SetTrkFit(th,ph);
430   return kTRUE;
431 }
432 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
433 void AliHMPIDReconHTA::FunMinPhot(Int_t &/* */,Double_t* /* */,Double_t &f,Double_t *par,Int_t iflag)
434 {
435 // Minimization function to find best track and thetaC parameters
436 // Arguments:    f = function value to minimize
437 //             par = list of parameter to find
438 //           iflag = flag status. See Minuit instructions
439 //   Returns:    none
440 //
441 // Note: it is necessary to call an instance of AlihMPIDParam. Not possible to use fParam
442 // because of the static instantiation of the function in Minuit
443   
444   AliHMPIDParam *pParam=AliHMPIDParam::Instance();
445   AliHMPIDReconHTA *pRecHTA=(AliHMPIDReconHTA*)gMinuit->GetObjectFit();
446   AliHMPIDRecon pRec;
447   Double_t sizeCh = 0.5*pParam->RadThick()+pParam->WinThick()+pParam->GapThick();
448   Double_t thTrk = par[0]; 
449   Double_t phTrk = par[1];
450   Double_t xrad = pRecHTA->MipX() - sizeCh*TMath::Tan(thTrk)*TMath::Cos(phTrk);
451   Double_t yrad = pRecHTA->MipY() - sizeCh*TMath::Tan(thTrk)*TMath::Sin(phTrk);
452   pRecHTA->SetRadXY(xrad,yrad);
453   pRec.SetTrack(xrad,yrad,thTrk,phTrk);
454
455   Double_t meanCkov =0;
456   Double_t meanCkov2=0;
457   Double_t thetaCer,phiCer;
458   Int_t nClAcc = 0;
459   Int_t nClTot=pRecHTA->NClu();
460     
461   for(Int_t i=0;i<nClTot;i++) {
462     if(!(pRecHTA->ClCk(i))) continue;
463     pRec.FindPhotCkov(pRecHTA->XClu(i),pRecHTA->YClu(i),thetaCer,phiCer);  
464     meanCkov  += thetaCer;
465     meanCkov2 += thetaCer*thetaCer;
466     nClAcc++;
467   }
468   if(nClAcc==0) {f=999;return;}
469   meanCkov  /=(Double_t)nClAcc;
470   meanCkov2 /=(Double_t)nClAcc;
471   Double_t rms = TMath::Sqrt(TMath::Abs(meanCkov2 - meanCkov*meanCkov));
472   f = rms/TMath::Sqrt((Double_t)nClAcc);
473   
474   if(iflag==3) {
475 /*
476     Printf("FunMinPhot before: photons candidates %i used %i",nClTot,nClAcc);
477     nClAcc = 0;
478     Double_t meanCkov1=0;
479     Double_t meanCkov2=0;
480     for(Int_t i=0;i<nClTot;i++) {
481       if(!(pRec->ClCk(i))) continue;
482       pRec->FindPhotCkov(pRec->XClu(i),pRec->YClu(i),thetaCer,phiCer);  
483       if(TMath::Abs(thetaCer-meanCkov)<2*rms) {
484         meanCkov1 += thetaCer;
485         meanCkov2 += thetaCer*thetaCer;
486         nClAcc++;
487       } else pRec->SetClCk(i,kFALSE);
488     }
489     meanCkov1/=nClAcc;
490     Double_t rms2 = (meanCkov2 - meanCkov*meanCkov*nClAcc)/nClAcc;
491     Printf("FunMinPhot after: photons candidates %i used %i thetaCer %f",nClTot,nClAcc,meanCkov1);
492     pRec->SetCkovFit(meanCkov1);
493     pRec->SetCkovSig2(rms2);
494     pRec->SetNClu(nClAcc);
495 */
496 //    Printf("FunMinPhot: reconstructed theta Cerenkov %f with %d photons",meanCkov,nClAcc);
497     pRecHTA->SetCkovFit(meanCkov);
498     pRecHTA->SetCkovSig2(rms*rms);
499     pRecHTA->SetNClu(nClAcc);
500   }
501 }//FunMinPhot()
502 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
503 void AliHMPIDReconHTA::InitDatabase()
504 {
505 // Construction a database of ring shapes on fly
506 //   Arguments: none
507 //   Returns  : none
508 //  N.B. fgDatabase points to a TH2I with x-min dist from MIP
509 //                                        y-dist from the ring of the MIP perpendicular to major axis
510 //        The content is the packed info of track theta and thetaC in degrees
511 //                        thetaC+1000*thTrk
512 //
513   TStopwatch timer;
514   timer.Start();
515   
516   AliInfo(Form("database HTA is being built.Please, wait..."));
517 //  
518   Double_t x[3]={0,0,0},y[3];
519
520   AliHMPIDRecon rec;
521
522   if(!fParam) fParam=AliHMPIDParam::Instance();
523   Double_t thetaMax = TMath::ACos(1./fParam->MeanIdxRad());
524   Double_t thTrkMax = 1./TMath::ASin(fParam->MeanIdxRad());    
525
526   Int_t nstepx = 1000;
527   Int_t nstepy = 1000;
528
529   TH2I *deconv = new TH2I("deconv","database;d1;d2;thC+1000*thTrk",500,0,50,150,0,15);
530   //
531   Double_t xrad = 0;
532   Double_t yrad = 0;
533   Double_t phTrk = 0;
534
535   for(Int_t i=0;i<nstepx;i++) {     //loop on thetaC
536     for(Int_t j=0;j<nstepy;j++) {   //loop on theta particle
537       Double_t thetaC = thetaMax/nstepx*((Double_t)i+0.5);
538       Double_t thTrk  = thTrkMax/nstepy*((Double_t)j+0.5);
539       //
540       //mip position
541       //
542       Double_t sizeCh = 0.5*fParam->RadThick()+fParam->WinThick()+fParam->GapThick();
543       Double_t xmip = xrad + sizeCh*TMath::Tan(thTrk)*TMath::Cos(phTrk);
544       Double_t ymip = yrad + sizeCh*TMath::Tan(thTrk)*TMath::Sin(phTrk);
545
546       Double_t dist1,dist2;
547       //
548       //first point at phi=0
549       //
550       rec.SetTrack(xrad,yrad,thTrk,phTrk);
551       TVector2 pos;
552       pos=rec.TracePhot(thetaC,0);
553
554       if(pos.X()==-999) {
555         dist1 = 0;            //open ring...anly the distance btw mip and point at 180 will be considered
556       } else {
557         x[0] = pos.X(); y[0] = pos.Y();
558         dist1   = TMath::Sqrt((x[0]-xmip)*(x[0]-xmip)+(y[0]-ymip)*(y[0]-ymip));
559       }
560       //
561       //second point at phi=180
562       //
563       rec.SetTrack(xrad,yrad,thTrk,phTrk);
564       pos=rec.TracePhot(thetaC,TMath::Pi());
565
566       if(pos.X()==-999) {Printf("it should not happens!Bye");return;}
567       x[1] = pos.X(); y[1] = pos.Y();
568       if((x[1]-xmip)*(x[0]-xmip)>0) continue; // to avoid circles out mips (for very low ThetaC)
569       dist2   = TMath::Sqrt((x[1]-xmip)*(x[1]-xmip)+(y[1]-ymip)*(y[1]-ymip));
570
571 //      Double_t distA = dist1+dist2;
572       Double_t distA = dist2;     // only the minimum: problem of acceptance
573       //
574       //second point at phi=90
575       //
576       rec.SetTrack(xrad,yrad,thTrk,phTrk);
577       pos=rec.TracePhot(thetaC,TMath::PiOver2());
578
579       if(pos.X()==-999) continue;
580       x[2] = pos.X(); y[2] = pos.Y();
581       Double_t distB   = TMath::Sqrt((x[2]-xmip)*(x[2]-xmip)+(y[2]-ymip)*(y[2]-ymip));
582 // compact the infos...      
583       Int_t compact = (Int_t)(thetaC*TMath::RadToDeg())+1000*(Int_t)(thTrk*TMath::RadToDeg());
584       Int_t bin = deconv->FindBin(distA,distB);
585       if(deconv->GetBinContent(bin)==0) deconv->Fill(distA,distB,compact);
586     }
587   }
588
589   FillZeroChan(deconv);
590   fgDatabase = deconv;
591
592   timer.Stop();
593   Double_t nSecs = timer.CpuTime();  
594   AliInfo(Form("database HTA successfully open in %3.1f sec.(CPU). Reconstruction is started.",nSecs));
595 }
596 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
597 void AliHMPIDReconHTA::FillZeroChan(TH2I *deconv)const
598 {
599   //If fills eventually channel without entries
600   //inthe histo "database" jyst interpolating the neighboring cells
601   // Arguments: histogram pointer of the database
602   //   Returns: none
603   //
604   Int_t nbinx = deconv->GetNbinsX();
605   Int_t nbiny = deconv->GetNbinsY();
606   for(Int_t i = 0;i<nbinx;i++) {
607     for(Int_t j = 0;j<nbiny;j++) {
608       if(deconv->GetBinContent(i,j) == 0) {
609         Int_t nXmin = i-1; Int_t nXmax=i+1;
610         Int_t nYmin = j-1; Int_t nYmax=j+1;
611         Int_t nc = 0;
612         Double_t meanC =0;
613         Double_t meanTrk =0;
614         for(Int_t ix=nXmin;ix<=nXmax;ix++) {
615           if(ix<0||ix>nbinx) continue;
616           for(Int_t iy=nYmin;iy<=nYmax;iy++) {
617             if(iy<0||iy>nbiny) continue;
618             meanC  +=  (Int_t)deconv->GetBinContent(ix,iy)%1000;
619             meanTrk+=  (Int_t)deconv->GetBinContent(ix,iy)/1000;
620             nc++;
621           }
622           meanC/=nc; meanTrk/=nc;
623           Int_t compact = (Int_t)meanC+1000*(Int_t)meanTrk;
624           if(compact>0)deconv->SetCellContent(i,j,compact);
625         }
626       }
627     }
628   }
629 }
630 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++