1 // ---- ROOT system ---
2 #include "TDirectory.h"
8 // --- AliRoot header files ---
9 #include "AliPHOSEvalRecPoint.h"
11 #include "AliPHOSGetter.h"
12 #include "AliPHOSRecCpvManager.h"
13 #include "AliPHOSRecEmcManager.h"
14 #include "AliPHOSDigitizer.h"
16 // --- Standard library ---
20 ClassImp(AliPHOSEvalRecPoint)
22 AliPHOSEvalRecPoint::AliPHOSEvalRecPoint(): AliPHOSCpvRecPoint()
30 AliPHOSEvalRecPoint::AliPHOSEvalRecPoint(Bool_t cpv, AliPHOSEvalRecPoint* parent) : AliPHOSCpvRecPoint()
37 TObjArray* wPool = (TObjArray*)GetWorkingPool();
39 cout<<" Couldn't find working pool. Exit."<<endl;
43 fParent = wPool->IndexOf((TObject*)parent);
44 fChi2Dof = parent->Chi2Dof();
55 // Add itself to working pool
56 AddToWorkingPool((TObject*)this);
60 AliPHOSEvalRecPoint::AliPHOSEvalRecPoint(Int_t i, Bool_t cpv) : AliPHOSCpvRecPoint()
66 AliPHOSEmcRecPoint* rp=0;
68 AliPHOSGetter* fGetter = AliPHOSGetter::GetInstance();
71 rp = (AliPHOSCpvRecPoint*)fGetter->CpvRecPoints()->At(i);
76 rp = (AliPHOSEmcRecPoint*)fGetter->EmcRecPoints()->At(i);
81 Int_t* Digits = rp->GetDigitsList();
82 Float_t* Energies = rp->GetEnergiesList();
83 Int_t nDigits = rp->GetMultiplicity();
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);
92 rp->GetLocalPosition(locpos);
95 // Add itself to working pool
96 AddToWorkingPool((TObject*)this);
100 AliPHOSClusterizer* AliPHOSEvalRecPoint::GetClusterizer()
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");
106 cout<<" Couldn't find Clusterizer. Exit."<<endl;
113 Bool_t AliPHOSEvalRecPoint::TooClose(AliPHOSRecPoint* pt) const
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);
123 if(dr<GetReconstructionManager()->MergeGammasMinDistanceCut())
130 Bool_t AliPHOSEvalRecPoint::NeedToSplit() const
135 void AliPHOSEvalRecPoint::DeleteParent()
140 void AliPHOSEvalRecPoint::UpdateWorkingPool()
143 Int_t i; //loop variable
145 for(i=0; i<this->InWorkingPool(); i++) {
146 AliPHOSEvalRecPoint* parent = (AliPHOSEvalRecPoint*)GetFromWorkingPool(i);
148 Int_t nChild = parent->HasChild(children);
149 for(Int_t iChild=0; iChild<nChild; iChild++)
151 ((AliPHOSEvalRecPoint*)children.At(iChild))->DeleteParent();
155 RemoveFromWorkingPool(parent);
161 for(i=0; i<this->InWorkingPool(); i++) {
162 AliPHOSEvalRecPoint* weak = (AliPHOSEvalRecPoint*)GetFromWorkingPool(i);
163 if (weak->KillWeakPoint()) delete weak;
166 for(i=0; i<this->InWorkingPool(); i++) {
167 AliPHOSEvalRecPoint* close = (AliPHOSEvalRecPoint*)GetFromWorkingPool(i);
168 close->MergeClosePoint();
171 for(i=0; i<this->InWorkingPool(); i++)
172 ((AliPHOSEvalRecPoint*)AliPHOSEvalRecPoint::GetFromWorkingPool(i))->SetIndexInList(i);
176 Int_t AliPHOSEvalRecPoint::HasChild(TObjArray& children)
178 for(Int_t iChild=0; iChild<InWorkingPool(); iChild++)
180 AliPHOSEvalRecPoint* child = (AliPHOSEvalRecPoint*)GetFromWorkingPool(iChild);
181 TObject* parent = (TObject*)child->Parent();
182 TObject* me = (TObject*)this;
183 if(parent==me) children.Add(child);
186 return children.GetEntries();
189 void AliPHOSEvalRecPoint::Init()
192 AliPHOSClusterizer* clusterizer = GetClusterizer();
194 cout<<" Cannot get clusterizer. Exit."<<endl;
198 TClonesArray* digits = AliPHOSGetter::GetInstance()->Digits();
202 LogWeight = clusterizer->GetEmcLogWeight();
205 LogWeight = clusterizer->GetCpvLogWeight();
208 EvalLocalPosition(LogWeight,digits); // evaluate initial position
212 void AliPHOSEvalRecPoint::MakeJob()
214 // Reconstruction algoritm implementation.
216 AliPHOSRecManager* recMng = GetReconstructionManager();
223 Int_t nChild = HasChild(children);
227 if(Chi2Dof()>recMng->OneGamChisqCut()) SplitMergedShowers();
230 for(Int_t iChild=0; iChild<nChild; iChild++) {
231 AliPHOSEvalRecPoint* child = (AliPHOSEvalRecPoint*)children.At(iChild);
232 child->EvaluatePosition();
234 if(child->Chi2Dof()>recMng->OneGamChisqCut())
235 child->SplitMergedShowers();
240 void AliPHOSEvalRecPoint::InitTwoGam(Float_t* gamma1, Float_t* gamma2)
242 //Compute start values for two gamma fit algorithm.
243 // gamma1[0], gamma1[1], gamma1[2] are Energy,x,z of the reconstructed "gammas".
246 TVector3 lpos; // start point choosing.
247 GetLocalPosition(lpos);
248 Float_t xx = lpos.Z();
249 Float_t yy = lpos.X();
250 Float_t E = GetEnergy();
252 cout<<" (x,z,E)[old] = ("<<yy<<","<<xx<<","<<E<<")"<<endl;
259 AliPHOSDigit * digit ;
260 Int_t nDigits = GetMultiplicity();
261 Int_t * Digits = GetDigitsList();
262 Float_t * Energies = GetEnergiesList();
274 AliPHOSGetter* fGetter = AliPHOSGetter::GetInstance();
275 const AliPHOSGeometry* fGeom = fGetter->PHOSGeometry();
277 Int_t iDigit; //loop variable
279 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++)
281 digit = (AliPHOSDigit*)fGetter->Digits()->At( Digits[iDigit] );
282 eDigit = Energies[iDigit];
283 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
284 fGeom->RelPosInModule(relid, iy, ix);
286 Float_t dx = ix - xx;
287 Float_t dy = iy - yy;
294 sqr = TMath::Sqrt(4.*exy*exy + (exx-eyy)*(exx-eyy));
295 Float_t euu = (exx+eyy+sqr)/2.;
296 if(sqr>1.e-10) cos2fi = (exx - eyy)/sqr;
297 Float_t cosfi = TMath::Sqrt((1.+cos2fi)/2.);
298 Float_t sinfi = TMath::Sqrt((1.-cos2fi)/2.);
299 if(exy<0) sinfi = -sinfi;
302 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++)
304 digit = (AliPHOSDigit*)fGetter->Digits()->At( Digits[iDigit] );
305 eDigit = Energies[iDigit];
306 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
307 fGeom->RelPosInModule(relid, iy, ix);
309 Float_t dx = ix - xx;
310 Float_t dy = iy - yy;
311 Float_t du = dx*cosfi + dy*sinfi;
312 eu3 += eDigit*du*du*du;
315 Float_t c = E*eu3*eu3/(euu*euu*euu)/2.;
317 if(eu3<0) sign = -1.;
318 Float_t r = 1.+ c + sign*TMath::Sqrt(2.*c + c*c);
320 if(TMath::Abs(r-1.)>0.1) u = eu3/euu*(r+1.)/(r-1.);
321 if(TMath::Abs(r-1.)<0.1) u = TMath::Sqrt(sqr/E/r)*(1.+r);
323 Float_t e2c = E/(1.+r);
325 Float_t u1 = -u/(1.+r);
328 Float_t x1c = xx+u1*cosfi;
329 Float_t y1c = yy+u1*sinfi;
330 Float_t x2c = xx+u2*cosfi;
331 Float_t y2c = yy+u2*sinfi;
333 // printf("e1c -> %f\n",e1c);
334 // printf("x1c -> %f\n",x1c);
335 // printf("y1c -> %f\n",y1c);
336 // printf("e2c -> %f\n",e2c);
337 // printf("x2c -> %f\n",x2c);
338 // printf("y2c -> %f\n",y2c);
350 void AliPHOSEvalRecPoint::TwoGam(Float_t* gamma1, Float_t* gamma2)
352 //Fitting algorithm for the two very closed gammas
353 //that merged into the cluster with one maximum.
354 // Starting points gamma1 and gamma2 must be provided before ( see InitTwoGam)
355 //Set chisquare of the fit.
358 Float_t st0 = GetReconstructionManager()->TwoGamInitialStep();
359 Float_t chmin = GetReconstructionManager()->TwoGamChisqMin();
360 Float_t emin = GetReconstructionManager()->TwoGamEmin();
361 Float_t stpmin = GetReconstructionManager()->TwoGamStepMin();
362 Int_t Niter = GetReconstructionManager()->TwoGamNumOfIterations(); // Number of iterations.
364 Float_t chisq = 100.; //Initial chisquare.
366 Int_t nadc = GetMultiplicity();
370 Int_t dof = nadc - 5;
372 Float_t chstop = chmin*dof;
382 Float_t ee1=0,xx1=0,yy1=0,ee2=0,xx2=0,yy2=0,cosi;
384 Float_t e1c = gamma1[0];
385 Float_t y1c = gamma1[1];
386 Float_t x1c = gamma1[2];
388 Float_t e2c = gamma2[0];
389 Float_t y2c = gamma2[1];
390 Float_t x2c = gamma2[2];
392 Float_t E = GetEnergy();
395 AliPHOSDigit * digit ;
396 Int_t nDigits = GetMultiplicity();
397 Int_t * Digits = GetDigitsList();
398 Float_t * Energies = GetEnergiesList();
404 AliPHOSGetter* fGetter = AliPHOSGetter::GetInstance();
405 const AliPHOSGeometry* fGeom = fGetter->PHOSGeometry();
407 for(Int_t Iter=0; Iter<Niter; Iter++)
416 for(Int_t iDigit = 0 ; iDigit < nDigits ; iDigit ++)
418 digit = (AliPHOSDigit*)fGetter->Digits()->At( Digits[iDigit] );
419 eDigit = Energies[iDigit];
420 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
421 fGeom->RelPosInModule(relid, iy, ix);
426 Float_t dx1 = x1c - ix;
427 Float_t dy1 = y1c - iy;
428 // cout<<" Mult "<<nDigits<<" dx1 "<<dx1<<" dy1 "<<dy1<<endl;
429 // AG(e1c,dx1,dy1,a1,gx1,gy1);
430 GetReconstructionManager()->AG(e1c,dx1,dy1,a1,gx1,gy1);
432 Float_t dx2 = x2c - ix;
433 Float_t dy2 = y2c - iy;
434 // cout<<" "<<" dx2 "<<dx2<<" dy2 "<<dy2<<endl;
435 // AG(e2c,dx2,dy2,a2,gx2,gy2);
436 GetReconstructionManager()->AG(e2c,dx2,dy2,a2,gx2,gy2);
439 // Float_t D = Const*A*(1. - A/E);
442 // // D = 9.; // ????
444 // Float_t da = A - eDigit;
445 // chisqc += da*da/D;
446 // Float_t dd = da/D;
447 // dd = dd*(2.-dd*Const*(1.-2.*A/E));
450 chisqc += GetReconstructionManager()->TwoGamChi2(A,eDigit,E,dd);
456 grec += (a1/e1c - a2/e2c)*E*dd;
460 Float_t grc = TMath::Sqrt(grx1c*grx1c + gry1c*gry1c + grx2c*grx2c + gry2c*gry2c);
461 if(grc<1.e-10) grc=1.e-10;
463 Float_t sc = 1. + chisqc/ch;
468 cosi = (grx1*grx1c + gry1*gry1c + grx2*grx2c + gry2*gry2c + gre*grec)/gr/grc;
469 st = st*sc/(1.4 - cosi);
491 Float_t step = st*gr;
494 cout<<" Iteration "<<Iter<<" dof "<<dof<<" chisq/dof "<<ch/dof<<" chstop/dof "<<chstop/dof<<" step " <<step<<" stpmin "<<stpmin<<endl;
500 Float_t de = st*gre*E;
504 if( (e1c>emin) && (e2c>emin) )
520 // if(ch>chisq*(nadc-2)-delch)
532 Float_t x1_new = yy1;
533 Float_t z1_new = xx1;
534 Float_t e1_new = ee1;
535 Float_t x2_new = yy2;
536 Float_t z2_new = xx2;
537 Float_t e2_new = ee2;
539 cout<<" (x,z,E)[1 fit] = ("<<x1_new<<","<<z1_new<<","<<e1_new<<")"<<endl;
540 cout<<" (x,z,E)[2 fit] = ("<<x2_new<<","<<z2_new<<","<<e2_new<<")"<<endl;
546 void AliPHOSEvalRecPoint::UnfoldTwoMergedPoints(Float_t* gamma1, Float_t* gamma2)
548 //Unfold TWO merged rec. points in the case when cluster has only one maximum,
549 //but it's fitting to the one gamma shower is too bad.
550 // Use TwoGam() to estimate the positions and energies of merged points.
556 Int_t* Digits = GetDigitsList();
557 Int_t nDigits = GetMultiplicity();
558 Float_t* Energies = GetEnergiesList();
559 Float_t* eFit = new Float_t[nDigits];
561 AliPHOSGetter* fGetter = AliPHOSGetter::GetInstance();
562 const AliPHOSGeometry* fGeom = fGetter->PHOSGeometry();
564 for(Int_t iDigit=0; iDigit<nDigits; iDigit++)
566 AliPHOSDigit* digit = (AliPHOSDigit*)fGetter->Digits()->At( Digits[iDigit] );
568 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
570 fGeom->RelPosInModule(relid, x, z);
573 for(Int_t iMax=0; iMax<Nmax; iMax++)
580 Float_t eMax = gamma[0];
581 Float_t xMax = gamma[1];
582 Float_t zMax = gamma[2];
584 Float_t dx = xMax - x;
585 Float_t dz = zMax - z;
587 GetReconstructionManager()->AG(eMax,dz,dx,Amp,gx,gy);
595 for( Int_t iMax=0; iMax<Nmax; iMax++)
602 Float_t eMax = gamma[0];
603 Float_t xMax = gamma[1];
604 Float_t zMax = gamma[2];
606 AliPHOSEvalRecPoint* newRP = new AliPHOSEvalRecPoint(IsCPV(),this);
607 TVector3 newpos(xMax,0,zMax);
608 newRP->SetLocalPosition(newpos);
610 for( Int_t iDigit = 0 ; iDigit < nDigits ; iDigit ++)
612 AliPHOSDigit* digit = (AliPHOSDigit*)fGetter->Digits()->At( Digits[iDigit] );
613 Float_t eDigit = Energies[iDigit];
615 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
617 fGeom->RelPosInModule(relid, ix, iz);
619 Float_t dx = xMax - ix;
620 Float_t dz = zMax - iz;
621 Float_t single_shower_gain,gxMax,gyMax;
622 GetReconstructionManager()->AG(eMax,dz,dx,single_shower_gain,gxMax,gyMax);
623 Float_t total_gain = eFit[iDigit];
624 Float_t ratio = single_shower_gain/total_gain;
625 // cout<<" ratio -> "<<ratio<<endl;
626 eDigit = eDigit*ratio;
627 newRP->AddDigit(*digit,eDigit);
630 cout<<"======= Split: daughter rec point "<<iMax<<" ================="<<endl;
642 void AliPHOSEvalRecPoint::EvaluatePosition()
644 // One gamma fit algorithm.
645 // Set chisq/dof of the fit.
646 // gamma1[0], gamma1[1], gamma1[2] are Energy,x,z of the reconstructed gamma, respectively.
648 Int_t nDigits = GetMultiplicity();
652 Int_t Niter = GetReconstructionManager()->OneGamNumOfIterations(); // number of iterations
653 Float_t St0 = GetReconstructionManager()->OneGamInitialStep(); // initial step
654 // const Float_t Stpmin = 0.005;
655 Float_t Stpmin = GetReconstructionManager()->OneGamStepMin();
656 Float_t chmin = GetReconstructionManager()->OneGamChisqMin();
664 GetLocalPosition(locpos);
665 Float_t E = GetEnergy();
666 Float_t xc = locpos.Z();
667 Float_t yc = locpos.X();
668 Float_t dx,dy,gx,gy,grxc,gryc;
670 Float_t chisq = 1.e+20;
674 Int_t dof = GetMultiplicity() - 2;
677 Float_t chstop = chmin*dof;
678 Float_t cosi,x1=0,y1=0;
681 AliPHOSGetter* fGetter = AliPHOSGetter::GetInstance();
682 const AliPHOSGeometry* fGeom = fGetter->PHOSGeometry();
684 for(Int_t Iter=0; Iter<Niter; Iter++)
690 for(Int_t iDigit = 0 ; iDigit < nDigits ; iDigit ++)
693 Float_t* Energies = GetEnergiesList();
694 Int_t* Digits = GetDigitsList();
695 digit = (AliPHOSDigit*)fGetter->Digits()->At( Digits[iDigit] );
696 eDigit = Energies[iDigit];
697 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
698 fGeom->RelPosInModule(relid, iy, ix);
707 GetReconstructionManager()->AG(E,dx,dy,A,gx,gy);
709 cout<<" (ix iy xc yc dx dy) "<<ix<<" "<<iy<<" "<<xc<<" "<<yc<<" "<<dx<<" "<<dy<<endl;
710 Float_t chi2dg = GetReconstructionManager()->OneGamChi2(A,eDigit,E,dd);
712 // Exclude digit with too large chisquare.
713 if(chi2dg > 10) { continue; }
721 Float_t grc = TMath::Sqrt(grxc*grxc + gryc*gryc);
722 if(grc<1.e-10) grc=1.e-10;
723 Float_t sc = 1. + chisqc/chisq;
724 cout<<" chisq: "<<chisq<<endl;
728 cosi = (grx*grxc + gry*gryc)/gr/grc;
729 st = st*sc/(1.4 - cosi);
741 Float_t step = st*gr;
744 cout<<" Iteration "<<Iter<<" dof "<<dof<<" chisq/dof "<<chisq/dof<<" chstop/dof "<<chstop/dof<<" step " <<step<<" Stpmin "<<Stpmin<<endl;
769 TVector3 newpos(gamma1[1],0,gamma1[2]);
770 //SetLocalPosition(newpos);
776 // RP->GetLocalPosition(pos);
777 // cout<<" (x,z)[old] = ("<<x_old<<","<<z_old<<")"<<endl;
778 // cout<<" (x,z)[new] = ("<<pos.X()<<","<<pos.Z()<<")"<<endl;
781 // cout<<" gamma[1]: "<<gamma1[1]<<" gamma1[2]: "<<gamma1[2]<<endl;
788 Bool_t AliPHOSEvalRecPoint::KillWeakPoint()
790 // Remove this point from further procession
791 // if it's energy is too small.
793 Float_t Thr0 = GetReconstructionManager()->KillGamMinEnergy();
795 if(GetEnergy()<Thr0) {
796 cout<<"+++++++ Killing this rec point ++++++++++"<<endl;
797 RemoveFromWorkingPool(this);
805 void AliPHOSEvalRecPoint::SplitMergedShowers()
807 // Split two very close rec. points.
809 if(GetMultiplicity()<2)
815 InitTwoGam(gamma1,gamma2); // start points for fit
816 TwoGam(gamma1,gamma2); // evaluate the positions and energies
817 UnfoldTwoMergedPoints(gamma1,gamma2); // make the job
822 void AliPHOSEvalRecPoint::MergeClosePoint()
825 AliPHOSGetter* fGetter = AliPHOSGetter::GetInstance();
826 // AliPHOSDigitizer* digitizer = fGetter->Digitizer();
827 // Float_t fPedestal = digitizer->GetPedestal(); // YVK 30.09.2001
828 // Float_t fSlope = digitizer->GetSlope();
830 for (Int_t i=0;i<InWorkingPool(); i++)
832 TObject* obj = GetFromWorkingPool(i);
833 if(!((TObject*)this)->IsEqual(obj))
835 AliPHOSRecPoint* rp = (AliPHOSRecPoint*)obj;
836 if(GetPHOSMod() == rp->GetPHOSMod())
840 cout<<"+++++++ Merging point 1: ++++++"<<endl;
842 cout<<"+++++++ and point 2: ++++++++++"<<endl;
843 ((AliPHOSEvalRecPoint*)rp)->Print("");
845 //merge two rec. points
848 this->GetLocalPosition(lpos1);
849 rp->GetLocalPosition(lpos2);
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();
858 for(Int_t iDigit=0; iDigit<rp->GetDigitsMultiplicity(); iDigit++)
860 AliPHOSDigit* digit = (AliPHOSDigit*)fGetter->Digits()->At(Digits[iDigit]);
861 Float_t eDigit = Energies[iDigit];
862 this->AddDigit(*digit,eDigit);
865 TVector3 newpos(new_x,0,new_z);
868 RemoveFromWorkingPool(rp);
871 cout<<"++++++ Resulting point: ++++++++"<<endl;
882 Int_t AliPHOSEvalRecPoint::UnfoldLocalMaxima()
884 // Make unfolding in the reconstruction point with several local maxima.
885 // Return the number of local maxima.
887 // if multiplicity less then 2 - nothing to unfold
888 if(GetMultiplicity()<2) return 1;
890 AliPHOSDigit * maxAt[1000];
891 Float_t maxAtEnergy[1000];
892 Float_t LocMaxCut, LogWeight;
897 // AliPHOSClusterizer* clusterizer = fGetter->Clusterizer("PHOSclu-v1");
898 AliPHOSClusterizer* clusterizer = GetClusterizer();
900 cout<<" Cannot get clusterizer. Exit."<<endl;
905 LocMaxCut = clusterizer->GetEmcLocalMaxCut();
906 LogWeight = clusterizer->GetEmcLogWeight();
909 LocMaxCut = clusterizer->GetCpvLocalMaxCut();
910 LogWeight = clusterizer->GetCpvLogWeight();
913 AliPHOSGetter* fGetter = AliPHOSGetter::GetInstance();
914 const AliPHOSGeometry* fGeom = fGetter->PHOSGeometry();
915 TClonesArray* digits = fGetter->Digits();
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;
921 Int_t* Digits = GetDigitsList();
922 Int_t nDigits = GetMultiplicity();
923 Float_t* Energies = GetEnergiesList();
924 Float_t* eFit = new Float_t[nDigits];
926 Int_t iDigit; //loop variable
928 for(iDigit=0; iDigit<nDigits; iDigit++)
933 for(iDigit=0; iDigit<nDigits; iDigit++)
936 AliPHOSDigit* digit = (AliPHOSDigit*)fGetter->Digits()->At( Digits[iDigit] );
937 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
939 fGeom->RelPosInModule(relid, x, z);
941 for(Int_t iMax=0; iMax<Nmax; iMax++)
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;
950 GetReconstructionManager()->AG(eMax,dz,dx,Amp,gx,gy);
957 for(Int_t iMax=0; iMax<Nmax; iMax++)
959 AliPHOSDigit* digitMax = maxAt[iMax];
960 fGeom->AbsToRelNumbering(digitMax->GetId(), relid) ;
961 fGeom->RelPosInModule(relid, xMax, zMax);
962 Float_t eMax = maxAtEnergy[iMax];
964 AliPHOSEvalRecPoint* newRP = new AliPHOSEvalRecPoint(IsCPV(),this);
965 newRP->AddDigit(*digitMax,maxAtEnergy[iMax]);
967 //Neighbous ( matrix 3x3 around the local maximum)
968 for(Int_t iDigit=0; iDigit<nDigits; iDigit++)
970 AliPHOSDigit* digit = (AliPHOSDigit*)fGetter->Digits()->At( Digits[iDigit] );
971 Float_t eDigit = Energies[iDigit];
972 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
974 fGeom->RelPosInModule(relid, ix, iz);
976 if(AreNeighbours(digitMax,digit))
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 cout<<" ratio -> "<<ratio<<endl;
985 eDigit = eDigit*ratio;
986 newRP->AddDigit(*digit,eDigit);
990 newRP->EvalLocalPosition(LogWeight,digits);
991 cout<<"======= Unfold: daughter rec point "<<iMax<<" ================="<<endl;
995 // RemoveFromWorkingPool(this);
1004 void AliPHOSEvalRecPoint::PrintPoint(Option_t* opt)
1006 AliPHOSCpvRecPoint::Print(opt);
1009 GetLocalPosition(lpos);
1011 cout<<" Chi2/dof = "<<Chi2Dof()<<endl;
1012 cout<<" Local (x,z) = ("<<lpos.X()<<","<<lpos.Z()<<") in module "<<GetPHOSMod()<<endl;
1013 // if(GetReconstructionManager())
1014 // cout<<" Distance of merger = "<<GetReconstructionManager()->MergeGammasMinDistanceCut()<<endl;
1017 AliPHOSRecManager* AliPHOSEvalRecPoint::GetReconstructionManager() const
1020 TFolder* aliceF = (TFolder*)gROOT->FindObjectAny("YSAlice");
1021 TFolder* wPoolF = (TFolder*)aliceF->FindObject("WhiteBoard/RecPoints/PHOS/SmP");
1022 AliPHOSRecManager* recMng = (AliPHOSRecManager*)wPoolF->FindObject("AliPHOSRecManager");
1024 cout<<" Couldn't find Reconstruction Manager. Exit."<<endl;
1032 AliPHOSRecPoint* AliPHOSEvalRecPoint::Parent()
1034 if(fParent<0) return NULL;
1036 return (AliPHOSRecPoint*)GetFromWorkingPool(fParent);
1040 Float_t AliPHOSEvalRecPoint::Chi2Dof() const
1045 const TObject* AliPHOSEvalRecPoint::GetWorkingPool()
1048 TFolder* aliceF = (TFolder*)gROOT->FindObjectAny("YSAlice");
1049 TFolder* wPoolF = (TFolder*)aliceF->FindObject("WhiteBoard/RecPoints/PHOS/SmP");
1050 TObject* wPool = wPoolF->FindObject("SmartPoints");
1052 cout<<" Couldn't find Working Pool. Exit."<<endl;
1059 void AliPHOSEvalRecPoint::AddToWorkingPool(TObject* obj)
1061 ((TObjArray*)GetWorkingPool())->Add(obj);
1064 TObject* AliPHOSEvalRecPoint::GetFromWorkingPool(Int_t i)
1066 // return fWorkingPool->At(i);
1067 return ((TObjArray*)GetWorkingPool())->At(i);
1070 Int_t AliPHOSEvalRecPoint::InWorkingPool()
1072 return ((TObjArray*)GetWorkingPool())->GetEntries();
1075 void AliPHOSEvalRecPoint::RemoveFromWorkingPool(TObject* obj)
1077 ((TObjArray*)GetWorkingPool())->Remove(obj);
1078 ((TObjArray*)GetWorkingPool())->Compress();
1081 void AliPHOSEvalRecPoint::PrintWorkingPool()
1083 ((TObjArray*)GetWorkingPool())->Print("");