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