HMPID/AliHMPIDReconHTA.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //////////////////////////////////////////////////////////////////////////
  17 //                                                                      //
  18 // AliHMPIDReconHTA                                                         //
  19 //                                                                      //
  20 // HMPID class to perfom pattern recognition based on Hough transfrom    //
  21 // for single chamber                                                   //
  22 //////////////////////////////////////////////////////////////////////////
  23
  24 #include "AliHMPIDReconHTA.h"//class header
  25 #include "AliHMPIDCluster.h" //CkovHiddenTrk()
  26 #include "AliHMPIDRecon.h"   //FunMinPhot()
  27 #include <TFile.h>           //Database()
  28 #include <TMinuit.h>         //FitFree()
  29 #include <TClonesArray.h>    //CkovHiddenTrk()
  30 #include <AliESDtrack.h>     //CkovHiddenTrk()
  31 #include <TH2F.h>            //InitDatabase()
  32 #include <TGraph.h>          //ShapeModel()
  33 #include <TSpline.h>         //ShapeModel()
  34 #include <TCanvas.h>         //ShapeModel()
  35 #include "TStopwatch.h"      //
  36
  37 Int_t AliHMPIDReconHTA::fgDB[500][150]={{75000*0}};
  38 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  39 AliHMPIDReconHTA::AliHMPIDReconHTA():
  40   TNamed("RichRec","RichPat"),
  41   fMipX(-999),
  42   fMipY(-999),
  43   fMipQ(-999),
  44   fRadX(-999),
  45   fRadY(-999),
  46   fIdxMip(0),
  47   fNClu(0),
  48   fXClu(0),
  49   fYClu(0),
  50   fPhiPhot(0),
  51   fThetaPhot(0),
  52   fClCk(0),
  53   fThTrkIn(-999),
  54   fPhTrkIn(-999),
  55   fThTrkFit(-999),
  56   fPhTrkFit(-999),
  57   fCkovFit(-999),
  58   fNCluFit(0),
  59   fCkovSig2(0),
  60   fFitStatus(0),
  61   fParam(AliHMPIDParam::Instance())
  62 {
  63 //..
  64 //hidden algorithm
  65 //..
  66   fParam->SetRefIdx(fParam->MeanIdxRad()); // initialization of ref index to a default one
  67   InitDatabase();
  68 }
  69 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  70 AliHMPIDReconHTA::~AliHMPIDReconHTA()
  71 {
  72   DeleteVars();
  73 }
  74 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  75 void AliHMPIDReconHTA::InitVars(Int_t n)
  76 {
  77 //..
  78 //Init some variables
  79 //..
  80   fXClu = new Double_t[n];
  81   fYClu = new Double_t[n];
  82   fPhiPhot = new Double_t[n];
  83   fThetaPhot = new Double_t[n];
  84   fClCk = new Bool_t[n];
  85   for(Int_t i=0;i<n;i++) fClCk[i] = kTRUE;
  86 //
  87 }
  88 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  89 void AliHMPIDReconHTA::DeleteVars()const
  90 {
  91 //..
  92 //Delete variables
  93 //..
  94   if(fXClu) delete [] fXClu;
  95   if(fYClu) delete [] fYClu;
  96   if(fPhiPhot) delete [] fPhiPhot;
  97   if(fThetaPhot) delete [] fThetaPhot;
  98   if(fClCk) delete [] fClCk;
  99 }
 100 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 101 Bool_t AliHMPIDReconHTA::CkovHiddenTrk(AliESDtrack *pTrk,TClonesArray *pCluLst,Int_t index, Double_t nmean)
 102 {
 103 // Pattern recognition method without any infos from tracking:HTA (Hidden Track Algorithm)...
 104 // The method finds in the chmuber the cluster with the highest charge
 105 // compatibile with a MIP, then the strategy is applied
 106 // Arguments:  pTrk     - pointer to ESD track
 107 //             pCluLs   - list of clusters for a given chamber
 108 //             pNmean   - pointer to ref. index
 109 //             pQthre   - pointer to qthre
 110 //   Returns:           - 0=ok,1=not fitted
 111
 112   AliHMPIDParam *pParam = AliHMPIDParam::Instance();
 113
 114   if(!CluPreFilter(pCluLst)) return kFALSE;
 115
 116   Int_t nCh=0;
 117   Int_t sizeClu=0;
 118
 119   fNClu = pCluLst->GetEntriesFast();
 120
 121   for (Int_t iClu=0;iClu<fNClu;iClu++){                                                       //clusters loop
 122     AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(iClu);                       //get pointer to current cluster
 123     fXClu[iClu] = pClu->X();fYClu[iClu] = pClu->Y();                                          //store x,y for fitting procedure
 124     fClCk[iClu] = kTRUE;                                                                      //all cluster are accepted at this stage to be reconstructed
 125
 126     if(iClu == index) {
 127
 128       fMipX = pClu->X();
 129       fMipY = pClu->Y();
 130       fMipQ = pClu->Q();
 131       sizeClu = pClu->Size();
 132       nCh = pClu->Ch();
 133       fClCk[index] = kFALSE;
 134       fIdxMip = index;
 135       AliDebug(1,Form(" MIP n. %i x %f y %f Q %f",iClu,pClu->X(),pClu->Y(),pClu->Q()));
 136     }
 137   }//clusters loop
 138
 139   pParam->SetRefIdx(nmean);
 140
 141   //
 142   Float_t xra,yra,th,ph; pTrk->GetHMPIDtrk(xra,yra,th,ph);
 143   AliDebug(1,Form(" simulated phiTRK %6.2f thetaTRK %6.2f",ph*TMath::RadToDeg(),th*TMath::RadToDeg()));
 144   //
 145
 146   if(!DoRecHiddenTrk()) {
 147     pTrk->SetHMPIDsignal(pParam->kNoPhotAccept);
 148     return kFALSE;
 149   }                                                                           //Do track and ring reconstruction,if problems returns 1
 150   AliDebug(1,Form("    fitted phi %6.2f ",fPhTrkFit*TMath::RadToDeg()));
 151
 152   pTrk->SetHMPIDtrk(fRadX,fRadY,fThTrkFit,fPhTrkFit);                                        //store track intersection info
 153   pTrk->SetHMPIDmip(fMipX,fMipY,(Int_t)fMipQ,NCluFit());                                     //store mip info + n. phots of the ring
 154   pTrk->SetHMPIDcluIdx(nCh,fIdxMip+1000*sizeClu);                                            //set cham number, index of cluster + cluster size
 155   pTrk->SetHMPIDsignal(fCkovFit);                                                            //find best Theta ckov for ring i.e. track
 156   pTrk->SetHMPIDchi2(fCkovSig2);                                                             //errors squared
 157   AliDebug(1,Form(" n clusters tot %i fitted to ring %i",fNClu,NCluFit()));
 158   for(Int_t i=0;i<fNClu;i++) {
 159     AliDebug(1,Form(" n.%3i  ThetaCer %8.3f PhiCer %8.3f check %i",i,fThetaPhot[i],fPhiPhot[i],fClCk[i]));
 160   }
 161   AliDebug(1,Form("CkovHiddenTrk: thetaC %f th %f ph %f",fCkovFit,fThTrkFit*TMath::RadToDeg(),fPhTrkFit*TMath::RadToDeg()));
 162
 163   return kTRUE;
 164
 165 }//CkovHiddenTrk()
 166 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 167 Bool_t AliHMPIDReconHTA::DoRecHiddenTrk()
 168 {
 169 // Pattern recognition method without any infos from tracking...
 170 // First a preclustering filter to avoid part of the noise
 171 // Then only ellipsed-rings are fitted (no possibility,
 172 // for the moment, to reconstruct very inclined tracks)
 173 // Finally a fitting with (th,ph) free, starting by very close values
 174 // previously evaluated.
 175 // Arguments:   none
 176 //   Returns:   none
 177   Double_t thTrkRec,phiTrkRec,thetaCRec;
 178
 179   if(!FindShape(thTrkRec,phiTrkRec,thetaCRec)) {
 180     AliDebug(1,Form(" FindShape failed...!"));
 181     return kFALSE;
 182   }
 183   AliDebug(1,Form(" FindShape accepted...!"));
 184
 185   if(!FitRing(thTrkRec,phiTrkRec)) {
 186     AliDebug(1,Form(" FitRing failed...!"));
 187     return kFALSE;
 188   }
 189   AliDebug(1,Form(" FitRing accepted...!"));
 190
 191   if(!UniformDistrib()) {
 192     AliDebug(1,Form(" UniformDistrib failed...!"));
 193     return kFALSE;
 194   }
 195   AliDebug(1,Form(" UniformDistrib accepted...!"));
 196
 197   return kTRUE;
 198 }//DoRecHiddenTrk()
 199 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 200 Bool_t AliHMPIDReconHTA::CluPreFilter(TClonesArray *pCluLst)
 201 {
 202 // Pre-filter of bkg clusters
 203 // Arguments:    pSluLst  -  List of the clusters for a given chamber
 204 //   Returns:    status   -  TRUE if filtering leaves enough photons, FALSE if not
 205 //
 206   Int_t nClusTot = pCluLst->GetEntriesFast();
 207   if(nClusTot<4||nClusTot>100) {
 208     return kFALSE;
 209   } else {
 210     InitVars(nClusTot);
 211     return kTRUE;
 212   }
 213 }
 214 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 215 Bool_t AliHMPIDReconHTA::FindShape(Double_t &thTrkRec,Double_t &phiTrkRec,Double_t &thetaCRec)
 216 {
 217 // Finds the estimates for phi and theta of the track and the ThetaCerenkov
 218 // by using a database of the shapes of the rings
 219 // Arguments:   none
 220 //   Returns:   thTrkRec  - estimate of theta track
 221 //              phiTrkRec - estimate of phi   track
 222 //              thetaCRec - estimate of ThetaCerenkov
 223 //              status    - TRUE if a good solution is found, FALSE if not
 224
 225   Double_t phiphot[1000];
 226   Double_t    dist[1000];
 227   Int_t     indphi[1000];
 228
 229   Bool_t status;
 230
 231   if(fNClu>1000) return kFALSE;  // too many clusters....
 232
 233 // Sort in phi angle...
 234   for(Int_t i=0;i<fNClu;i++) {
 235     if(!fClCk[i]) {
 236       AliDebug(1,Form(" n.%3i  xMIP    %8.3f yMIP %8.3f check %i",i,fMipX,fMipY,fClCk[i]));
 237       phiphot[i] = 999.;
 238       dist[i]    = 999.;
 239       continue;
 240     }
 241     phiphot[i] = (TMath::ATan2(fMipY-fYClu[i],fMipX-fXClu[i])+TMath::Pi())*TMath::RadToDeg();
 242     dist[i]=TMath::Sqrt((fMipX-fXClu[i])*(fMipX-fXClu[i])+(fMipY-fYClu[i])*(fMipY-fYClu[i]));
 243     AliDebug(1,Form(" n.%3i  phiphot %8.3f dist %8.3f check %i",i,phiphot[i],dist[i],fClCk[i]));
 244   }
 245
 246   TMath::Sort(fNClu,phiphot,indphi,kFALSE);
 247
 248 // Purify with a truncated mean;
 249   Int_t np=0;
 250   Double_t dMean  = 0;
 251   Double_t dMean2 = 0;
 252   for(Int_t i=0;i<fNClu;i++) {
 253     if(!fClCk[indphi[i]]) continue;                                                  // Check if a good photon candidate or not
 254     dMean +=dist[indphi[i]];
 255     dMean2+=dist[indphi[i]]*dist[indphi[i]];
 256     np++;
 257   }
 258
 259   if(np>0){
 260    dMean  /=(Double_t)np;
 261    dMean2 /=(Double_t)np;}
 262   Double_t rms = TMath::Sqrt(dMean2 - dMean*dMean);
 263
 264   for(Int_t i=0;i<fNClu;i++) {
 265     if(!fClCk[indphi[i]]) continue;                                                  // Check if a good photon candidate or not
 266     if(TMath::Abs(dMean-dist[indphi[i]]) > 1.5*rms) {
 267       fClCk[indphi[i]] = kFALSE;
 268       continue;
 269     }
 270   }
 271
 272   AliDebug(1,"Purification of photons...");
 273 //
 274 //  purify vectors for good photon candidates
 275 //
 276   Int_t npeff=0;
 277   Double_t *phiphotP = new Double_t[fNClu+1];
 278   Double_t *distP    = new Double_t[fNClu+1];
 279   for(Int_t i=0;i<fNClu;i++) {
 280     AliDebug(1,Form(" n. %3i phiphot %8.3f dist %8.3f check %i",i,phiphot[indphi[i]],dist[indphi[i]],fClCk[indphi[i]]));
 281     if(!fClCk[indphi[i]]) continue;                                                  // Check if a good photon candidate or not
 282     phiphotP[npeff] = phiphot[indphi[i]];
 283     distP[npeff]    = dist[indphi[i]];
 284     npeff++;
 285   }
 286
 287   if(npeff<3) {
 288     AliDebug(1,Form("FindShape failed: no enough photons = %i...",npeff));
 289     delete [] phiphotP;
 290     delete [] distP;
 291     return kFALSE;
 292   }
 293
 294   Double_t xA,xB;
 295   status = kFALSE;
 296
 297   if(!ShapeModel(npeff,phiphotP,distP,xA,xB,phiTrkRec)) {AliDebug(1,Form("ShapeModel failed            ")); return kFALSE;}
 298
 299 //  if(xA > 50 || xB > 15)                                {AliDebug(1,Form("xA and xB failed out of range")); return kFALSE;}
 300
 301   Int_t binxDB,binyDB;
 302   Int_t compactDB=-1;
 303
 304   if(xA > xB)  {                                        //geometrically not possible, try to recover on a mean values...
 305
 306     FindBinDB(xA,xA,binxDB,binyDB);
 307     if(binxDB<0 || binyDB<0)                              {AliDebug(1,Form("bin < 0 ! failed             ")); return kFALSE;}
 308     Int_t compactDB1 = CompactDB(binxDB,binyDB);
 309     FindBinDB(xB,xB,binxDB,binyDB);
 310     if(binxDB<0 || binyDB<0)                              {AliDebug(1,Form("bin < 0 ! failed             ")); return kFALSE;}
 311     Int_t compactDB2 = CompactDB(binxDB,binyDB);
 312     Double_t thetaCRec1 =  (Double_t)(compactDB1%1000);
 313     Double_t thetaCRec2 =  (Double_t)(compactDB2%1000);
 314     Double_t thTrkRec1  =  (Double_t)(compactDB1/1000);
 315     Double_t thTrkRec2  =  (Double_t)(compactDB2/1000);
 316     thetaCRec = 0.5*(thetaCRec1+thetaCRec2);
 317     thTrkRec  = 0.5*( thTrkRec1+ thTrkRec2);
 318
 319   } else {
 320
 321     FindBinDB(xA,xB,binxDB,binyDB);
 322     if(binxDB<0 || binyDB<0) {AliDebug(1,Form("bin < 0 ! failed             ")); return kFALSE;}
 323
 324     compactDB = CompactDB(binxDB,binyDB);
 325
 326     if(compactDB<0)                                       {AliDebug(1,Form("compact< 0! failed           ")); return kFALSE;}
 327     //
 328     //
 329     thetaCRec =  (Double_t)(compactDB%1000);
 330     thTrkRec  =  (Double_t)(compactDB/1000);
 331
 332   }
 333
 334   AliDebug(1,Form(" CompactDB %i thTrkRec %8.3f thetaCRec %8.3f ",compactDB,thTrkRec,thetaCRec));
 335
 336   phiTrkRec *= TMath::DegToRad();
 337   thTrkRec  *= TMath::DegToRad();
 338   thetaCRec *= TMath::DegToRad();
 339
 340   status = kTRUE;
 341
 342   delete [] phiphotP;
 343   delete [] distP;
 344
 345   return status;
 346 }
 347 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 348 Bool_t AliHMPIDReconHTA::ShapeModel(Int_t np,Double_t *phiphot,Double_t *dist,Double_t &xA,Double_t &xB,Double_t &phiStart)
 349 {
 350 // Find a Spline curve to define dist. vs. phi angle
 351 // in order to estimate the phi of the track
 352 // Arguments:   np     - # points corresponding to # photon candidates
 353 //             dist    - distance of each photon from MIP
 354 //             phiphot - phi of the photon in the DRS
 355 //   Returns:  xA      - min. distance from MIP
 356 //             xB      - dist. from mip perpedicular to the major axis
 357 //             phiStart- estimate of the track phi
 358
 359   TGraph *phigr = new TGraph(np,phiphot,dist);
 360   phiStart = FindSimmPhi();
 361
 362   Double_t phiStart1 = phiStart;
 363   if(phiStart1 > 360) phiStart1 -= 360;
 364   Double_t phiStart2 = phiStart+90;
 365   if(phiStart2 > 360) phiStart2 -= 360;
 366   xA = phigr->Eval(phiStart1);
 367   xB = phigr->Eval(phiStart2);
 368   //---
 369   AliDebug(1,Form("phiStart %f phiStart1 %f phiStart2 %f ",phiStart,phiStart1,phiStart2));
 370   AliDebug(1,Form("xA  %f xB  %f",xA,xB));
 371
 372   phiStart += 180;
 373   if(phiStart>360) phiStart-=360;
 374
 375   return kTRUE;
 376 }
 377 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 378 Double_t AliHMPIDReconHTA::VertParab(Double_t x1,Double_t y1,Double_t x2, Double_t y2, Double_t x3, Double_t y3)const
 379 {
 380 // It uses parabola from 3 points to evaluate the x-coord of the parab
 381 // Arguments:    xi,yi - points
 382 //   Returns:    x-coord of the vertex
 383
 384   Double_t a = ((x1-x3)*(y1-y2)-(x1-x2)*(y1-y3))/((x1*x1-x2*x2)*(x1-x3)-(x1*x1-x3*x3)*(x1-x2));
 385   Double_t b = (y1-y2 - a*(x1*x1-x2*x2))/(x1-x2);
 386 //  Double_t c = y1 - a*x1*x1-b*x1;
 387   return -0.5*b/a;
 388 }
 389 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 390 Bool_t AliHMPIDReconHTA::FitRing(Double_t thTrkRec,Double_t phiTrkRec)
 391 {
 392   Double_t th = thTrkRec;
 393   Double_t ph = phiTrkRec;
 394
 395   FitFree(th,ph);
 396   while(FitStatus()) {
 397     th = fThTrkFit;
 398     ph = fPhTrkFit;
 399     FitFree(th,ph);
 400   }
 401   return kTRUE;
 402 }
 403 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 404 Bool_t AliHMPIDReconHTA::FitFree(Double_t thTrkRec,Double_t phiTrkRec)
 405 {
 406 // Fit performed by minimizing RMS/sqrt(n) of the
 407 // photons reconstructed. First phi is fixed and theta
 408 // is fouond, then (th,ph) of the track
 409 // as free parameters
 410 // Arguments:    PhiRec phi of the track
 411 //   Returns:    none
 412
 413   TMinuit *pMinuit = new TMinuit(2);
 414   pMinuit->mncler();
 415   gMinuit->SetObjectFit((TObject*)this);  gMinuit->SetFCN(AliHMPIDReconHTA::FunMinPhot);  //set fit function
 416   Double_t aArg=-1,parStep,parLow,parHigh;     Int_t iErrFlg;                 //tmp vars for TMinuit
 417   Double_t d1,d2,d3;
 418   TString sName;
 419   Double_t th,ph;
 420
 421   pMinuit->mnexcm("SET PRI",&aArg,1,iErrFlg);                                          //suspend all printout from TMinuit
 422   pMinuit->mnexcm("SET NOW",&aArg,0,iErrFlg);
 423
 424   if(thTrkRec==0) thTrkRec = 3.*TMath::DegToRad();    // not to start from the edge...
 425
 426   AliDebug(1,Form(" Minimization STARTED with phiTRK %6.2f thetaTRK %6.2f",phiTrkRec,thTrkRec));
 427
 428   pMinuit->mnparm(0," thTrk  ",thTrkRec ,parStep=0.01,parLow=0,parHigh=TMath::PiOver4(),iErrFlg);
 429   pMinuit->mnparm(1," phiTrk ",phiTrkRec,parStep=0.01,parLow=0,parHigh=TMath::TwoPi(),iErrFlg);
 430
 431   pMinuit->FixParameter(1);
 432   pMinuit->mnexcm("SIMPLEX" ,&aArg,0,iErrFlg);
 433   pMinuit->mnexcm("MIGRAD"  ,&aArg,0,iErrFlg);
 434   pMinuit->Release(1);
 435   pMinuit->mnexcm("MIGRAD"  ,&aArg,0,iErrFlg);
 436
 437   pMinuit->mnpout(0,sName,th,d1,d2,d3,iErrFlg);
 438   pMinuit->mnpout(1,sName,ph,d1,d2,d3,iErrFlg);
 439
 440   Double_t f,par[2];
 441   Double_t *grad=0x0;
 442   par[0] = th;par[1] = ph;
 443   pMinuit->Eval(2,grad,f,par,3);
 444
 445   SetTrkFit(th,ph);
 446   return kTRUE;
 447 }
 448 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 449 void AliHMPIDReconHTA::FunMinPhot(Int_t &/* */,Double_t* /* */,Double_t &f,Double_t *par,Int_t iflag)
 450 {
 451 // Minimization function to find best track and thetaC parameters
 452 // Arguments:    f = function value to minimize
 453 //             par = list of parameter to find
 454 //           iflag = flag status. See Minuit instructions
 455 //   Returns:    none
 456 //
 457 // Note: it is necessary to call an instance of AlihMPIDParam. Not possible to use fParam
 458 // because of the static instantiation of the function in Minuit
 459
 460   AliHMPIDParam *pParam=AliHMPIDParam::Instance();
 461   AliHMPIDReconHTA *pRecHTA=(AliHMPIDReconHTA*)gMinuit->GetObjectFit();
 462   AliHMPIDRecon pRec;
 463   Double_t sizeCh = 0.5*pParam->RadThick()+pParam->WinThick()+pParam->GapThick();
 464   Double_t thTrk = par[0];
 465   Double_t phTrk = par[1];
 466   Double_t xrad = pRecHTA->MipX() - sizeCh*TMath::Tan(thTrk)*TMath::Cos(phTrk);
 467   Double_t yrad = pRecHTA->MipY() - sizeCh*TMath::Tan(thTrk)*TMath::Sin(phTrk);
 468   pRecHTA->SetRadXY(xrad,yrad);
 469   pRec.SetTrack(xrad,yrad,thTrk,phTrk);
 470
 471   Double_t meanCkov =0;
 472   Double_t meanCkov2=0;
 473   Double_t thetaCer,phiCer;
 474   Int_t nClAcc = 0;
 475   Int_t nClTot=pRecHTA->NClu();
 476
 477   for(Int_t i=0;i<nClTot;i++) {
 478
 479     if(pRecHTA->IdxMip() == i) {
 480       pRecHTA->SetPhotAngles(i,999.,999.);
 481       continue;
 482     }
 483
 484     if(!(pRecHTA->ClCk(i))) continue;
 485
 486     Bool_t status = pRec.FindPhotCkov(pRecHTA->XClu(i),pRecHTA->YClu(i),thetaCer,phiCer);
 487     if(!status) {
 488       pRecHTA->SetPhotAngles(i,999.,999.);
 489       continue;
 490     }
 491     pRecHTA->SetPhotAngles(i,thetaCer,phiCer);
 492     meanCkov  += thetaCer;
 493     meanCkov2 += thetaCer*thetaCer;
 494     nClAcc++;
 495   }
 496
 497   if(nClAcc==0) {f=999;return;}
 498   meanCkov  /=(Double_t)nClAcc;
 499   meanCkov2 /=(Double_t)nClAcc;
 500   Double_t rms = TMath::Sqrt(TMath::Abs(meanCkov2 - meanCkov*meanCkov));
 501   f = rms/TMath::Sqrt((Double_t)nClAcc);
 502
 503   if(iflag==3) {
 504     Int_t nClAccStep1 = nClAcc;
 505     nClAcc = 0;
 506     Double_t meanCkov1=0;
 507     Double_t meanCkov3=0;
 508     for(Int_t i=0;i<nClTot;i++) {
 509
 510       if(!(pRecHTA->ClCk(i))) continue;
 511       thetaCer = pRecHTA->PhotTheta(i);
 512       if(TMath::Abs(thetaCer-meanCkov)<2*rms) {
 513         meanCkov1 += thetaCer;
 514         meanCkov3 += thetaCer*thetaCer;
 515         nClAcc++;
 516       } else pRecHTA->SetClCk(i,kFALSE);
 517     }
 518
 519     if(nClAcc<3) {
 520       pRecHTA->SetFitStatus(kFALSE);
 521       pRecHTA->SetCkovFit(meanCkov);
 522       pRecHTA->SetCkovSig2(rms*rms);
 523       pRecHTA->SetNCluFit(nClAccStep1);
 524       return;
 525     }
 526
 527     meanCkov1/=nClAcc;
 528     Double_t rms2 = (meanCkov3 - meanCkov*meanCkov*nClAcc)/nClAcc;
 529
 530     // get a logger instance
 531     // what for??
 532     //AliLog::GetRootLogger();
 533
 534     if(nClAcc!=nClAccStep1) pRecHTA->SetFitStatus(kTRUE); else pRecHTA->SetFitStatus(kFALSE);
 535
 536     pRecHTA->SetCkovFit(meanCkov1);
 537     pRecHTA->SetCkovSig2(rms2);
 538     pRecHTA->SetNCluFit(nClAcc);
 539   }
 540 }//FunMinPhot()
 541 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 542 void AliHMPIDReconHTA::InitDatabase()
 543 {
 544 // Construction a database of ring shapes on fly
 545 //   Arguments: none
 546 //   Returns  : none
 547 //  N.B. fgDB is the distance with x-min from MIP
 548 //                                 y-dist from the ring of the MIP perpendicular to major axis
 549 //        The content is the packed info of track theta and thetaC in degrees
 550 //                        thetaC+1000*thTrk
 551 //
 552 //  TFile *pout = new TFile("./database.root","recreate");
 553
 554   static Bool_t isDone = kFALSE;
 555
 556   TStopwatch timer;
 557
 558   if(isDone) {
 559    return;
 560   }
 561
 562   if(!isDone) {
 563     timer.Start();
 564   }
 565
 566   AliInfo(Form("database HTA is being built.Please, wait..."));
 567 //
 568   Double_t x[3]={0,0,0},y[3];
 569
 570   AliHMPIDRecon rec;
 571
 572   if(!fParam) fParam=AliHMPIDParam::Instance();
 573   Double_t thetaMax = TMath::ACos(1./fParam->MeanIdxRad());
 574   Double_t thTrkMax = 1./TMath::ASin(fParam->MeanIdxRad());
 575
 576   Int_t nstepx = 1000;
 577   Int_t nstepy = 1000;
 578
 579   //
 580   Double_t xrad = 0;
 581   Double_t yrad = 0;
 582   Double_t phTrk = 0;
 583
 584   for(Int_t i=0;i<nstepx;i++) {     //loop on thetaC
 585     for(Int_t j=0;j<nstepy;j++) {   //loop on theta particle
 586       Double_t thetaC = thetaMax/nstepx*((Double_t)i+0.5);
 587       Double_t thTrk  = thTrkMax/nstepy*((Double_t)j+0.5);
 588       //
 589       //mip position
 590       //
 591       Double_t sizeCh = 0.5*fParam->RadThick()+fParam->WinThick()+fParam->GapThick();
 592       Double_t xmip = xrad + sizeCh*TMath::Tan(thTrk)*TMath::Cos(phTrk);
 593       Double_t ymip = yrad + sizeCh*TMath::Tan(thTrk)*TMath::Sin(phTrk);
 594
 595       Double_t dist1,dist2;
 596       //
 597       //first point at phi=0
 598       //
 599       rec.SetTrack(xrad,yrad,thTrk,phTrk);
 600       TVector2 pos;
 601       pos=rec.TracePhot(thetaC,0);
 602
 603       if(pos.X()==-999) {
 604         dist1 = 0;            //open ring...only the distance btw mip and point at 180 will be considered
 605       } else {
 606         x[0] = pos.X(); y[0] = pos.Y();
 607         dist1   = TMath::Sqrt((x[0]-xmip)*(x[0]-xmip)+(y[0]-ymip)*(y[0]-ymip));
 608       }
 609       //
 610       //second point at phi=180
 611       //
 612       rec.SetTrack(xrad,yrad,thTrk,phTrk);
 613       pos=rec.TracePhot(thetaC,TMath::Pi());
 614
 615       if(pos.X()==-999) {AliDebug(1,Form("it should not happens!Bye"));return;}
 616       x[1] = pos.X(); y[1] = pos.Y();
 617       if((x[1]-xmip)*(x[0]-xmip)>0) continue; // to avoid circles out mips (for very low ThetaC)
 618       dist2   = TMath::Sqrt((x[1]-xmip)*(x[1]-xmip)+(y[1]-ymip)*(y[1]-ymip));
 619
 620 //      Double_t distA = dist1+dist2;
 621       Double_t distA = dist2;     // only the minimum: problem of acceptance
 622       //
 623       //second point at phi=90
 624       //
 625       rec.SetTrack(xrad,yrad,thTrk,phTrk);
 626       pos=rec.TracePhot(thetaC,TMath::PiOver2());
 627
 628       if(pos.X()==-999) continue;
 629       x[2] = pos.X(); y[2] = pos.Y();
 630       Double_t distB   = TMath::Sqrt((x[2]-xmip)*(x[2]-xmip)+(y[2]-ymip)*(y[2]-ymip));
 631 // compact the infos...
 632       Int_t compact = (Int_t)(thetaC*TMath::RadToDeg())+1000*(Int_t)(thTrk*TMath::RadToDeg());
 633       Int_t binxDB,binyDB;
 634       FindBinDB(distA,distB,binxDB,binyDB);
 635       if(fgDB[binxDB][binyDB]==0) fgDB[binxDB][binyDB] = compact;
 636     }
 637   }
 638
 639   FillZeroChan();
 640
 641   if(!isDone) {
 642     timer.Stop();
 643     Double_t nSecs = timer.CpuTime();
 644     AliInfo(Form("database HTA successfully open in %3.1f sec.(CPU). Reconstruction is started.",nSecs));
 645     isDone = kTRUE;
 646   }
 647
 648 //  pout->Write();
 649 //  pout->Close();
 650
 651 }//InitDatabase()
 652 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 653 void AliHMPIDReconHTA::FillZeroChan()const
 654 {
 655   //If fills eventually channel without entries
 656   //inthe histo "database" jyst interpolating the neighboring cells
 657   // Arguments: histogram pointer of the database
 658   //   Returns: none
 659   //
 660   const Int_t nxDB = 500;
 661   const Int_t nyDB = 150;
 662
 663   for(Int_t i = 0;i<nxDB;i++) {
 664     for(Int_t j = 0;j<nyDB;j++) {
 665       if(fgDB[i][j] == 0) {
 666         fgDB[i][j] = -1;
 667         Int_t nXmin = i-1; Int_t nXmax=i+1;
 668         Int_t nYmin = j-1; Int_t nYmax=j+1;
 669         Int_t nc = 0;
 670         Double_t meanC =0;
 671         Double_t meanTrk =0;
 672         for(Int_t ix=nXmin;ix<=nXmax;ix++) {
 673           if(ix<0||ix>=nxDB) continue;
 674           for(Int_t iy=nYmin;iy<=nYmax;iy++) {
 675             if(iy<0||iy>=nyDB) continue;
 676             meanC  +=  (Int_t)(fgDB[ix][iy]%1000);
 677             meanTrk+=  (Int_t)(fgDB[ix][iy]/1000);
 678             nc++;
 679           }
 680           meanC/=nc; meanTrk/=nc;
 681           Int_t compact = (Int_t)meanC+1000*(Int_t)meanTrk;
 682           if(compact>0)fgDB[i][j] = compact;
 683         }
 684       }
 685     }
 686   }
 687 }//FillZeroChan()
 688 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 689 Int_t AliHMPIDReconHTA::r2(Double_t *coef, Double_t &x1, Double_t &x2)
 690 {
 691   //2nd deg. equation
 692   //solution
 693   // Arguments: coef 2 1 0: ax^2+bx+c=0
 694   //   Returns: n. of solutions
 695   //            x1= 1st sol
 696   //            x2= 2nd sol
 697   Double_t a,b,c;
 698   a = coef[2];
 699   b = coef[1];
 700   c = coef[0];
 701   Double_t delta = b*b-4*a*c;
 702   if(delta<0) {return 0;}
 703   if(delta==0) {
 704     x1=x2=-b/(2*a);
 705     return 1;
 706   }
 707   if(a==0) {
 708     x1 = -c/b;
 709     return 1;
 710   }
 711   // delta>0
 712   x1 = (-b+TMath::Sqrt(delta))/(2*a);
 713   x2 = (-b-TMath::Sqrt(delta))/(2*a);
 714   return 2;
 715 }//r2()
 716 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 717
 718 Double_t AliHMPIDReconHTA::FindSimmPhi()
 719 {
 720 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 721 // Reconstruction of phiTRK angle with two methods (in switching)
 722 //
 723 // - least square method (for closed rings)
 724 // - by minimum distance mip-photon (for open rings)
 725
 726   Float_t coeff1ord=0;     Float_t coeff2ord=0;     Float_t coeff0ord=0;
 727   Float_t xrotsumm =0;     Float_t yrotsumm =0;     Float_t xx =0;
 728   Float_t yy =0;           Float_t xy =0;
 729   Double_t xmin=0;
 730   Double_t ymin=0;
 731
 732   Int_t np=0;
 733
 734   Double_t distMin = 999.;
 735
 736   for(Int_t i=0;i<fNClu;i++) {
 737     if(!fClCk[i]) continue;
 738     np++;
 739     xrotsumm+=fXClu[i];         // summ xi
 740     yrotsumm+=fYClu[i];         // summ yi
 741     xx+=fXClu[i]*fXClu[i];      // summ xixi
 742     yy+=fYClu[i]*fYClu[i];      // summ yiyi
 743     xy+=fXClu[i]*fYClu[i];      // summ yixi
 744     Double_t dist2= (fXClu[i]-fMipX)*(fXClu[i]-fMipX)+(fYClu[i]-fMipY)*(fYClu[i]-fMipY);
 745     if(dist2<distMin) {
 746       distMin = dist2;
 747       xmin = fXClu[i];
 748       ymin = fYClu[i];
 749     }
 750   }
 751
 752   Double_t AngM = TMath::ATan2(ymin-fMipY,(xmin-fMipX))*TMath::RadToDeg();
 753   if (AngM<0) AngM+=360;
 754
 755   AliDebug(1,Form(" Simple angle prediction with MIN phi = %f",AngM));
 756
 757   //_____calc. met min quadr using effective distance _________________________________________________
 758
 759   if(np>0){
 760     coeff2ord = xy-xrotsumm*yrotsumm/np;
 761     coeff1ord = yrotsumm*yrotsumm/np - xrotsumm*xrotsumm/np - yy + xx;
 762   }
 763   coeff0ord = -coeff2ord;
 764
 765   AliDebug(1,Form(" a = %f b = %f c = %f",coeff2ord,coeff1ord,coeff0ord));
 766
 767   Double_t m1=0, m2=0;  Double_t n1=0, n2=0;
 768                             // c           // b         // a
 769   Double_t coeff[3]={coeff0ord,coeff1ord,coeff2ord};
 770
 771   r2(coeff,m1,m2);
 772
 773   if(np>0){
 774     n1=(yrotsumm-m1*xrotsumm)/np;
 775     n2=(yrotsumm-m2*xrotsumm)/np;}
 776   // 2 solutions.................
 777
 778   // negative angles solved...
 779
 780   Double_t d1 =(1/(m1*m1+1))*(yy+m1*m1*xx+np*n1*n1-2*m1*xy-2*n1*yrotsumm+2*m1*n1*xrotsumm);
 781   Double_t d2 =(1/(m2*m2+1))*(yy+m2*m2*xx+np*n2*n2-2*m2*xy-2*n2*yrotsumm+2*m2*n2*xrotsumm);
 782
 783   AliDebug(1,Form(" predicted distance d1 %f for angle %6.2f d2 %f for angle %6.2f",d1,TMath::ATan(m1)*TMath::RadToDeg(),
 784                                                                                     d2,TMath::ATan(m2)*TMath::RadToDeg()));
 785   Double_t mMin;
 786   if(d1 > d2) mMin = m2; else mMin = m1;
 787
 788   Double_t phiTrk=0;
 789   //
 790   if(ymin <  fMipY && xmin >  fMipX)  {phiTrk  =  TMath::ATan(mMin)*TMath::RadToDeg()+180;}
 791   if(ymin >  fMipY && xmin <  fMipX)  {phiTrk  =  TMath::ATan(mMin)*TMath::RadToDeg()+360;}
 792   if(ymin >  fMipY && xmin >  fMipX)  {phiTrk  =  TMath::ATan(mMin)*TMath::RadToDeg()+180;}
 793   if(ymin <  fMipY && xmin <  fMipX)  {phiTrk  =  TMath::ATan(mMin)*TMath::RadToDeg();}
 794   if(ymin == fMipY && xmin >  fMipX)  {phiTrk  =  TMath::ATan(mMin)*TMath::RadToDeg()+180;}
 795   if(ymin == fMipY && xmin <  fMipX)  {phiTrk  =  TMath::ATan(mMin)*TMath::RadToDeg();}
 796   if(ymin <  fMipY && xmin == fMipX)  {phiTrk  =   90;}
 797   if(ymin >  fMipY && xmin == fMipX)  {phiTrk  =  270;}
 798
 799   //  ------------------------- choose the best-----------------------
 800
 801
 802   if( AngM-40 <=  phiTrk  &&  AngM+40 >= phiTrk)   return phiTrk;  else return AngM;
 803 }
 804 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 805 void AliHMPIDReconHTA::FindBinDB(Double_t x,Double_t y,Int_t &binX,Int_t &binY)
 806 {
 807   const Int_t nxDB = 500;
 808   const Int_t nyDB = 150;
 809   const Double_t xlowDB =  0;
 810   const Double_t xhigDB = 50;
 811   const Double_t ylowDB =  0;
 812   const Double_t yhigDB = 30;
 813
 814   binX = -1;
 815   binY = -1;
 816   if(x<xlowDB && x>xhigDB &&
 817      y<ylowDB && y>yhigDB)    return;
 818   binX = Int_t((x-xlowDB)/(xhigDB-xlowDB)*nxDB);
 819   binY = Int_t((y-ylowDB)/(yhigDB-ylowDB)*nyDB);
 820   if(binX>=nxDB || binY>=nyDB) {
 821     binX = -1;
 822     binY = -1;
 823   }
 824
 825 }
 826 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 827 Bool_t AliHMPIDReconHTA::UniformDistrib()
 828 {
 829   AliHMPIDParam *pParam=AliHMPIDParam::Instance();
 830   AliHMPIDRecon pRec;
 831
 832   Double_t sizeCh = 0.5*pParam->RadThick()+pParam->WinThick()+pParam->GapThick();
 833   Double_t xrad = MipX() - sizeCh*TMath::Tan(fThTrkFit)*TMath::Cos(fPhTrkFit);
 834   Double_t yrad = MipY() - sizeCh*TMath::Tan(fThTrkFit)*TMath::Sin(fPhTrkFit);
 835
 836   pRec.SetTrack(xrad,yrad,fThTrkFit,fPhTrkFit);
 837
 838   Int_t npeff=0;
 839   Int_t nPhotPerBin = 4;
 840   for(Int_t i=0;i<fNClu;i++) {
 841     if(!ClCk(i)) continue;
 842     npeff++;
 843   }
 844
 845   Int_t nPhiBins = npeff/nPhotPerBin+1;
 846   if(nPhiBins<=1) nPhiBins = 2;
 847
 848   Double_t *iPhiBin;
 849   iPhiBin = new Double_t[nPhiBins];
 850
 851   for(Int_t i=0;i<nPhiBins;i++) iPhiBin[i] =0;
 852
 853   for(Int_t i=0;i<fNClu;i++) {
 854     if(!ClCk(i)) continue;
 855     Double_t phiCer = PhotPhi(i);
 856
 857     Double_t PhiProva = phiCer*TMath::RadToDeg();
 858     if(PhiProva<0) PhiProva+= 360;
 859     Int_t index = (Int_t)((Float_t)nPhiBins*PhiProva/360.);
 860     iPhiBin[index]++;
 861    }
 862
 863    Double_t chi2 = 0;
 864    for(Int_t i=0;i<nPhiBins;i++) {
 865      Double_t theo = (Double_t)npeff/(Double_t)nPhiBins;
 866      if(theo==0) continue;
 867      chi2+= (iPhiBin[i] - theo)*(iPhiBin[i] - theo)/theo;
 868    }
 869
 870     delete [] iPhiBin;
 871
 872     Double_t prob = TMath::Prob(chi2, nPhiBins-1);
 873     AliDebug(1,Form(" Probability for uniform distrib: %6f.3 %s",prob,(prob<0.05) ? "rejected" : "accepted"));
 874     if(prob<0.05) return kFALSE;
 875     return kTRUE;
 876
 877  }
 878 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++