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