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