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