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