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