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