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