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 "TStopwatch.h"      //
  35
  36 //TH2F* AliHMPIDReconHTA::fgDatabase = 0x0;
  37 TH2F* AliHMPIDReconHTA::fgDatabase = new TH2F("deconv","database;d1;d2;thC+1000*thTrk",500,0,50,150,0,15);
  38 Int_t AliHMPIDReconHTA::fgDB[501][51]={25551*0};
  39 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  40 AliHMPIDReconHTA::AliHMPIDReconHTA():
  41   TTask("RichRec","RichPat"),
  42   fMipX(-999),
  43   fMipY(-999),
  44   fMipQ(-999),
  45   fRadX(-999),
  46   fRadY(-999),
  47   fIdxMip(0),
  48   fNClu(0),
  49   fXClu(0),
  50   fYClu(0),
  51   fClCk(0),
  52   fThTrkIn(-999),
  53   fPhTrkIn(-999),
  54   fThTrkFit(-999),
  55   fPhTrkFit(-999),
  56   fCkovFit(-999),
  57   fCkovSig2(0),
  58   fParam(AliHMPIDParam::Instance())
  59 {
  60 //..
  61 //hidden algorithm
  62 //..
  63   fParam->SetRefIdx(fParam->MeanIdxRad()); // initialization of ref index to a default one
  64   InitDatabase();
  65 }
  66 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  67 AliHMPIDReconHTA::~AliHMPIDReconHTA()
  68 {
  69   DeleteVars();
  70 }
  71 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  72 void AliHMPIDReconHTA::InitVars(Int_t n)
  73 {
  74 //..
  75 //Init some variables
  76 //..
  77   fXClu = new Double_t[n];
  78   fYClu = new Double_t[n];
  79   fClCk = new Bool_t[n];
  80   for(Int_t i=0;i<n;i++) fClCk[i] = kTRUE;
  81 //
  82 }
  83 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  84 void AliHMPIDReconHTA::DeleteVars()const
  85 {
  86 //..
  87 //Delete variables
  88 //..
  89   if(fXClu) delete fXClu;
  90   if(fYClu) delete fYClu;
  91   if(fClCk) delete fClCk;
  92 }
  93 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  94 Bool_t AliHMPIDReconHTA::CkovHiddenTrk(AliESDtrack *pTrk,TClonesArray *pCluLst,Int_t index, Double_t nmean)
  95 {
  96 // Pattern recognition method without any infos from tracking:HTA (Hidden Track Algorithm)...
  97 // The method finds in the chmuber the cluster with the highest charge
  98 // compatibile with a MIP, then the strategy is applied
  99 // Arguments:  pTrk     - pointer to ESD track
 100 //             pCluLs   - list of clusters for a given chamber
 101 //             pNmean   - pointer to ref. index
 102 //             pQthre   - pointer to qthre
 103 //   Returns:           - 0=ok,1=not fitted
 104
 105   AliHMPIDParam *pParam = AliHMPIDParam::Instance();
 106
 107   if(!CluPreFilter(pCluLst)) return kFALSE;
 108
 109   Int_t nCh=0;
 110   Int_t sizeClu=0;
 111
 112   fNClu = pCluLst->GetEntriesFast();
 113
 114   for (Int_t iClu=0;iClu<fNClu;iClu++){                                                       //clusters loop
 115     AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(iClu);                       //get pointer to current cluster
 116     fXClu[iClu] = pClu->X();fYClu[iClu] = pClu->Y();                                          //store x,y for fitting procedure
 117     fClCk[iClu] = kTRUE;                                                                      //all cluster are accepted at this stage to be reconstructed
 118
 119     if(iClu == index) {
 120
 121       fMipX = pClu->X();
 122       fMipY = pClu->Y();
 123       fMipQ = pClu->Q();
 124       sizeClu = pClu->Size();
 125       nCh = pClu->Ch();
 126       fClCk[index] = kFALSE;
 127       fIdxMip = index;
 128  //    Printf(" n. %d x %f y %f Q %f",iClu,pClu->X(),pClu->Y(),pClu->Q());
 129     }
 130   }//clusters loop
 131
 132   pParam->SetRefIdx(nmean);
 133
 134   //
 135   Float_t xra,yra,th,ph; pTrk->GetHMPIDtrk(xra,yra,th,ph);
 136 //  Printf(" simulated phi %6.2f ",ph*TMath::RadToDeg());
 137   //
 138
 139   if(!DoRecHiddenTrk()) {
 140     pTrk->SetHMPIDsignal(pParam->kNoPhotAccept);
 141     return kFALSE;
 142   }                                                                           //Do track and ring reconstruction,if problems returns 1
 143 //  Printf("    fitted phi %6.2f ",fPhTrkFit*TMath::RadToDeg());
 144
 145   pTrk->SetHMPIDtrk(fRadX,fRadY,fThTrkFit,fPhTrkFit);                                        //store track intersection info
 146   pTrk->SetHMPIDmip(fMipX,fMipY,(Int_t)fMipQ,fNClu);                                         //store mip info
 147   pTrk->SetHMPIDcluIdx(nCh,fIdxMip+1000*sizeClu);                                            //set cham number, index of cluster + cluster size
 148   pTrk->SetHMPIDsignal(fCkovFit);                                                            //find best Theta ckov for ring i.e. track
 149   pTrk->SetHMPIDchi2(fCkovSig2);                                                             //errors squared
 150 //    Printf(" n clusters tot %i accepted %i",pCluLst->GetEntriesFast(),fNClu);
 151 //    Printf("CkovHiddenTrk: thetaC %f th %f ph %f",fCkovFit,fThTrkFit,fPhTrkFit);
 152
 153   return kTRUE;
 154
 155 }//CkovHiddenTrk()
 156 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 157 Bool_t AliHMPIDReconHTA::DoRecHiddenTrk()
 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)) {
 170 //    Printf("FindShape failed...!");
 171     return kFALSE;
 172   }
 173
 174   if(!FitFree(thTrkRec,phiTrkRec)) {
 175 //    Printf("FitFree failed...!");
 176     return kFALSE;
 177   }
 178
 179   return kTRUE;
 180 }//DoRecHiddenTrk()
 181 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 182 Bool_t AliHMPIDReconHTA::CluPreFilter(TClonesArray *pCluLst)
 183 {
 184 // Pre-filter of bkg clusters
 185 // Arguments:    pSluLst  -  List of the clusters for a given chamber
 186 //   Returns:    status   -  TRUE if filtering leaves enough photons, FALSE if not
 187 //
 188   Int_t nClusTot = pCluLst->GetEntriesFast();
 189   if(nClusTot<4||nClusTot>100) {
 190     return kFALSE;
 191   } else {
 192     InitVars(nClusTot);
 193     return kTRUE;
 194   }
 195 }
 196 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 197 Bool_t AliHMPIDReconHTA::FindShape(Double_t &thTrkRec,Double_t &phiTrkRec,Double_t &thetaCRec)
 198 {
 199 // Finds the estimates for phi and theta of the track and the ThetaCerenkov
 200 // by using a database of the shapes of the rings
 201 // Arguments:   none
 202 //   Returns:   thTrkRec  - estimate of theta track
 203 //              phiTrkRec - estimate of phi   track
 204 //              thetaCRec - estimate of ThetaCerenkov
 205 //              status    - TRUE if a good solution is found, FALSE if not
 206
 207   Double_t *phiphot = new Double_t[fNClu];
 208   Double_t *dist    = new Double_t[fNClu];
 209   Int_t    *indphi  = new    Int_t[fNClu];
 210
 211   Bool_t status;
 212
 213 // Sort in phi angle...
 214 //  Printf(" mipX %f mipy %f",fMipX,fMipY);
 215   for(Int_t i=0;i<fNClu;i++) {
 216     if(!fClCk[i]) {
 217       phiphot[i] = 999.;
 218       dist[i] = 999.;
 219       continue;
 220     }
 221     phiphot[i] = (TMath::ATan2(fMipY-fYClu[i],fMipX-fXClu[i])+TMath::Pi())*TMath::RadToDeg();
 222     dist[i]=TMath::Sqrt((fMipX-fXClu[i])*(fMipX-fXClu[i])+(fMipY-fYClu[i])*(fMipY-fYClu[i]));
 223 //    Printf(" n.%3i  phiphot %f dist %f check %i",i,phiphot[i],dist[i],fClCk[i]);
 224   }
 225
 226   TMath::Sort(fNClu,phiphot,indphi,kFALSE);
 227
 228 // Purify with a truncated mean;
 229   Int_t np=0;
 230   Double_t dMean  = 0;
 231   Double_t dMean2 = 0;
 232   for(Int_t i=0;i<fNClu;i++) {
 233     if(!fClCk[indphi[i]]) continue;                                                  // Check if a good photon candidate or not
 234     dMean +=dist[indphi[i]];
 235     dMean2+=dist[indphi[i]]*dist[indphi[i]];
 236     np++;
 237   }
 238
 239   dMean  /=(Double_t)np;
 240   dMean2 /=(Double_t)np;
 241   Double_t rms = TMath::Sqrt(dMean2 - dMean*dMean);
 242
 243   for(Int_t i=0;i<fNClu;i++) {
 244     if(!fClCk[indphi[i]]) continue;                                                  // Check if a good photon candidate or not
 245     if(TMath::Abs(dMean-dist[indphi[i]]) > 2*rms) {
 246       fClCk[indphi[i]] = kFALSE;
 247       continue;
 248     }
 249   }
 250
 251 //
 252 //  purify vectors for good photon candidates
 253 //
 254   Int_t npeff=0;
 255   Double_t *phiphotP = new Double_t[fNClu+1];
 256   Double_t *distP    = new Double_t[fNClu+1];
 257   for(Int_t i=0;i<fNClu;i++) {
 258     if(!fClCk[indphi[i]]) continue;                                                  // Check if a good photon candidate or not
 259     phiphotP[npeff] = phiphot[indphi[i]];
 260     distP[npeff]    = dist[indphi[i]];
 261 //    Printf("n. %2i phi %f dist %f",npeff,phiphotP[npeff],distP[npeff]);
 262     npeff++;
 263   }
 264
 265   delete [] phiphot;
 266   delete [] dist;
 267   delete [] indphi;
 268
 269   if(npeff<3) {
 270     Printf("FindShape failed: no enough photons = %i...",npeff);
 271     delete [] phiphotP;
 272     delete [] distP;
 273     return kFALSE;
 274   }
 275
 276 //  for(Int_t i=0;i<npeff;i++) {Printf(" n. %d phiphot %f dist %f",i,phiphotP[i],distP[i]);}
 277
 278   Double_t xA,xB;
 279   status = kFALSE;
 280
 281   if(!ShapeModel(npeff,phiphotP,distP,xA,xB,phiTrkRec)) {Printf("ShapeModel failed            "); return kFALSE;}
 282
 283 //  if(xA > 50 || xB > 15)                                {Printf("xA and xB failed out of range"); return kFALSE;}
 284
 285   Int_t bin = fgDatabase->FindBin(xA,xB);
 286   if(bin<=0)                                            {Printf("bin < 0 ! failed             "); return kFALSE;}
 287
 288   Int_t compact = (Int_t)fgDatabase->GetBinContent(bin);
 289
 290
 291   if(compact<0)                                         {Printf("compact< 0! failed           "); return kFALSE;}
 292   //
 293   Int_t binxDB,binyDB;
 294   FindBinDB(xA,xB,binxDB,binyDB);
 295   Int_t compactDB;
 296   if(binxDB>0 && binyDB>0 ) compactDB = fgDB[binxDB][binyDB];
 297   //
 298   Printf("compact %i compactDB %i",compact,compactDB);
 299
 300   thetaCRec =  (Double_t)(compact%1000);
 301   thTrkRec  =  (Double_t)(compact/1000);
 302
 303   thTrkRec  *= TMath::DegToRad();
 304   thetaCRec *= TMath::DegToRad();
 305
 306   status = kTRUE;
 307
 308   delete [] phiphotP;
 309   delete [] distP;
 310
 311   return status;
 312 }
 313 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 314 Bool_t AliHMPIDReconHTA::ShapeModel(Int_t np,Double_t *phiphot,Double_t *dist,Double_t &xA,Double_t &xB,Double_t &phiStart)
 315 {
 316 // Find a Spline curve to define dist. vs. phi angle
 317 // in order to estimate the phi of the track
 318 // Arguments:   np     - # points corresponding to # photon candidates
 319 //             dist    - distance of each photon from MIP
 320 //             phiphot - phi of the photon in the DRS
 321 //   Returns:  xA      - min. distance from MIP
 322 //             xB      - dist. from mip perpedicular to the major axis
 323 //             phiStart- estimate of the track phi
 324
 325   TGraph *phigr = new TGraph(np,phiphot,dist);
 326   TSpline3 *sphi = new TSpline3("sphi",phigr);
 327   if(!sphi) {Printf("Spline not created!Bye.");return kFALSE;}
 328
 329   Int_t locMin = TMath::LocMin(np,dist);
 330   Int_t locMax = TMath::LocMax(np,dist);
 331
 332   Double_t minX = phiphot[locMin];
 333 //  Double_t minY =    dist[locMin];
 334   Double_t maxX = phiphot[locMax];
 335 //  Double_t maxY =    dist[locMax];
 336
 337   Int_t ip[3] = {-1,0,1};
 338   if(locMin==0   ) {ip[0]= 0;ip[1]= 1;ip[2]= 2;}
 339   if(locMin==np-1) {ip[0]=-2;ip[1]=-1;ip[2]= 0;}
 340
 341   Double_t minXf = VertParab(phiphot[locMin+ip[0]],dist[locMin+ip[0]],
 342                              phiphot[locMin+ip[1]],dist[locMin+ip[1]],
 343                              phiphot[locMin+ip[2]],dist[locMin+ip[2]]);
 344   if(minXf< phiphot[locMin+ip[0]] || minXf > phiphot[locMin+ip[2]]) minXf = minX;
 345
 346   ip[0]=-1;ip[1]=0;ip[2]=1;
 347   if(locMax==0   ) {ip[0]= 0;ip[1]= 1;ip[2]= 2;}
 348   if(locMax==np-1) {ip[0]=-2;ip[1]=-1;ip[2]= 0;}
 349
 350   Double_t maxXf = VertParab(phiphot[locMax+ip[0]],dist[locMax+ip[0]],
 351                              phiphot[locMax+ip[1]],dist[locMax+ip[1]],
 352                              phiphot[locMax+ip[2]],dist[locMax+ip[2]]);
 353   if(maxXf< phiphot[locMax+ip[0]] || maxXf > phiphot[locMax+ip[2]]) maxXf = maxX;
 354
 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
 373   phiStart*=TMath::DegToRad();
 374   //----
 375   Double_t phitest = FindSimmPhi();
 376   phiStart = phitest;
 377   //---
 378 //  Printf("phiStart %f phiTest %f",phiStart*TMath::RadToDeg(),phitest*TMath::RadToDeg());
 379
 380   return kTRUE;
 381 }
 382 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 383 Double_t AliHMPIDReconHTA::VertParab(Double_t x1,Double_t y1,Double_t x2, Double_t y2, Double_t x3, Double_t y3)const
 384 {
 385 // It uses parabola from 3 points to evaluate the x-coord of the parab
 386 // Arguments:    xi,yi - points
 387 //   Returns:    x-coord of the vertex
 388
 389   Double_t a = ((x1-x3)*(y1-y2)-(x1-x2)*(y1-y3))/((x1*x1-x2*x2)*(x1-x3)-(x1*x1-x3*x3)*(x1-x2));
 390   Double_t b = (y1-y2 - a*(x1*x1-x2*x2))/(x1-x2);
 391 //  Double_t c = y1 - a*x1*x1-b*x1;
 392   return -0.5*b/a;
 393 }
 394 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 395 Bool_t AliHMPIDReconHTA::FitFree(Double_t thTrkRec,Double_t phiTrkRec)
 396 {
 397 // Fit performed by minimizing RMS/sqrt(n) of the
 398 // photons reconstructed. First phi is fixed and theta
 399 // is fouond, then (th,ph) of the track
 400 // as free parameters
 401 // Arguments:    PhiRec phi of the track
 402 //   Returns:    none
 403
 404   TMinuit *pMinuit = new TMinuit(2);
 405   pMinuit->mncler();
 406   gMinuit->SetObjectFit((TObject*)this);  gMinuit->SetFCN(AliHMPIDReconHTA::FunMinPhot);  //set fit function
 407   Double_t aArg=-1,parStep,parLow,parHigh;     Int_t iErrFlg;                 //tmp vars for TMinuit
 408   Double_t d1,d2,d3;
 409   TString sName;
 410   Double_t th,ph;
 411
 412   pMinuit->mnexcm("SET PRI",&aArg,1,iErrFlg);                                          //suspend all printout from TMinuit
 413   pMinuit->mnexcm("SET NOW",&aArg,0,iErrFlg);
 414
 415   if(thTrkRec==0) thTrkRec = 3.*TMath::DegToRad();    // not to start from the edge...
 416
 417   pMinuit->mnparm(0," thTrk  ",thTrkRec ,parStep=0.01,parLow=0,parHigh=TMath::PiOver4(),iErrFlg);
 418   pMinuit->mnparm(1," phiTrk ",phiTrkRec,parStep=0.01,parLow=0,parHigh=TMath::TwoPi(),iErrFlg);
 419
 420   pMinuit->FixParameter(1);
 421   pMinuit->mnexcm("SIMPLEX" ,&aArg,0,iErrFlg);
 422   pMinuit->mnexcm("MIGRAD"  ,&aArg,0,iErrFlg);
 423   pMinuit->Release(1);
 424   pMinuit->mnexcm("MIGRAD"  ,&aArg,0,iErrFlg);
 425
 426   pMinuit->mnpout(0,sName,th,d1,d2,d3,iErrFlg);
 427   pMinuit->mnpout(1,sName,ph,d1,d2,d3,iErrFlg);
 428
 429   Double_t f,par[2];
 430   Double_t *grad=0x0;
 431   par[0] = th;par[1] = ph;
 432   pMinuit->Eval(2,grad,f,par,3);
 433
 434   SetTrkFit(th,ph);
 435   return kTRUE;
 436 }
 437 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 438 void AliHMPIDReconHTA::FunMinPhot(Int_t &/* */,Double_t* /* */,Double_t &f,Double_t *par,Int_t iflag)
 439 {
 440 // Minimization function to find best track and thetaC parameters
 441 // Arguments:    f = function value to minimize
 442 //             par = list of parameter to find
 443 //           iflag = flag status. See Minuit instructions
 444 //   Returns:    none
 445 //
 446 // Note: it is necessary to call an instance of AlihMPIDParam. Not possible to use fParam
 447 // because of the static instantiation of the function in Minuit
 448
 449   AliHMPIDParam *pParam=AliHMPIDParam::Instance();
 450   AliHMPIDReconHTA *pRecHTA=(AliHMPIDReconHTA*)gMinuit->GetObjectFit();
 451   AliHMPIDRecon pRec;
 452   Double_t sizeCh = 0.5*pParam->RadThick()+pParam->WinThick()+pParam->GapThick();
 453   Double_t thTrk = par[0];
 454   Double_t phTrk = par[1];
 455   Double_t xrad = pRecHTA->MipX() - sizeCh*TMath::Tan(thTrk)*TMath::Cos(phTrk);
 456   Double_t yrad = pRecHTA->MipY() - sizeCh*TMath::Tan(thTrk)*TMath::Sin(phTrk);
 457   pRecHTA->SetRadXY(xrad,yrad);
 458   pRec.SetTrack(xrad,yrad,thTrk,phTrk);
 459
 460   Double_t meanCkov =0;
 461   Double_t meanCkov2=0;
 462   Double_t thetaCer,phiCer;
 463   Int_t nClAcc = 0;
 464   Int_t nClTot=pRecHTA->NClu();
 465
 466   for(Int_t i=0;i<nClTot;i++) {
 467     if(!(pRecHTA->ClCk(i))) continue;
 468     pRec.FindPhotCkov(pRecHTA->XClu(i),pRecHTA->YClu(i),thetaCer,phiCer);
 469     meanCkov  += thetaCer;
 470     meanCkov2 += thetaCer*thetaCer;
 471     nClAcc++;
 472   }
 473   if(nClAcc==0) {f=999;return;}
 474   meanCkov  /=(Double_t)nClAcc;
 475   meanCkov2 /=(Double_t)nClAcc;
 476   Double_t rms = TMath::Sqrt(TMath::Abs(meanCkov2 - meanCkov*meanCkov));
 477   f = rms/TMath::Sqrt((Double_t)nClAcc);
 478
 479   if(iflag==3) {
 480 /*
 481     Printf("FunMinPhot before: photons candidates %i used %i",nClTot,nClAcc);
 482     nClAcc = 0;
 483     Double_t meanCkov1=0;
 484     Double_t meanCkov2=0;
 485     for(Int_t i=0;i<nClTot;i++) {
 486       if(!(pRec->ClCk(i))) continue;
 487       pRec->FindPhotCkov(pRec->XClu(i),pRec->YClu(i),thetaCer,phiCer);
 488       if(TMath::Abs(thetaCer-meanCkov)<2*rms) {
 489         meanCkov1 += thetaCer;
 490         meanCkov2 += thetaCer*thetaCer;
 491         nClAcc++;
 492       } else pRec->SetClCk(i,kFALSE);
 493     }
 494     meanCkov1/=nClAcc;
 495     Double_t rms2 = (meanCkov2 - meanCkov*meanCkov*nClAcc)/nClAcc;
 496     Printf("FunMinPhot after: photons candidates %i used %i thetaCer %f",nClTot,nClAcc,meanCkov1);
 497     pRec->SetCkovFit(meanCkov1);
 498     pRec->SetCkovSig2(rms2);
 499     pRec->SetNClu(nClAcc);
 500 */
 501 //    Printf("FunMinPhot: reconstructed theta Cerenkov %f with %d photons",meanCkov,nClAcc);
 502     pRecHTA->SetCkovFit(meanCkov);
 503     pRecHTA->SetCkovSig2(rms*rms);
 504     pRecHTA->SetNClu(nClAcc);
 505   }
 506 }//FunMinPhot()
 507 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 508 void AliHMPIDReconHTA::InitDatabase()
 509 {
 510 // Construction a database of ring shapes on fly
 511 //   Arguments: none
 512 //   Returns  : none
 513 //  N.B. fgDatabase points to a TH2I with x-min dist from MIP
 514 //                                        y-dist from the ring of the MIP perpendicular to major axis
 515 //        The content is the packed info of track theta and thetaC in degrees
 516 //                        thetaC+1000*thTrk
 517 //
 518 //  TFile *pout = new TFile("./database.root","recreate");
 519
 520   TStopwatch timer;
 521   timer.Start();
 522
 523
 524 //  if(!fgDatabase) fgDatabase = new TH2F("deconv","database;d1;d2;thC+1000*thTrk",500,0,50,150,0,15);
 525   if(fgDatabase->GetEntries()>=1) {
 526    AliInfo("HTA database already built. ");
 527    return;
 528   }
 529   AliInfo(Form("database HTA is being built.Please, wait..."));
 530 //
 531   Double_t x[3]={0,0,0},y[3];
 532
 533   AliHMPIDRecon rec;
 534
 535   if(!fParam) fParam=AliHMPIDParam::Instance();
 536   Double_t thetaMax = TMath::ACos(1./fParam->MeanIdxRad());
 537   Double_t thTrkMax = 1./TMath::ASin(fParam->MeanIdxRad());
 538
 539   Int_t nstepx = 1000;
 540   Int_t nstepy = 1000;
 541
 542 //  TH2F *fgDatabase = new TH2F("deconv","database;d1;d2;thC+1000*thTrk",500,0,50,150,0,15);
 543   //
 544   Double_t xrad = 0;
 545   Double_t yrad = 0;
 546   Double_t phTrk = 0;
 547
 548   for(Int_t i=0;i<nstepx;i++) {     //loop on thetaC
 549     for(Int_t j=0;j<nstepy;j++) {   //loop on theta particle
 550       Double_t thetaC = thetaMax/nstepx*((Double_t)i+0.5);
 551       Double_t thTrk  = thTrkMax/nstepy*((Double_t)j+0.5);
 552       //
 553       //mip position
 554       //
 555       Double_t sizeCh = 0.5*fParam->RadThick()+fParam->WinThick()+fParam->GapThick();
 556       Double_t xmip = xrad + sizeCh*TMath::Tan(thTrk)*TMath::Cos(phTrk);
 557       Double_t ymip = yrad + sizeCh*TMath::Tan(thTrk)*TMath::Sin(phTrk);
 558
 559       Double_t dist1,dist2;
 560       //
 561       //first point at phi=0
 562       //
 563       rec.SetTrack(xrad,yrad,thTrk,phTrk);
 564       TVector2 pos;
 565       pos=rec.TracePhot(thetaC,0);
 566
 567       if(pos.X()==-999) {
 568         dist1 = 0;            //open ring...only the distance btw mip and point at 180 will be considered
 569       } else {
 570         x[0] = pos.X(); y[0] = pos.Y();
 571         dist1   = TMath::Sqrt((x[0]-xmip)*(x[0]-xmip)+(y[0]-ymip)*(y[0]-ymip));
 572       }
 573       //
 574       //second point at phi=180
 575       //
 576       rec.SetTrack(xrad,yrad,thTrk,phTrk);
 577       pos=rec.TracePhot(thetaC,TMath::Pi());
 578
 579       if(pos.X()==-999) {Printf("it should not happens!Bye");return;}
 580       x[1] = pos.X(); y[1] = pos.Y();
 581       if((x[1]-xmip)*(x[0]-xmip)>0) continue; // to avoid circles out mips (for very low ThetaC)
 582       dist2   = TMath::Sqrt((x[1]-xmip)*(x[1]-xmip)+(y[1]-ymip)*(y[1]-ymip));
 583
 584 //      Double_t distA = dist1+dist2;
 585       Double_t distA = dist2;     // only the minimum: problem of acceptance
 586       //
 587       //second point at phi=90
 588       //
 589       rec.SetTrack(xrad,yrad,thTrk,phTrk);
 590       pos=rec.TracePhot(thetaC,TMath::PiOver2());
 591
 592       if(pos.X()==-999) continue;
 593       x[2] = pos.X(); y[2] = pos.Y();
 594       Double_t distB   = TMath::Sqrt((x[2]-xmip)*(x[2]-xmip)+(y[2]-ymip)*(y[2]-ymip));
 595 // compact the infos...
 596       Int_t compact = (Int_t)(thetaC*TMath::RadToDeg())+1000*(Int_t)(thTrk*TMath::RadToDeg());
 597       Int_t bin = fgDatabase->FindBin(distA,distB);
 598       if(fgDatabase->GetBinContent(bin)==0) fgDatabase->Fill(distA,distB,compact);
 599       Int_t binxDB,binyDB;
 600       FindBinDB(distA,distB,binxDB,binyDB);
 601       fgDB[binxDB][binyDB] = compact;
 602     }
 603   }
 604
 605   FillZeroChan();
 606 //  fgDatabase = deconv;
 607
 608   timer.Stop();
 609   Double_t nSecs = timer.CpuTime();
 610   AliInfo(Form("database HTA successfully open in %3.1f sec.(CPU). Reconstruction is started.",nSecs));
 611
 612 //  pout->Write();
 613 //  pout->Close();
 614
 615 }//InitDatabase()
 616 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 617 void AliHMPIDReconHTA::FillZeroChan()const
 618 {
 619   //If fills eventually channel without entries
 620   //inthe histo "database" jyst interpolating the neighboring cells
 621   // Arguments: histogram pointer of the database
 622   //   Returns: none
 623   //
 624   Int_t nbinx = fgDatabase->GetNbinsX();
 625   Int_t nbiny = fgDatabase->GetNbinsY();
 626   for(Int_t i = 0;i<nbinx;i++) {
 627     for(Int_t j = 0;j<nbiny;j++) {
 628       if(fgDatabase->GetBinContent(i,j) == 0) {
 629         fgDatabase->SetCellContent(i,j,-1);
 630         Int_t nXmin = i-1; Int_t nXmax=i+1;
 631         Int_t nYmin = j-1; Int_t nYmax=j+1;
 632         Int_t nc = 0;
 633         Double_t meanC =0;
 634         Double_t meanTrk =0;
 635         for(Int_t ix=nXmin;ix<=nXmax;ix++) {
 636           if(ix<0||ix>nbinx) continue;
 637           for(Int_t iy=nYmin;iy<=nYmax;iy++) {
 638             if(iy<0||iy>nbiny) continue;
 639             meanC  +=  (Int_t)fgDatabase->GetBinContent(ix,iy)%1000;
 640             meanTrk+=  (Int_t)fgDatabase->GetBinContent(ix,iy)/1000;
 641             nc++;
 642           }
 643           meanC/=nc; meanTrk/=nc;
 644           Int_t compact = (Int_t)meanC+1000*(Int_t)meanTrk;
 645           if(compact>0)fgDatabase->SetCellContent(i,j,compact);
 646         }
 647       }
 648     }
 649   }
 650 }//FillZeroChan()
 651 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 652 Int_t AliHMPIDReconHTA::r2(Double_t *coef, Double_t &x1, Double_t &x2)
 653 {
 654   //2nd deg. equation
 655   //solution
 656   // Arguments: coef 2 1 0: ax^2+bx+c=0
 657   //   Returns: n. of solutions
 658   //            x1= 1st sol
 659   //            x2= 2nd sol
 660   Double_t a,b,c;
 661   a = coef[2];
 662   b = coef[1];
 663   c = coef[0];
 664   Double_t delta = b*b-4*a*c;
 665   if(delta<0) {return 0;}
 666   if(delta==0) {
 667     x1=x2=-b/(2*a);
 668     return 1;
 669   }
 670   // delta>0
 671   x1 = (-b+TMath::Sqrt(delta))/(2*a);
 672   x2 = (-b-TMath::Sqrt(delta))/(2*a);
 673   return 2;
 674 }//r2()
 675 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 676
 677 Double_t AliHMPIDReconHTA::FindSimmPhi()
 678 {
 679 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 680 // RESTITUISCE  IN  OUTPUT  IL  VALORE  FINALE  DELL'ANGOLO  RICOSTRUITO
 681 // CON  I  DUE  METODI:
 682 // - metodo dei minimi quadrati con le distanze effettive;................................(PER RING CHIUSI)
 683 // - metodo della determin della pendenza individuando la distanza minima mip-fotone;.....(PER RING APERTI)
 684
 685   Float_t coeff1ord=0;     Float_t coeff2ord=0;     Float_t coeff0ord=0;
 686   Float_t xrotsumm =0;     Float_t yrotsumm =0;     Float_t xx =0;
 687   Float_t yy =0;           Float_t xy =0;           Float_t yx =0;
 688   Double_t xmin=0;
 689   Double_t ymin=0;
 690
 691   Int_t np=0;
 692
 693   Double_t distMin = 999.;
 694
 695   for(Int_t i=0;i<fNClu;i++) {
 696     if(!fClCk[i]) continue;
 697     np++;
 698     xrotsumm+=fXClu[i];         // summ xi
 699     yrotsumm+=fYClu[i];         // summ yi
 700     xx+=fXClu[i]*fXClu[i];      // summ xixi
 701     yy+=fYClu[i]*fYClu[i];      // summ yiyi
 702     xy+=fXClu[i]*fYClu[i];      // summ yixi
 703     Double_t dist2= (fXClu[i]-fMipX)*(fXClu[i]-fMipX)+(fYClu[i]-fMipY)*(fYClu[i]-fMipY);
 704     if(dist2<distMin) {
 705       distMin = dist2;
 706       xmin = fXClu[i];
 707       ymin = fYClu[i];
 708     }
 709   }
 710
 711   Double_t AngM=0;
 712   if(ymin <  fMipY && xmin >  fMipX)  {AngM  =  TMath::ATan((ymin-fMipY)/(xmin-fMipX))*TMath::RadToDeg()+180;}
 713   if(ymin >  fMipY && xmin <  fMipX)  {AngM  =  TMath::ATan((ymin-fMipY)/(xmin-fMipX))*TMath::RadToDeg()+360;}
 714   if(ymin >  fMipY && xmin >  fMipX)  {AngM  =  TMath::ATan((ymin-fMipY)/(xmin-fMipX))*TMath::RadToDeg()+180;}
 715   if(ymin <  fMipY && xmin <  fMipX)  {AngM  =  TMath::ATan((ymin-fMipY)/(xmin-fMipX))*TMath::RadToDeg();}
 716   if(ymin == fMipY && xmin >  fMipX)  {AngM  =  TMath::ATan((ymin-fMipY)/(xmin-fMipX))*TMath::RadToDeg()+180;}
 717   if(ymin == fMipY && xmin <  fMipX)  {AngM  =  TMath::ATan((ymin-fMipY)/(xmin-fMipX))*TMath::RadToDeg();}
 718   if(ymin <  fMipY && xmin == fMipX)  {AngM  =  90;}
 719   if(ymin >  fMipY && xmin == fMipX)  {AngM  =  270;}
 720
 721   //_____calc. met min quadr using effective distance _________________________________________________
 722
 723   coeff2ord= xy-xrotsumm*yrotsumm/np;
 724   coeff1ord= yrotsumm*yrotsumm/np - xrotsumm*xrotsumm/np - yy + xx;
 725   coeff0ord= xrotsumm*yrotsumm/np -  yx;
 726   Double_t m1=0, m2=0;  Double_t n1=0, n2=0;
 727                             // c           // b         // a
 728   Double_t coeff[3]={coeff0ord,coeff1ord,coeff2ord};
 729
 730   r2(coeff,m1,m2);
 731
 732   n1=(yrotsumm-m1*xrotsumm)/np;
 733   n2=(yrotsumm-m2*xrotsumm)/np;
 734   // 2 solutions.................
 735
 736   Double_t PhiTrk1= TMath::ATan(m1);
 737   Double_t PhiTrk2= TMath::ATan(m2);
 738
 739   // negative angles solved...
 740
 741   Double_t PhiTrk1Positive=0;
 742
 743   if(PhiTrk1<0)    PhiTrk1Positive= PhiTrk1 + 180*TMath::DegToRad();
 744   if(PhiTrk1>=0)   PhiTrk1Positive= PhiTrk1;
 745
 746   Double_t d1 =(1/(m1*m1+1))*(yy+m1*m1*xx+np*n1*n1-2*m1*xy-2*n1*yrotsumm+2*m1*n1*xrotsumm);
 747   Double_t d2 =(1/(m2*m2+1))*(yy+m2*m2*xx+np*n2*n2-2*m2*xy-2*n2*yrotsumm+2*m2*n2*xrotsumm);
 748
 749   Double_t mMin;
 750   if(d1 > d2) mMin = m2; else mMin = m1;
 751
 752   Double_t PhiTrk = TMath::ATan(mMin)*TMath::RadToDeg();
 753   Double_t PhiTrkPositive=0;
 754   //
 755   if(ymin <  fMipY && xmin >  fMipX)  {PhiTrkPositive  =  TMath::ATan(mMin)*TMath::RadToDeg()+180;}
 756   if(ymin >  fMipY && xmin <  fMipX)  {PhiTrkPositive  =  TMath::ATan(mMin)*TMath::RadToDeg()+360;}
 757   if(ymin >  fMipY && xmin >  fMipX)  {PhiTrkPositive  =  TMath::ATan(mMin)*TMath::RadToDeg()+180;}
 758   if(ymin <  fMipY && xmin <  fMipX)  {PhiTrkPositive  =  TMath::ATan(mMin)*TMath::RadToDeg();}
 759   if(ymin == fMipY && xmin >  fMipX)  {PhiTrkPositive  =  TMath::ATan(mMin)*TMath::RadToDeg()+180;}
 760   if(ymin == fMipY && xmin <  fMipX)  {PhiTrkPositive  =  TMath::ATan(mMin)*TMath::RadToDeg();}
 761   if(ymin <  fMipY && xmin == fMipX)  {PhiTrkPositive  =   90;}
 762   if(ymin >  fMipY && xmin == fMipX)  {PhiTrkPositive  =  270;}
 763
 764   //  ------------------------- choose the best-----------------------
 765
 766
 767   Double_t PhiTrkFinal=0;
 768   if( AngM-40 <=  PhiTrkPositive  &&  AngM+40 >= PhiTrkPositive)   PhiTrkFinal = PhiTrkPositive;  else PhiTrkFinal = AngM;
 769
 770   return PhiTrkFinal*TMath::DegToRad();
 771 }
 772 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 773 void AliHMPIDReconHTA::FindBinDB(Double_t x,Double_t y,Int_t &binX,Int_t &binY)
 774 {
 775   const Int_t nxDB = 500;
 776   const Double_t xlowDB =  0;
 777   const Double_t xhigDB = 50;
 778   const Double_t ylowDB =  0;
 779   const Double_t yhigDB = 15;
 780
 781   binX = -1;
 782   binY = -1;
 783   if(x<xlowDB && x>xlowDB &&
 784      y<ylowDB && y>ylowDB)    return;
 785   binX = (Int_t)((x-xlowDB)/(xhigDB-xlowDB));
 786   binY = (Int_t)((y-ylowDB)/(yhigDB-ylowDB));
 787 }