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