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