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