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