PHOS/AliPHOSEvalRecPoint.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 /* $Id$ */
  17
  18 //*-- Author: Boris Polichtchouk, IHEP
  19 //////////////////////////////////////////////////////////////////////////////
  20 // Reconstructed point operations for the Clusterization class for IHEP reconstruction.
  21 // Performs clusterization (collects neighbouring active cells)
  22 // It differs from AliPHOSClusterizerv1 in neighbour definition only
  23
  24 // ---- ROOT system ---
  25 #include <TMath.h>
  26 #include <TDirectory.h>
  27 #include <TBranch.h>
  28 #include <TTree.h>
  29 #include <TROOT.h>
  30 #include <TFolder.h>
  31
  32 // --- AliRoot header files ---
  33 #include "AliLog.h"
  34 #include "AliConfig.h"
  35 #include "AliPHOSDigit.h"
  36 #include "AliPHOSClusterizer.h"
  37 #include "AliPHOSEvalRecPoint.h"
  38 #include "AliRun.h"
  39 #include "AliPHOSLoader.h"
  40 #include "AliPHOSRecCpvManager.h"
  41 #include "AliPHOSRecEmcManager.h"
  42 #include "AliPHOSDigitizer.h"
  43 #include "AliPHOSGeometry.h"
  44
  45 // --- Standard library ---
  46
  47
  48 ClassImp(AliPHOSEvalRecPoint)
  49
  50 AliPHOSEvalRecPoint::AliPHOSEvalRecPoint() :
  51   fIsEmc(kFALSE),
  52   fIsCpv(kTRUE),
  53   fParent(-333),
  54   fChi2Dof(-111),
  55   fEventFolderName(AliConfig::GetDefaultEventFolderName())
  56 {
  57   // default ctor
  58 }
  59
  60 AliPHOSEvalRecPoint::AliPHOSEvalRecPoint(Bool_t cpv, AliPHOSEvalRecPoint* parent) :
  61   fIsEmc(kFALSE),
  62   fIsCpv(kFALSE),
  63   fParent(-333),
  64   fChi2Dof(-111),
  65   fEventFolderName("")
  66 {
  67   // ctor
  68   fParent=-333;
  69   fChi2Dof=-111;
  70
  71 //    fParent=parent;
  72   TObjArray* wPool = (TObjArray*)GetWorkingPool();
  73   if(!wPool) {
  74     Fatal("ctor", "Couldn't find working pool") ;
  75   }
  76
  77   fParent = wPool->IndexOf((TObject*)parent);
  78   fChi2Dof = parent->Chi2Dof();
  79
  80   if(cpv) {
  81     fIsEmc = kFALSE;
  82     fIsCpv = kTRUE;
  83   }
  84   else {
  85     fIsEmc = kTRUE;
  86     fIsCpv = kFALSE;
  87   }
  88
  89   // Add itself to working pool
  90   AddToWorkingPool((TObject*)this);
  91
  92 }
  93
  94 AliPHOSEvalRecPoint::AliPHOSEvalRecPoint(Int_t i, Bool_t cpv) :
  95   fIsEmc(kFALSE),
  96   fIsCpv(kFALSE),
  97   fParent(-333),
  98   fChi2Dof(-111),
  99   fEventFolderName(AliConfig::GetDefaultEventFolderName())
 100 {
 101   // ctor
 102   AliPHOSEmcRecPoint* rp=0;
 103
 104   AliPHOSLoader* fLoader = AliPHOSLoader::GetPHOSLoader(fEventFolderName);
 105
 106   if(cpv) {
 107     rp = (AliPHOSCpvRecPoint*)fLoader->CpvRecPoints()->At(i);
 108     fIsEmc = kFALSE;
 109     fIsCpv = kTRUE;
 110   }
 111   else {
 112     rp = (AliPHOSEmcRecPoint*)fLoader->EmcRecPoints()->At(i);
 113     fIsEmc = kTRUE;
 114     fIsCpv = kFALSE;
 115   }
 116
 117   Int_t* digits = rp->GetDigitsList();
 118   Float_t* energies = rp->GetEnergiesList();
 119   Int_t nDigits = rp->GetMultiplicity();
 120
 121   for(Int_t iDigit=0; iDigit<nDigits; iDigit++) {
 122     AliPHOSDigit* digit = (AliPHOSDigit*)fLoader->Digits()->At( digits[iDigit] );
 123     Float_t eDigit = energies[iDigit];
 124     this->AddDigit(*digit,eDigit);
 125   }
 126
 127   TVector3 locpos;
 128   rp->GetLocalPosition(locpos);
 129   fLocPos = locpos;
 130
 131   // Add itself to working pool
 132   AddToWorkingPool((TObject*)this);
 133
 134 }
 135
 136 AliPHOSClusterizer* AliPHOSEvalRecPoint::GetClusterizer()
 137 {
 138   // returns clusterizer task
 139   TFolder* aliceF = (TFolder*)gROOT->FindObjectAny("YSAlice");
 140   TFolder* wPoolF = (TFolder*)aliceF->FindObject("WhiteBoard/RecPoints/PHOS/SmP");
 141   AliPHOSClusterizer* clu = (AliPHOSClusterizer*)wPoolF->FindObject("PHOS:clu-v1");
 142   if(!clu) {
 143     Fatal("GetClusterizer", "Couldn't find Clusterizer") ;
 144   }
 145
 146   return clu;
 147 }
 148
 149 Bool_t AliPHOSEvalRecPoint::TooClose(AliPHOSRecPoint* pt) const
 150 {
 151   // check if a rec.point pt is too close to this one
 152   TVector3 herPos;
 153   TVector3 myPos;
 154   pt->GetLocalPosition(herPos);
 155   this->GetLocalPosition(myPos);
 156   Float_t dx = herPos.X() - myPos.X();
 157   Float_t dz = herPos.Z() - myPos.Z();
 158   Float_t dr = TMath::Sqrt(dx*dx + dz*dz);
 159
 160   if(dr<GetReconstructionManager()->MergeGammasMinDistanceCut())
 161     return kTRUE;
 162   else
 163     return kFALSE;
 164
 165 }
 166
 167 Bool_t AliPHOSEvalRecPoint::NeedToSplit() const
 168 {
 169   // rec.point needs to be split
 170   return kFALSE;
 171 }
 172
 173 void AliPHOSEvalRecPoint::DeleteParent()
 174 {
 175   // delete parent
 176   fParent=-333;
 177 }
 178
 179 void AliPHOSEvalRecPoint::UpdateWorkingPool()
 180 {
 181   // update pool of rec.points
 182   Int_t i; //loop variable
 183
 184   for(i=0; i<this->InWorkingPool(); i++) {
 185      AliPHOSEvalRecPoint* parent = (AliPHOSEvalRecPoint*)GetFromWorkingPool(i);
 186      TObjArray children;
 187      Int_t nChild = parent->HasChild(children);
 188      for(Int_t iChild=0; iChild<nChild; iChild++)
 189        {
 190          ((AliPHOSEvalRecPoint*)children.At(iChild))->DeleteParent();
 191        }
 192
 193      if(nChild) {
 194        RemoveFromWorkingPool(parent);
 195        delete parent;
 196      }
 197
 198   }
 199
 200   for(i=0; i<this->InWorkingPool(); i++) {
 201     AliPHOSEvalRecPoint* weak = (AliPHOSEvalRecPoint*)GetFromWorkingPool(i);
 202     if (weak->KillWeakPoint()) delete weak;
 203   }
 204
 205   for(i=0; i<this->InWorkingPool(); i++) {
 206     AliPHOSEvalRecPoint* close = (AliPHOSEvalRecPoint*)GetFromWorkingPool(i);
 207     close->MergeClosePoint();
 208   }
 209
 210   for(i=0; i<this->InWorkingPool(); i++)
 211     ((AliPHOSEvalRecPoint*)AliPHOSEvalRecPoint::GetFromWorkingPool(i))->SetIndexInList(i);
 212
 213 }
 214
 215 Int_t AliPHOSEvalRecPoint::HasChild(TObjArray& children)
 216 {
 217   // returns the number of children
 218   for(Int_t iChild=0; iChild<InWorkingPool(); iChild++)
 219     {
 220       AliPHOSEvalRecPoint* child = (AliPHOSEvalRecPoint*)GetFromWorkingPool(iChild);
 221       TObject* parent = (TObject*)child->Parent();
 222       TObject* me = (TObject*)this;
 223       if(parent==me) children.Add(child);
 224     }
 225
 226   return children.GetEntries();
 227 }
 228
 229 void AliPHOSEvalRecPoint::Init()
 230 {
 231   // initialization
 232   AliPHOSClusterizer* clusterizer = GetClusterizer();
 233   if(!clusterizer) {
 234     Fatal("Init", "Cannot get clusterizer") ;
 235   }
 236
 237   TClonesArray* digits = AliPHOSLoader::GetPHOSLoader(fEventFolderName)->Digits();
 238   Float_t logWeight=0;
 239
 240   if(this->IsEmc()) {
 241     logWeight = clusterizer->GetEmcLogWeight();
 242   }
 243   else {
 244     logWeight = clusterizer->GetCpvLogWeight();
 245   }
 246
 247   TVector3 vtx(0.,0.,0.) ;
 248   TVector3 inc(0.,0.,0.) ;
 249   EvalLocalPosition(logWeight,vtx,digits,inc); // evaluate initial position
 250 }
 251
 252
 253 void AliPHOSEvalRecPoint::MakeJob()
 254 {
 255   // Reconstruction algoritm implementation.
 256
 257   AliPHOSRecManager* recMng = GetReconstructionManager();
 258
 259   Init();
 260
 261   UnfoldLocalMaxima();
 262
 263   TObjArray children;
 264   Int_t nChild = HasChild(children);
 265
 266   if(!nChild) {
 267     EvaluatePosition();
 268     if(Chi2Dof()>recMng->OneGamChisqCut()) SplitMergedShowers();
 269   }
 270
 271   for(Int_t iChild=0; iChild<nChild; iChild++) {
 272     AliPHOSEvalRecPoint* child = (AliPHOSEvalRecPoint*)children.At(iChild);
 273     child->EvaluatePosition();
 274
 275     if(child->Chi2Dof()>recMng->OneGamChisqCut())
 276       child->SplitMergedShowers();
 277   }
 278
 279 }
 280
 281 void AliPHOSEvalRecPoint::InitTwoGam(Float_t* gamma1, Float_t* gamma2)
 282 {
 283   //Compute start values for two gamma fit algorithm.
 284   // gamma1[0], gamma1[1], gamma1[2] are Energy,x,z of the reconstructed "gammas".
 285
 286
 287   TVector3 lpos; // start point choosing.
 288   GetLocalPosition(lpos);
 289   Float_t xx = lpos.Z();
 290   Float_t yy = lpos.X();
 291   Float_t e = GetEnergy();
 292
 293   AliInfo(Form("(x,z,e)[old] = (%f, %f, %f)", yy, xx, e)) ;
 294
 295 //    xx = XY(xx/e);
 296 //    yy = XY(yy/e);
 297
 298
 299   Float_t eDigit ;
 300   AliPHOSDigit * digit ;
 301   Int_t nDigits = GetMultiplicity();
 302   Int_t * digits = GetDigitsList();
 303   Float_t * energies = GetEnergiesList();
 304
 305   Float_t ix ;
 306   Float_t iy ;
 307   Int_t relid[4] ;
 308
 309   Float_t exx = 0;
 310   Float_t eyy = 0;
 311   Float_t exy = 0;
 312   Float_t sqr;
 313   Float_t cos2fi = 1.;
 314
 315   AliPHOSLoader* fLoader = AliPHOSLoader::GetPHOSLoader(fEventFolderName);
 316   AliPHOSGeometry * phosgeom =  AliPHOSGeometry::GetInstance() ;
 317
 318   Int_t iDigit; //loop variable
 319
 320   for(iDigit = 0 ; iDigit < nDigits ; iDigit ++)
 321     {
 322       digit = (AliPHOSDigit*)fLoader->Digits()->At( digits[iDigit] );
 323       eDigit = energies[iDigit];
 324       phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
 325       phosgeom->RelPosInModule(relid, iy, ix);
 326
 327       Float_t dx =  ix - xx;
 328       Float_t dy =  iy - yy;
 329       exx += eDigit*dx*dx;
 330       eyy += eDigit*dy*dy;
 331       exy += eDigit*dx*dy;
 332
 333     }
 334
 335   sqr = TMath::Sqrt(4.*exy*exy + (exx-eyy)*(exx-eyy));
 336   Float_t euu = (exx+eyy+sqr)/2.;
 337   if(sqr>1.e-10) cos2fi = (exx - eyy)/sqr;
 338   Float_t cosfi = TMath::Sqrt((1.+cos2fi)/2.);
 339   Float_t sinfi = TMath::Sqrt((1.-cos2fi)/2.);
 340   if(exy<0) sinfi = -sinfi;
 341
 342   Float_t eu3 = 0;
 343   for(iDigit = 0 ; iDigit < nDigits ; iDigit ++)
 344     {
 345       digit = (AliPHOSDigit*)fLoader->Digits()->At( digits[iDigit] );
 346       eDigit = energies[iDigit];
 347       phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
 348       phosgeom->RelPosInModule(relid, iy, ix);
 349
 350       Float_t dx =  ix - xx;
 351       Float_t dy =  iy - yy;
 352       Float_t du = dx*cosfi + dy*sinfi;
 353       eu3 += eDigit*du*du*du;
 354     }
 355
 356   Float_t c = e*eu3*eu3/(euu*euu*euu)/2.;
 357   Float_t sign = 1.;
 358   if(eu3<0) sign = -1.;
 359   Float_t r = 1.+ c + sign*TMath::Sqrt(2.*c + c*c);
 360   Float_t u = 0;
 361   if(TMath::Abs(r-1.)>0.1) u = eu3/euu*(r+1.)/(r-1.);
 362   if(TMath::Abs(r-1.)<0.1) u = TMath::Sqrt(sqr/e/r)*(1.+r);
 363
 364   Float_t e2c = e/(1.+r);
 365   Float_t e1c = e-e2c;
 366   Float_t u1 = -u/(1.+r);
 367   Float_t u2 = u+u1;
 368
 369   Float_t x1c = xx+u1*cosfi;
 370   Float_t y1c = yy+u1*sinfi;
 371   Float_t x2c = xx+u2*cosfi;
 372   Float_t y2c = yy+u2*sinfi;
 373
 374 //    printf("e1c -> %f\n",e1c);
 375 //    printf("x1c -> %f\n",x1c);
 376 //    printf("y1c -> %f\n",y1c);
 377 //    printf("e2c -> %f\n",e2c);
 378 //    printf("x2c -> %f\n",x2c);
 379 //    printf("y2c -> %f\n",y2c);
 380
 381   gamma1[0] = e1c;
 382   gamma1[1] = y1c;
 383   gamma1[2] = x1c;
 384
 385   gamma2[0] = e2c;
 386   gamma2[1] = y2c;
 387   gamma2[2] = x2c;
 388
 389 }
 390
 391 void AliPHOSEvalRecPoint::TwoGam(Float_t* gamma1, Float_t* gamma2)
 392 {
 393   //Fitting algorithm for the two very closed gammas
 394   //that merged into the cluster with one maximum.
 395   // Starting points gamma1 and gamma2 must be provided before ( see InitTwoGam)
 396   //Set chisquare of the fit.
 397
 398
 399   Float_t st0 = GetReconstructionManager()->TwoGamInitialStep();
 400   Float_t chmin = GetReconstructionManager()->TwoGamChisqMin();
 401   Float_t emin = GetReconstructionManager()->TwoGamEmin();
 402   Float_t stpmin = GetReconstructionManager()->TwoGamStepMin();
 403   Int_t nIter = GetReconstructionManager()->TwoGamNumOfIterations(); // Number of iterations.
 404
 405   Float_t chisq = 100.; //Initial chisquare.
 406
 407   Int_t nadc = GetMultiplicity();
 408   if(nadc<3)
 409       fChi2Dof= -111.;
 410
 411   Int_t dof = nadc - 5;
 412   if(dof<1) dof=1;
 413   Float_t chstop = chmin*dof;
 414
 415   Float_t ch = 1.e+20;
 416   Float_t st = st0;
 417   Float_t grx1 = 0.;
 418   Float_t gry1 = 0.;
 419   Float_t grx2 = 0.;
 420   Float_t gry2 = 0.;
 421   Float_t gre = 0.;
 422   Float_t gr = 1.;
 423   Float_t ee1=0,xx1=0,yy1=0,ee2=0,xx2=0,yy2=0,cosi;
 424
 425   Float_t e1c = gamma1[0];
 426   Float_t y1c = gamma1[1];
 427   Float_t x1c = gamma1[2];
 428
 429   Float_t e2c = gamma2[0];
 430   Float_t y2c = gamma2[1];
 431   Float_t x2c = gamma2[2];
 432
 433   Float_t e = GetEnergy();
 434
 435   Float_t eDigit ;
 436   AliPHOSDigit * digit ;
 437   Int_t nDigits = GetMultiplicity();
 438   Int_t * digits = GetDigitsList();
 439   Float_t * energies = GetEnergiesList();
 440
 441   Float_t ix ;
 442   Float_t iy ;
 443   Int_t relid[4] ;
 444
 445   AliPHOSLoader* fLoader = AliPHOSLoader::GetPHOSLoader(fEventFolderName);
 446   AliPHOSGeometry * phosgeom =  AliPHOSGeometry::GetInstance() ;
 447
 448   for(Int_t iter=0; iter<nIter; iter++)
 449     {
 450       Float_t chisqc = 0.;
 451       Float_t grx1c = 0.;
 452       Float_t gry1c = 0.;
 453       Float_t grx2c = 0.;
 454       Float_t gry2c = 0.;
 455       Float_t grec = 0.;
 456
 457       for(Int_t iDigit = 0 ; iDigit < nDigits ; iDigit ++)
 458         {
 459           digit = (AliPHOSDigit*)fLoader->Digits()->At( digits[iDigit] );
 460           eDigit = energies[iDigit];
 461           phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
 462           phosgeom->RelPosInModule(relid, iy, ix);
 463
 464           Float_t a1,gx1,gy1;
 465           Float_t a2,gx2,gy2;
 466
 467           Float_t dx1 =  x1c - ix;
 468           Float_t dy1 =  y1c - iy;
 469
 470           GetReconstructionManager()->AG(e1c,dx1,dy1,a1,gx1,gy1);
 471
 472           Float_t dx2 =  x2c - ix;
 473           Float_t dy2 =  y2c - iy;
 474
 475           GetReconstructionManager()->AG(e2c,dx2,dy2,a2,gx2,gy2);
 476
 477           Float_t a = a1+a2;
 478 //        Float_t D = Const*a*(1. - a/e);
 479 //        if(D<0) D=0;
 480
 481 //  //            D = 9.; // ????
 482
 483 //        Float_t da = a - eDigit;
 484 //        chisqc += da*da/D;
 485 //        Float_t dd = da/D;
 486 //        dd = dd*(2.-dd*Const*(1.-2.*a/e));
 487
 488           Float_t dd;
 489           chisqc += GetReconstructionManager()->TwoGamChi2(a,eDigit,e,dd);
 490
 491           grx1c += gx1*dd;
 492           gry1c += gy1*dd;
 493           grx2c += gx2*dd;
 494           gry2c += gy2*dd;
 495           grec += (a1/e1c - a2/e2c)*e*dd;
 496
 497         }
 498
 499       Float_t grc = TMath::Sqrt(grx1c*grx1c + gry1c*gry1c + grx2c*grx2c + gry2c*gry2c);
 500       if(grc<1.e-10) grc=1.e-10;
 501
 502       Float_t sc = 1. + chisqc/ch;
 503       st = st/sc;
 504
 505       if(chisqc>ch)
 506         goto loop20;
 507       cosi = (grx1*grx1c + gry1*gry1c + grx2*grx2c + gry2*gry2c + gre*grec)/gr/grc;
 508       st = st*sc/(1.4 - cosi);
 509
 510       ee1 = e1c;
 511       xx1 = x1c;
 512       yy1 = y1c;
 513       ee2 = e2c;
 514       xx2 = x2c;
 515       yy2 = y2c;
 516
 517       ch = chisqc;
 518
 519       if(ch<chstop)
 520         goto loop101;
 521
 522       grx1 = grx1c;
 523       gry1 = gry1c;
 524       grx2 = grx2c;
 525       gry2 = gry2c;
 526       gre = grec;
 527       gr = grc;
 528
 529     loop20: ;
 530       Float_t step = st*gr;
 531
 532       AliInfo(Form("Iteration %d dof %d chisq/dof %f chstop/dof %f step %d stpmin %d",
 533            iter, dof, ch/dof, chstop/dof, step, stpmin)) ;
 534
 535
 536       if(step<stpmin)
 537         goto loop101;
 538
 539       Float_t de = st*gre*e;
 540       e1c = ee1 - de;
 541       e2c = ee2 + de;
 542
 543       if( (e1c>emin) && (e2c>emin) )
 544         goto loop25;
 545       st = st/2.;
 546       goto loop20;
 547
 548     loop25: ;
 549       x1c = xx1 - st*grx1;
 550       y1c = yy1 - st*gry1;
 551       x2c = xx2 - st*grx2;
 552       y2c = yy2 - st*gry2;
 553
 554
 555     }
 556
 557  loop101:
 558
 559 //    if(ch>chisq*(nadc-2)-delch)
 560 //      return ch/dof;
 561
 562   chisq = ch/dof;
 563   gamma1[0] = ee1;
 564   gamma1[1] = yy1;
 565   gamma1[2] = xx1;
 566
 567   gamma2[0] = ee2;
 568   gamma2[1] = yy2;
 569   gamma2[2] = xx2;
 570
 571   Float_t x1New = yy1;
 572   Float_t z1New = xx1;
 573   Float_t e1New = ee1;
 574   Float_t x2New = yy2;
 575   Float_t z2New = xx2;
 576   Float_t e2New = ee2;
 577
 578   TString message ;
 579   message  = "     (x,z,e)[1 fit] = (%f, %f, %f)\n" ;
 580   message  = "     (x,z,e)[2 fit] = (%f, %f, %f)\n" ;
 581   AliInfo(Form(message.Data(),
 582        x1New, z1New, e1New,
 583        x2New, z2New, e2New)) ;
 584
 585   fChi2Dof = chisq;
 586
 587 }
 588
 589 void AliPHOSEvalRecPoint::UnfoldTwoMergedPoints(Float_t* gamma1, Float_t* gamma2)
 590 {
 591   //Unfold TWO merged rec. points in the case when cluster has only one maximum,
 592   //but it's fitting to the one gamma shower is too bad.
 593   // Use TwoGam() to estimate the positions and energies of merged points.
 594
 595
 596   Int_t nMax = 2;
 597   Float_t* gamma;
 598
 599   Int_t* digits = GetDigitsList();
 600   Int_t nDigits = GetMultiplicity();
 601   Float_t* energies = GetEnergiesList();
 602   Float_t* eFit = new Float_t[nDigits];
 603
 604   AliPHOSLoader* fLoader = AliPHOSLoader::GetPHOSLoader(fEventFolderName);
 605   AliPHOSGeometry * phosgeom =  AliPHOSGeometry::GetInstance() ;
 606
 607   for(Int_t iDigit=0; iDigit<nDigits; iDigit++)
 608     {
 609       AliPHOSDigit* digit = (AliPHOSDigit*)fLoader->Digits()->At( digits[iDigit] );
 610       Int_t relid[4] ;
 611       phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
 612       Float_t x,z;
 613       phosgeom->RelPosInModule(relid, x, z);
 614
 615       Float_t gain = 0.;
 616       for(Int_t iMax=0; iMax<nMax; iMax++)
 617         {
 618           if(iMax==0)
 619             gamma = gamma1;
 620           else
 621             gamma = gamma2;
 622
 623           Float_t eMax = gamma[0];
 624           Float_t xMax = gamma[1];
 625           Float_t zMax = gamma[2];
 626
 627           Float_t dx = xMax - x;
 628           Float_t dz = zMax - z;
 629           Float_t amp,gx,gy;
 630           GetReconstructionManager()->AG(eMax,dz,dx,amp,gx,gy);
 631           gain += amp;
 632         }
 633
 634       eFit[iDigit] = gain;
 635     }
 636
 637
 638   for( Int_t iMax=0; iMax<nMax; iMax++)
 639     {
 640       if(iMax==0)
 641         gamma = gamma1;
 642       else
 643         gamma = gamma2;
 644
 645       Float_t eMax = gamma[0];
 646       Float_t xMax = gamma[1];
 647       Float_t zMax = gamma[2];
 648
 649       AliPHOSEvalRecPoint* newRP = new AliPHOSEvalRecPoint(IsCPV(),this);
 650       TVector3 newpos(xMax,0,zMax);
 651       newRP->SetLocalPosition(newpos);
 652
 653       for( Int_t iDigit = 0 ; iDigit < nDigits ; iDigit ++)
 654         {
 655           AliPHOSDigit* digit = (AliPHOSDigit*)fLoader->Digits()->At( digits[iDigit] );
 656           Float_t eDigit = energies[iDigit];
 657           Int_t relid[4] ;
 658           phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
 659           Float_t ix,iz;
 660           phosgeom->RelPosInModule(relid, ix, iz);
 661
 662           Float_t dx = xMax - ix;
 663           Float_t dz = zMax - iz;
 664           Float_t singleShowerGain,gxMax,gyMax;
 665           GetReconstructionManager()->AG(eMax,dz,dx,singleShowerGain,gxMax,gyMax);
 666           Float_t totalGain = eFit[iDigit];
 667           Float_t ratio = singleShowerGain/totalGain;
 668           eDigit = eDigit*ratio;
 669           newRP->AddDigit(*digit,eDigit);
 670         }
 671       AliInfo(Form("======= Split: daughter rec point %d =================",
 672                    iMax)) ;
 673       newRP->Print();
 674
 675     }
 676
 677   delete[] eFit;
 678
 679
 680 }
 681
 682
 683 void AliPHOSEvalRecPoint::EvaluatePosition()
 684 {
 685   // One gamma fit algorithm.
 686   // Set chisq/dof of the fit.
 687   // gamma1[0], gamma1[1], gamma1[2] are Energy,x,z of the reconstructed gamma, respectively.
 688
 689   Int_t nDigits = GetMultiplicity();
 690   if(nDigits<2)
 691     return;
 692
 693   Int_t nIter = GetReconstructionManager()->OneGamNumOfIterations(); // number of iterations
 694   Float_t st0 = GetReconstructionManager()->OneGamInitialStep(); // initial step
 695 //    const Float_t stpMin = 0.005;
 696   Float_t stpMin = GetReconstructionManager()->OneGamStepMin();
 697   Float_t chmin =  GetReconstructionManager()->OneGamChisqMin();
 698
 699   TVector3 locpos;
 700   AliPHOSDigit* digit;
 701   Float_t eDigit;
 702   Int_t relid[4] ;
 703   Float_t ix, iy;
 704
 705   GetLocalPosition(locpos);
 706   Float_t e = GetEnergy();
 707   Float_t xc = locpos.Z();
 708   Float_t yc = locpos.X();
 709   Float_t dx,dy,gx,gy,grxc,gryc;
 710   Float_t st = st0;
 711   Float_t chisq = 1.e+20;
 712   Float_t gr = 1.;
 713   Float_t grx = 0.;
 714   Float_t gry = 0.;
 715   Int_t dof = GetMultiplicity() - 2;
 716   if(dof<1)
 717     dof = 1;
 718   Float_t chstop = chmin*dof;
 719   Float_t cosi,x1=0,y1=0;
 720   Float_t chisqc;
 721
 722   AliPHOSLoader* fLoader = AliPHOSLoader::GetPHOSLoader(fEventFolderName);
 723   AliPHOSGeometry * phosgeom =  AliPHOSGeometry::GetInstance() ;
 724
 725   for(Int_t iter=0; iter<nIter; iter++)
 726     {
 727
 728       chisqc = 0.;
 729       grxc = 0.;
 730       gryc = 0.;
 731       for(Int_t iDigit = 0 ; iDigit < nDigits ; iDigit ++)
 732         {
 733
 734           Float_t* energies = GetEnergiesList();
 735           Int_t* digits = GetDigitsList();
 736           digit = (AliPHOSDigit*)fLoader->Digits()->At( digits[iDigit] );
 737           eDigit = energies[iDigit];
 738           phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
 739           phosgeom->RelPosInModule(relid, iy, ix);
 740
 741           dx =  xc - ix;
 742           dy =  yc - iy;
 743
 744           if(!dx) dx=dy;
 745           if(!dy) dy=dx;
 746
 747           Float_t a;
 748           GetReconstructionManager()->AG(e,dx,dy,a,gx,gy);
 749           Float_t dd;
 750           AliInfo(Form("  (ix iy  xc yc  dx dy)   %f %f %f %f %f %f",
 751                        ix, iy, xc, yc, dx, dy)) ;
 752           Float_t chi2dg = GetReconstructionManager()->OneGamChi2(a,eDigit,e,dd);
 753
 754           // Exclude digit with too large chisquare.
 755           if(chi2dg > 10) {  continue; }
 756
 757           chisqc += chi2dg;
 758           grxc += gx*dd;
 759           gryc += gy*dd;
 760
 761         }
 762
 763       Float_t grc = TMath::Sqrt(grxc*grxc + gryc*gryc);
 764       if(grc<1.e-10) grc=1.e-10;
 765       Float_t sc = 1. + chisqc/chisq;
 766        AliInfo(Form(" chisq: %f", chisq)) ;
 767       st = st/sc;
 768       if(chisqc>chisq)
 769         goto loop20;
 770       cosi = (grx*grxc + gry*gryc)/gr/grc;
 771       st = st*sc/(1.4 - cosi);
 772       x1 = xc;
 773       y1 = yc;
 774       grx = grxc;
 775       gry = gryc;
 776       gr = grc;
 777       chisq = chisqc;
 778
 779       if(chisq<chstop)
 780         goto loop101;
 781
 782     loop20: ;
 783       Float_t step = st*gr;
 784
 785       AliInfo(Form(" Iteration %d dof %d chisq/dof %f chstop/dof %f step %d stpMin %d",
 786            iter, dof, chisq/dof, chisq/dof, chstop/dof, step, stpMin)) ;
 787
 788
 789       if(step<stpMin)
 790         goto loop101;
 791
 792       if(step>1.)
 793         st = 1./gr;
 794       xc = x1 - st*grx;
 795       yc = y1 - st*gry;
 796     }
 797
 798
 799  loop101:
 800   chisq = chisq/dof;
 801 //    if(chisq>Chsqcut)
 802 //      {
 803 //  //        TwoGam();
 804 //      }
 805
 806   Float_t gamma1[3];
 807   gamma1[0] = e;
 808   gamma1[1] = y1;
 809   gamma1[2] = x1;
 810
 811   TVector3 newpos(gamma1[1],0,gamma1[2]);
 812   //SetLocalPosition(newpos);
 813
 814   fLocPos=newpos;
 815   fAmp=e;
 816
 817 //    TVector3 pos;
 818 //    RP->GetLocalPosition(pos);
 819
 820
 821
 822   fChi2Dof = chisq;
 823
 824 }
 825
 826
 827 Bool_t AliPHOSEvalRecPoint::KillWeakPoint()
 828 {
 829   // Remove this point from further procession
 830   // if it's energy is too small.
 831
 832   Float_t thr0 = GetReconstructionManager()->KillGamMinEnergy();
 833
 834   if(GetEnergy()<thr0) {
 835     AliInfo(Form("+++++++ Killing this rec point ++++++++++")) ;
 836     RemoveFromWorkingPool(this);
 837     return kTRUE;
 838   }
 839
 840   return kFALSE;
 841 }
 842
 843
 844 void AliPHOSEvalRecPoint::SplitMergedShowers()
 845 {
 846   // Split two very close rec. points.
 847
 848   if(GetMultiplicity()<2)
 849     return;
 850
 851   Float_t gamma1[3];
 852   Float_t gamma2[3];
 853
 854   InitTwoGam(gamma1,gamma2); // start points for fit
 855   TwoGam(gamma1,gamma2); // evaluate the positions and energies
 856   UnfoldTwoMergedPoints(gamma1,gamma2); // make the job
 857
 858 }
 859
 860
 861 void AliPHOSEvalRecPoint::MergeClosePoint()
 862 {
 863   // merge rec.points if they are too close
 864   AliPHOSLoader* fLoader = AliPHOSLoader::GetPHOSLoader(fEventFolderName);
 865 //    AliPHOSDigitizer* digitizer = fLoader->Digitizer();
 866 //    Float_t fPedestal = digitizer->GetPedestal();  // YVK 30.09.2001
 867 //    Float_t fSlope = digitizer->GetSlope();
 868
 869   for (Int_t i=0;i<InWorkingPool(); i++)
 870     {
 871       TObject* obj = GetFromWorkingPool(i);
 872       if(!((TObject*)this)->IsEqual(obj))
 873         {
 874           AliPHOSRecPoint* rp = (AliPHOSRecPoint*)obj;
 875           if(GetPHOSMod() == rp->GetPHOSMod())
 876             {
 877               if(TooClose(rp))
 878                 {
 879                   AliInfo(Form("+++++++ Merging point 1: ++++++")) ;
 880                   this->Print();
 881                   AliInfo(Form("+++++++ and point 2: ++++++++++")) ;
 882                   ((AliPHOSEvalRecPoint*)rp)->Print();
 883
 884                   //merge two rec. points
 885                   TVector3 lpos1;
 886                   TVector3 lpos2;
 887                   this->GetLocalPosition(lpos1);
 888                   rp->GetLocalPosition(lpos2);
 889
 890                   Int_t* digits = rp->GetDigitsList();
 891                   Float_t dE = rp->GetEnergy()/(rp->GetEnergy()+this->GetEnergy());
 892                   Float_t newX = lpos1.X()*dE + lpos2.X()*(1.-dE);
 893                   Float_t newZ = lpos1.Z()*dE + lpos2.Z()*(1.-dE);
 894                   Float_t newE = rp->GetEnergy()+this->GetEnergy();
 895                   Float_t* energies = ((AliPHOSEmcRecPoint*)rp)->GetEnergiesList();
 896
 897                   for(Int_t iDigit=0; iDigit<rp->GetDigitsMultiplicity(); iDigit++)
 898                     {
 899                       AliPHOSDigit* digit = (AliPHOSDigit*)fLoader->Digits()->At(digits[iDigit]);
 900                       Float_t eDigit = energies[iDigit];
 901                       this->AddDigit(*digit,eDigit);
 902                     }
 903
 904                   TVector3 newpos(newX,0,newZ);
 905                   fLocPos = newpos;
 906                   fAmp = newE;
 907                   RemoveFromWorkingPool(rp);
 908                   delete rp;
 909
 910                   AliInfo(Form("++++++ Resulting point: ++++++++")) ;
 911                   this->Print();
 912
 913                   break;
 914                 }
 915             }
 916         }
 917     }
 918
 919 }
 920
 921 Int_t AliPHOSEvalRecPoint::UnfoldLocalMaxima()
 922 {
 923   // Make unfolding in the reconstruction point with several local maxima.
 924   // Return the number of local maxima.
 925
 926   // if multiplicity less then 2 - nothing to unfold
 927   if(GetMultiplicity()<2) return 1;
 928
 929   AliPHOSDigit * maxAt[1000];
 930   Float_t maxAtEnergy[1000];
 931   Float_t locMaxCut, logWeight;
 932   Int_t relid[4] ;
 933   Float_t xMax;
 934   Float_t zMax;
 935
 936   AliPHOSClusterizer* clusterizer = GetClusterizer();
 937   if(!clusterizer) {
 938     AliFatal(Form("Cannot get clusterizer")) ;
 939   }
 940
 941   if(this->IsEmc()) {
 942     locMaxCut = clusterizer->GetEmcLocalMaxCut();
 943     logWeight = clusterizer->GetEmcLogWeight();
 944   }
 945   else {
 946     locMaxCut = clusterizer->GetCpvLocalMaxCut();
 947     logWeight = clusterizer->GetCpvLogWeight();
 948   }
 949
 950   AliPHOSLoader* fLoader = AliPHOSLoader::GetPHOSLoader(fEventFolderName);
 951   AliPHOSGeometry * phosgeom =  AliPHOSGeometry::GetInstance() ;
 952   TClonesArray* digits = fLoader->Digits();
 953
 954   // if number of local maxima less then 2 - nothing to unfold
 955   Int_t nMax = GetNumberOfLocalMax(maxAt,maxAtEnergy,locMaxCut,digits);
 956   if(nMax<2) return nMax;
 957
 958   Int_t* digitsList = GetDigitsList();
 959   Int_t nDigits = GetMultiplicity();
 960   Float_t* energies = GetEnergiesList();
 961   Float_t* eFit = new Float_t[nDigits];
 962
 963   Int_t iDigit; //loop variable
 964
 965   for(iDigit=0; iDigit<nDigits; iDigit++)
 966     {
 967       eFit[iDigit] = 0.;
 968     }
 969
 970   for(iDigit=0; iDigit<nDigits; iDigit++)
 971     {
 972
 973       AliPHOSDigit* digit = (AliPHOSDigit*)fLoader->Digits()->At( digitsList[iDigit] );
 974       phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
 975       Float_t x,z;
 976       phosgeom->RelPosInModule(relid, x, z);
 977
 978       for(Int_t iMax=0; iMax<nMax; iMax++)
 979         {
 980           AliPHOSDigit* digitMax = maxAt[iMax];
 981           Float_t eMax = maxAtEnergy[iMax];
 982           phosgeom->AbsToRelNumbering(digitMax->GetId(), relid) ;
 983           phosgeom->RelPosInModule(relid, xMax, zMax);
 984           Float_t dx = xMax - x;
 985           Float_t dz = zMax - z;
 986           Float_t amp,gx,gy;
 987           GetReconstructionManager()->AG(eMax,dz,dx,amp,gx,gy);
 988 //        amp = amp + 0.5;
 989           eFit[iDigit] += amp;
 990         }
 991     }
 992
 993
 994   for(Int_t iMax=0; iMax<nMax; iMax++)
 995     {
 996       AliPHOSDigit* digitMax = maxAt[iMax];
 997       phosgeom->AbsToRelNumbering(digitMax->GetId(), relid) ;
 998       phosgeom->RelPosInModule(relid, xMax, zMax);
 999       Float_t eMax = maxAtEnergy[iMax];
1000
1001       AliPHOSEvalRecPoint* newRP = new AliPHOSEvalRecPoint(IsCPV(),this);
1002       newRP->AddDigit(*digitMax,maxAtEnergy[iMax]);
1003
1004       //Neighbous ( matrix 3x3 around the local maximum)
1005       for(iDigit=0; iDigit<nDigits; iDigit++)
1006         {
1007           AliPHOSDigit* digit = (AliPHOSDigit*)fLoader->Digits()->At( digitsList[iDigit] );
1008           Float_t eDigit = energies[iDigit];
1009           phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
1010           Float_t ix,iz;
1011           phosgeom->RelPosInModule(relid, ix, iz);
1012
1013           if(AreNeighbours(digitMax,digit))
1014             {
1015               Float_t dx = xMax - ix;
1016               Float_t dz = zMax - iz;
1017               Float_t singleShowerGain,gxMax,gyMax;
1018               GetReconstructionManager()->AG(eMax,dz,dx,singleShowerGain,gxMax,gyMax);
1019               Float_t totalGain = eFit[iDigit];
1020               Float_t ratio = singleShowerGain/totalGain;
1021               AliInfo(Form(" ratio -> %f", ratio)) ;
1022               eDigit = eDigit*ratio;
1023               newRP->AddDigit(*digit,eDigit);
1024             }
1025         }
1026
1027       TVector3 vtx(0.,0.,0.) ;
1028       TVector3 inc(0.,0.,0.) ;
1029       newRP->EvalLocalPosition(logWeight,vtx,digits,inc);
1030       AliInfo(Form("======= Unfold: daughter rec point %d =================",
1031                    iMax)) ;
1032       newRP->Print();
1033     }
1034
1035 //    RemoveFromWorkingPool(this);
1036
1037   delete[] eFit;
1038
1039   return nMax;
1040
1041 }
1042
1043 void AliPHOSEvalRecPoint::PrintPoint()
1044 {
1045   // print rec.point to stdout
1046   AliPHOSCpvRecPoint::Print();
1047
1048   TVector3 lpos;
1049   GetLocalPosition(lpos);
1050
1051   TString message ;
1052   message  = "       Chi2/dof = %f" ;
1053   message += "       Local (x,z) = (%f, %f) in module %d" ;
1054   AliInfo(Form(message.Data(), Chi2Dof(), lpos.X(), lpos.Z(), GetPHOSMod())) ;
1055
1056 }
1057
1058 AliPHOSRecManager* AliPHOSEvalRecPoint::GetReconstructionManager() const
1059 {
1060   // returns reconstruction manager
1061   TFolder* aliceF = (TFolder*)gROOT->FindObjectAny("YSAlice");
1062   TFolder* wPoolF = (TFolder*)aliceF->FindObject("WhiteBoard/RecPoints/PHOS/SmP");
1063   AliPHOSRecManager* recMng = (AliPHOSRecManager*)wPoolF->FindObject("AliPHOSRecManager");
1064   if(!recMng) {
1065     AliFatal(Form("Couldn't find Reconstruction Manager")) ;
1066   }
1067
1068   return recMng;
1069 }
1070
1071
1072 AliPHOSRecPoint* AliPHOSEvalRecPoint::Parent()
1073 {
1074   // returns a parent
1075   if(fParent<0) return NULL;
1076   else
1077     return (AliPHOSRecPoint*)GetFromWorkingPool(fParent);
1078 //    return fParent;
1079 }
1080
1081 Float_t AliPHOSEvalRecPoint::Chi2Dof() const
1082 {
1083   // returns a chi^2 per degree of freedom
1084   return fChi2Dof;
1085 }
1086
1087 const TObject* AliPHOSEvalRecPoint::GetWorkingPool()
1088 {
1089   // returns a pool of rec.points
1090   TFolder* aliceF = (TFolder*)gROOT->FindObjectAny("YSAlice");
1091   TFolder* wPoolF = (TFolder*)aliceF->FindObject("WhiteBoard/RecPoints/PHOS/SmP");
1092   TObject* wPool = wPoolF->FindObject("SmartPoints");
1093   if(!wPool) {
1094     AliFatal(Form("Couldn't find Working Pool")) ;
1095   }
1096
1097   return wPool;
1098 }
1099
1100 void AliPHOSEvalRecPoint::AddToWorkingPool(TObject* obj)
1101 {
1102   // add a rec.point to a pool
1103   ((TObjArray*)GetWorkingPool())->Add(obj);
1104 }
1105
1106 TObject* AliPHOSEvalRecPoint::GetFromWorkingPool(Int_t i)
1107 {
1108   // return a rec.point from a pool at an index i
1109 //    return fWorkingPool->At(i);
1110   return ((TObjArray*)GetWorkingPool())->At(i);
1111 }
1112
1113 Int_t AliPHOSEvalRecPoint::InWorkingPool()
1114 {
1115   // return the number of rec.points in a pool
1116   return ((TObjArray*)GetWorkingPool())->GetEntries();
1117 }
1118
1119 void AliPHOSEvalRecPoint::RemoveFromWorkingPool(TObject* obj)
1120 {
1121   // remove a rec.point from a pool
1122   ((TObjArray*)GetWorkingPool())->Remove(obj);
1123   ((TObjArray*)GetWorkingPool())->Compress();
1124 }
1125
1126 void AliPHOSEvalRecPoint::PrintWorkingPool()
1127 {
1128   // print pool of rec.points to stdout
1129   ((TObjArray*)GetWorkingPool())->Print();
1130 }