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