PHOS/AliPHOSEvalRecPoint.cxx

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