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