1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
19 //*-- Author: Boris Polichtchouk, IHEP
20 //////////////////////////////////////////////////////////////////////////////
21 // Reconstructed point operations for the Clusterization class for IHEP reconstruction.
22 // Performs clusterization (collects neighbouring active cells)
23 // It differs from AliPHOSClusterizerv1 in neighbour definition only
25 // ---- ROOT system ---
26 #include "TDirectory.h"
32 // --- AliRoot header files ---
33 #include "AliConfig.h"
34 #include "AliPHOSDigit.h"
35 #include "AliPHOSClusterizer.h"
36 #include "AliPHOSEvalRecPoint.h"
38 #include "AliPHOSLoader.h"
39 #include "AliPHOSRecCpvManager.h"
40 #include "AliPHOSRecEmcManager.h"
41 #include "AliPHOSDigitizer.h"
42 #include "AliPHOSGeometry.h"
44 // --- Standard library ---
47 ClassImp(AliPHOSEvalRecPoint)
49 AliPHOSEvalRecPoint::AliPHOSEvalRecPoint(): fEventFolderName(AliConfig::GetDefaultEventFolderName())
58 AliPHOSEvalRecPoint::AliPHOSEvalRecPoint(Bool_t cpv, AliPHOSEvalRecPoint* parent) : AliPHOSCpvRecPoint()
65 TObjArray* wPool = (TObjArray*)GetWorkingPool();
67 Fatal("ctor", "Couldn't find working pool") ;
70 fParent = wPool->IndexOf((TObject*)parent);
71 fChi2Dof = parent->Chi2Dof();
82 // Add itself to working pool
83 AddToWorkingPool((TObject*)this);
87 AliPHOSEvalRecPoint::AliPHOSEvalRecPoint(Int_t i, Bool_t cpv) : fEventFolderName(AliConfig::GetDefaultEventFolderName())
93 AliPHOSEmcRecPoint* rp=0;
95 AliPHOSLoader* fLoader = AliPHOSLoader::GetPHOSLoader(fEventFolderName);
98 rp = (AliPHOSCpvRecPoint*)fLoader->CpvRecPoints()->At(i);
103 rp = (AliPHOSEmcRecPoint*)fLoader->EmcRecPoints()->At(i);
108 Int_t* digits = rp->GetDigitsList();
109 Float_t* energies = rp->GetEnergiesList();
110 Int_t nDigits = rp->GetMultiplicity();
112 for(Int_t iDigit=0; iDigit<nDigits; iDigit++) {
113 AliPHOSDigit* digit = (AliPHOSDigit*)fLoader->Digits()->At( digits[iDigit] );
114 Float_t eDigit = energies[iDigit];
115 this->AddDigit(*digit,eDigit);
119 rp->GetLocalPosition(locpos);
122 // Add itself to working pool
123 AddToWorkingPool((TObject*)this);
127 AliPHOSClusterizer* AliPHOSEvalRecPoint::GetClusterizer()
129 // returns clusterizer task
130 TFolder* aliceF = (TFolder*)gROOT->FindObjectAny("YSAlice");
131 TFolder* wPoolF = (TFolder*)aliceF->FindObject("WhiteBoard/RecPoints/PHOS/SmP");
132 AliPHOSClusterizer* clu = (AliPHOSClusterizer*)wPoolF->FindObject("PHOS:clu-v1");
134 Fatal("GetClusterizer", "Couldn't find Clusterizer") ;
140 Bool_t AliPHOSEvalRecPoint::TooClose(AliPHOSRecPoint* pt) const
142 // check if a rec.point pt is too close to this one
145 pt->GetLocalPosition(herPos);
146 this->GetLocalPosition(myPos);
147 Float_t dx = herPos.X() - myPos.X();
148 Float_t dz = herPos.Z() - myPos.Z();
149 Float_t dr = TMath::Sqrt(dx*dx + dz*dz);
151 if(dr<GetReconstructionManager()->MergeGammasMinDistanceCut())
158 Bool_t AliPHOSEvalRecPoint::NeedToSplit() const
160 // rec.point needs to be split
164 void AliPHOSEvalRecPoint::DeleteParent()
170 void AliPHOSEvalRecPoint::UpdateWorkingPool()
172 // update pool of rec.points
173 Int_t i; //loop variable
175 for(i=0; i<this->InWorkingPool(); i++) {
176 AliPHOSEvalRecPoint* parent = (AliPHOSEvalRecPoint*)GetFromWorkingPool(i);
178 Int_t nChild = parent->HasChild(children);
179 for(Int_t iChild=0; iChild<nChild; iChild++)
181 ((AliPHOSEvalRecPoint*)children.At(iChild))->DeleteParent();
185 RemoveFromWorkingPool(parent);
191 for(i=0; i<this->InWorkingPool(); i++) {
192 AliPHOSEvalRecPoint* weak = (AliPHOSEvalRecPoint*)GetFromWorkingPool(i);
193 if (weak->KillWeakPoint()) delete weak;
196 for(i=0; i<this->InWorkingPool(); i++) {
197 AliPHOSEvalRecPoint* close = (AliPHOSEvalRecPoint*)GetFromWorkingPool(i);
198 close->MergeClosePoint();
201 for(i=0; i<this->InWorkingPool(); i++)
202 ((AliPHOSEvalRecPoint*)AliPHOSEvalRecPoint::GetFromWorkingPool(i))->SetIndexInList(i);
206 Int_t AliPHOSEvalRecPoint::HasChild(TObjArray& children)
208 // returns the number of children
209 for(Int_t iChild=0; iChild<InWorkingPool(); iChild++)
211 AliPHOSEvalRecPoint* child = (AliPHOSEvalRecPoint*)GetFromWorkingPool(iChild);
212 TObject* parent = (TObject*)child->Parent();
213 TObject* me = (TObject*)this;
214 if(parent==me) children.Add(child);
217 return children.GetEntries();
220 void AliPHOSEvalRecPoint::Init()
223 AliPHOSClusterizer* clusterizer = GetClusterizer();
225 Fatal("Init", "Cannot get clusterizer") ;
228 TClonesArray* digits = AliPHOSLoader::GetPHOSLoader(fEventFolderName)->Digits();
232 logWeight = clusterizer->GetEmcLogWeight();
235 logWeight = clusterizer->GetCpvLogWeight();
238 EvalLocalPosition(logWeight,digits); // evaluate initial position
242 void AliPHOSEvalRecPoint::MakeJob()
244 // Reconstruction algoritm implementation.
246 AliPHOSRecManager* recMng = GetReconstructionManager();
253 Int_t nChild = HasChild(children);
257 if(Chi2Dof()>recMng->OneGamChisqCut()) SplitMergedShowers();
260 for(Int_t iChild=0; iChild<nChild; iChild++) {
261 AliPHOSEvalRecPoint* child = (AliPHOSEvalRecPoint*)children.At(iChild);
262 child->EvaluatePosition();
264 if(child->Chi2Dof()>recMng->OneGamChisqCut())
265 child->SplitMergedShowers();
270 void AliPHOSEvalRecPoint::InitTwoGam(Float_t* gamma1, Float_t* gamma2)
272 //Compute start values for two gamma fit algorithm.
273 // gamma1[0], gamma1[1], gamma1[2] are Energy,x,z of the reconstructed "gammas".
276 TVector3 lpos; // start point choosing.
277 GetLocalPosition(lpos);
278 Float_t xx = lpos.Z();
279 Float_t yy = lpos.X();
280 Float_t e = GetEnergy();
282 Info("InitTwoGam", "(x,z,e)[old] = (%f, %f, %f)", yy, xx, e) ;
289 AliPHOSDigit * digit ;
290 Int_t nDigits = GetMultiplicity();
291 Int_t * digits = GetDigitsList();
292 Float_t * energies = GetEnergiesList();
304 AliPHOSLoader* fLoader = AliPHOSLoader::GetPHOSLoader(fEventFolderName);
305 const AliPHOSGeometry* fGeom = fLoader->PHOSGeometry();
307 Int_t iDigit; //loop variable
309 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++)
311 digit = (AliPHOSDigit*)fLoader->Digits()->At( digits[iDigit] );
312 eDigit = energies[iDigit];
313 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
314 fGeom->RelPosInModule(relid, iy, ix);
316 Float_t dx = ix - xx;
317 Float_t dy = iy - yy;
324 sqr = TMath::Sqrt(4.*exy*exy + (exx-eyy)*(exx-eyy));
325 Float_t euu = (exx+eyy+sqr)/2.;
326 if(sqr>1.e-10) cos2fi = (exx - eyy)/sqr;
327 Float_t cosfi = TMath::Sqrt((1.+cos2fi)/2.);
328 Float_t sinfi = TMath::Sqrt((1.-cos2fi)/2.);
329 if(exy<0) sinfi = -sinfi;
332 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++)
334 digit = (AliPHOSDigit*)fLoader->Digits()->At( digits[iDigit] );
335 eDigit = energies[iDigit];
336 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
337 fGeom->RelPosInModule(relid, iy, ix);
339 Float_t dx = ix - xx;
340 Float_t dy = iy - yy;
341 Float_t du = dx*cosfi + dy*sinfi;
342 eu3 += eDigit*du*du*du;
345 Float_t c = e*eu3*eu3/(euu*euu*euu)/2.;
347 if(eu3<0) sign = -1.;
348 Float_t r = 1.+ c + sign*TMath::Sqrt(2.*c + c*c);
350 if(TMath::Abs(r-1.)>0.1) u = eu3/euu*(r+1.)/(r-1.);
351 if(TMath::Abs(r-1.)<0.1) u = TMath::Sqrt(sqr/e/r)*(1.+r);
353 Float_t e2c = e/(1.+r);
355 Float_t u1 = -u/(1.+r);
358 Float_t x1c = xx+u1*cosfi;
359 Float_t y1c = yy+u1*sinfi;
360 Float_t x2c = xx+u2*cosfi;
361 Float_t y2c = yy+u2*sinfi;
363 // printf("e1c -> %f\n",e1c);
364 // printf("x1c -> %f\n",x1c);
365 // printf("y1c -> %f\n",y1c);
366 // printf("e2c -> %f\n",e2c);
367 // printf("x2c -> %f\n",x2c);
368 // printf("y2c -> %f\n",y2c);
380 void AliPHOSEvalRecPoint::TwoGam(Float_t* gamma1, Float_t* gamma2)
382 //Fitting algorithm for the two very closed gammas
383 //that merged into the cluster with one maximum.
384 // Starting points gamma1 and gamma2 must be provided before ( see InitTwoGam)
385 //Set chisquare of the fit.
388 Float_t st0 = GetReconstructionManager()->TwoGamInitialStep();
389 Float_t chmin = GetReconstructionManager()->TwoGamChisqMin();
390 Float_t emin = GetReconstructionManager()->TwoGamEmin();
391 Float_t stpmin = GetReconstructionManager()->TwoGamStepMin();
392 Int_t nIter = GetReconstructionManager()->TwoGamNumOfIterations(); // Number of iterations.
394 Float_t chisq = 100.; //Initial chisquare.
396 Int_t nadc = GetMultiplicity();
400 Int_t dof = nadc - 5;
402 Float_t chstop = chmin*dof;
412 Float_t ee1=0,xx1=0,yy1=0,ee2=0,xx2=0,yy2=0,cosi;
414 Float_t e1c = gamma1[0];
415 Float_t y1c = gamma1[1];
416 Float_t x1c = gamma1[2];
418 Float_t e2c = gamma2[0];
419 Float_t y2c = gamma2[1];
420 Float_t x2c = gamma2[2];
422 Float_t e = GetEnergy();
425 AliPHOSDigit * digit ;
426 Int_t nDigits = GetMultiplicity();
427 Int_t * digits = GetDigitsList();
428 Float_t * energies = GetEnergiesList();
434 AliPHOSLoader* fLoader = AliPHOSLoader::GetPHOSLoader(fEventFolderName);
435 const AliPHOSGeometry* fGeom = fLoader->PHOSGeometry();
437 for(Int_t iter=0; iter<nIter; iter++)
446 for(Int_t iDigit = 0 ; iDigit < nDigits ; iDigit ++)
448 digit = (AliPHOSDigit*)fLoader->Digits()->At( digits[iDigit] );
449 eDigit = energies[iDigit];
450 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
451 fGeom->RelPosInModule(relid, iy, ix);
456 Float_t dx1 = x1c - ix;
457 Float_t dy1 = y1c - iy;
458 // Info("TwoGam", "Mult %d dx1 %f dy1 %f", nDigits, dx1, dy1) ;
459 // AG(e1c,dx1,dy1,a1,gx1,gy1);
460 GetReconstructionManager()->AG(e1c,dx1,dy1,a1,gx1,gy1);
462 Float_t dx2 = x2c - ix;
463 Float_t dy2 = y2c - iy;
464 // Info("TwoGam", " dx2 %f dy2 %f", dx2, dy2) ;
465 // AG(e2c,dx2,dy2,a2,gx2,gy2);
466 GetReconstructionManager()->AG(e2c,dx2,dy2,a2,gx2,gy2);
469 // Float_t D = Const*a*(1. - a/e);
472 // // D = 9.; // ????
474 // Float_t da = a - eDigit;
475 // chisqc += da*da/D;
476 // Float_t dd = da/D;
477 // dd = dd*(2.-dd*Const*(1.-2.*a/e));
480 chisqc += GetReconstructionManager()->TwoGamChi2(a,eDigit,e,dd);
486 grec += (a1/e1c - a2/e2c)*e*dd;
490 Float_t grc = TMath::Sqrt(grx1c*grx1c + gry1c*gry1c + grx2c*grx2c + gry2c*gry2c);
491 if(grc<1.e-10) grc=1.e-10;
493 Float_t sc = 1. + chisqc/ch;
498 cosi = (grx1*grx1c + gry1*gry1c + grx2*grx2c + gry2*gry2c + gre*grec)/gr/grc;
499 st = st*sc/(1.4 - cosi);
521 Float_t step = st*gr;
523 Info("TwoGam", "Iteration %d dof %d chisq/dof %f chstop/dof %f step %d stpmin %d",
524 iter, dof, ch/dof, chstop/dof, step, stpmin) ;
530 Float_t de = st*gre*e;
534 if( (e1c>emin) && (e2c>emin) )
550 // if(ch>chisq*(nadc-2)-delch)
570 message = " (x,z,e)[1 fit] = (%f, %f, %f)\n" ;
571 message = " (x,z,e)[2 fit] = (%f, %f, %f)\n" ;
572 Info("TwoGam", message.Data(),
574 x2New, z2New, e2New) ;
580 void AliPHOSEvalRecPoint::UnfoldTwoMergedPoints(Float_t* gamma1, Float_t* gamma2)
582 //Unfold TWO merged rec. points in the case when cluster has only one maximum,
583 //but it's fitting to the one gamma shower is too bad.
584 // Use TwoGam() to estimate the positions and energies of merged points.
590 Int_t* digits = GetDigitsList();
591 Int_t nDigits = GetMultiplicity();
592 Float_t* energies = GetEnergiesList();
593 Float_t* eFit = new Float_t[nDigits];
595 AliPHOSLoader* fLoader = AliPHOSLoader::GetPHOSLoader(fEventFolderName);
596 const AliPHOSGeometry* fGeom = fLoader->PHOSGeometry();
598 for(Int_t iDigit=0; iDigit<nDigits; iDigit++)
600 AliPHOSDigit* digit = (AliPHOSDigit*)fLoader->Digits()->At( digits[iDigit] );
602 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
604 fGeom->RelPosInModule(relid, x, z);
607 for(Int_t iMax=0; iMax<nMax; iMax++)
614 Float_t eMax = gamma[0];
615 Float_t xMax = gamma[1];
616 Float_t zMax = gamma[2];
618 Float_t dx = xMax - x;
619 Float_t dz = zMax - z;
621 GetReconstructionManager()->AG(eMax,dz,dx,amp,gx,gy);
629 for( Int_t iMax=0; iMax<nMax; iMax++)
636 Float_t eMax = gamma[0];
637 Float_t xMax = gamma[1];
638 Float_t zMax = gamma[2];
640 AliPHOSEvalRecPoint* newRP = new AliPHOSEvalRecPoint(IsCPV(),this);
641 TVector3 newpos(xMax,0,zMax);
642 newRP->SetLocalPosition(newpos);
644 for( Int_t iDigit = 0 ; iDigit < nDigits ; iDigit ++)
646 AliPHOSDigit* digit = (AliPHOSDigit*)fLoader->Digits()->At( digits[iDigit] );
647 Float_t eDigit = energies[iDigit];
649 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
651 fGeom->RelPosInModule(relid, ix, iz);
653 Float_t dx = xMax - ix;
654 Float_t dz = zMax - iz;
655 Float_t singleShowerGain,gxMax,gyMax;
656 GetReconstructionManager()->AG(eMax,dz,dx,singleShowerGain,gxMax,gyMax);
657 Float_t totalGain = eFit[iDigit];
658 Float_t ratio = singleShowerGain/totalGain;
659 eDigit = eDigit*ratio;
660 newRP->AddDigit(*digit,eDigit);
662 Info("UnfoldTwoMergedPoints", "======= Split: daughter rec point %d =================", iMax) ;
673 void AliPHOSEvalRecPoint::EvaluatePosition()
675 // One gamma fit algorithm.
676 // Set chisq/dof of the fit.
677 // gamma1[0], gamma1[1], gamma1[2] are Energy,x,z of the reconstructed gamma, respectively.
679 Int_t nDigits = GetMultiplicity();
683 Int_t nIter = GetReconstructionManager()->OneGamNumOfIterations(); // number of iterations
684 Float_t st0 = GetReconstructionManager()->OneGamInitialStep(); // initial step
685 // const Float_t stpMin = 0.005;
686 Float_t stpMin = GetReconstructionManager()->OneGamStepMin();
687 Float_t chmin = GetReconstructionManager()->OneGamChisqMin();
695 GetLocalPosition(locpos);
696 Float_t e = GetEnergy();
697 Float_t xc = locpos.Z();
698 Float_t yc = locpos.X();
699 Float_t dx,dy,gx,gy,grxc,gryc;
701 Float_t chisq = 1.e+20;
705 Int_t dof = GetMultiplicity() - 2;
708 Float_t chstop = chmin*dof;
709 Float_t cosi,x1=0,y1=0;
712 AliPHOSLoader* fLoader = AliPHOSLoader::GetPHOSLoader(fEventFolderName);
713 const AliPHOSGeometry* fGeom = fLoader->PHOSGeometry();
715 for(Int_t iter=0; iter<nIter; iter++)
721 for(Int_t iDigit = 0 ; iDigit < nDigits ; iDigit ++)
724 Float_t* energies = GetEnergiesList();
725 Int_t* digits = GetDigitsList();
726 digit = (AliPHOSDigit*)fLoader->Digits()->At( digits[iDigit] );
727 eDigit = energies[iDigit];
728 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
729 fGeom->RelPosInModule(relid, iy, ix);
738 GetReconstructionManager()->AG(e,dx,dy,a,gx,gy);
740 Info("EvaluatePosition", " (ix iy xc yc dx dy) %f %f %f %f %f %f", ix, iy, xc, yc, dx, dy) ;
741 Float_t chi2dg = GetReconstructionManager()->OneGamChi2(a,eDigit,e,dd);
743 // Exclude digit with too large chisquare.
744 if(chi2dg > 10) { continue; }
752 Float_t grc = TMath::Sqrt(grxc*grxc + gryc*gryc);
753 if(grc<1.e-10) grc=1.e-10;
754 Float_t sc = 1. + chisqc/chisq;
755 Info("EvaluatePosition", " chisq: %f", chisq) ;
759 cosi = (grx*grxc + gry*gryc)/gr/grc;
760 st = st*sc/(1.4 - cosi);
772 Float_t step = st*gr;
774 Info("EvaluatePosition", " Iteration %d dof %d chisq/dof %f chstop/dof %f step %d stpMin %d",
775 iter, dof, chisq/dof, chisq/dof, chstop/dof, step, stpMin) ;
800 TVector3 newpos(gamma1[1],0,gamma1[2]);
801 //SetLocalPosition(newpos);
807 // RP->GetLocalPosition(pos);
816 Bool_t AliPHOSEvalRecPoint::KillWeakPoint()
818 // Remove this point from further procession
819 // if it's energy is too small.
821 Float_t thr0 = GetReconstructionManager()->KillGamMinEnergy();
823 if(GetEnergy()<thr0) {
824 Info("KillWeakPoint", "+++++++ Killing this rec point ++++++++++") ;
825 RemoveFromWorkingPool(this);
833 void AliPHOSEvalRecPoint::SplitMergedShowers()
835 // Split two very close rec. points.
837 if(GetMultiplicity()<2)
843 InitTwoGam(gamma1,gamma2); // start points for fit
844 TwoGam(gamma1,gamma2); // evaluate the positions and energies
845 UnfoldTwoMergedPoints(gamma1,gamma2); // make the job
850 void AliPHOSEvalRecPoint::MergeClosePoint()
852 // merge rec.points if they are too close
853 AliPHOSLoader* fLoader = AliPHOSLoader::GetPHOSLoader(fEventFolderName);
854 // AliPHOSDigitizer* digitizer = fLoader->Digitizer();
855 // Float_t fPedestal = digitizer->GetPedestal(); // YVK 30.09.2001
856 // Float_t fSlope = digitizer->GetSlope();
858 for (Int_t i=0;i<InWorkingPool(); i++)
860 TObject* obj = GetFromWorkingPool(i);
861 if(!((TObject*)this)->IsEqual(obj))
863 AliPHOSRecPoint* rp = (AliPHOSRecPoint*)obj;
864 if(GetPHOSMod() == rp->GetPHOSMod())
868 Info("MergeClosePoint", "+++++++ Merging point 1: ++++++") ;
870 Info("MergeClosePoint", "+++++++ and point 2: ++++++++++") ;
871 ((AliPHOSEvalRecPoint*)rp)->Print();
873 //merge two rec. points
876 this->GetLocalPosition(lpos1);
877 rp->GetLocalPosition(lpos2);
879 Int_t* digits = rp->GetDigitsList();
880 Float_t dE = rp->GetEnergy()/(rp->GetEnergy()+this->GetEnergy());
881 Float_t newX = lpos1.X()*dE + lpos2.X()*(1.-dE);
882 Float_t newZ = lpos1.Z()*dE + lpos2.Z()*(1.-dE);
883 Float_t newE = rp->GetEnergy()+this->GetEnergy();
884 Float_t* energies = ((AliPHOSEmcRecPoint*)rp)->GetEnergiesList();
886 for(Int_t iDigit=0; iDigit<rp->GetDigitsMultiplicity(); iDigit++)
888 AliPHOSDigit* digit = (AliPHOSDigit*)fLoader->Digits()->At(digits[iDigit]);
889 Float_t eDigit = energies[iDigit];
890 this->AddDigit(*digit,eDigit);
893 TVector3 newpos(newX,0,newZ);
896 RemoveFromWorkingPool(rp);
899 Info("MergeClosePoint", "++++++ Resulting point: ++++++++") ;
910 Int_t AliPHOSEvalRecPoint::UnfoldLocalMaxima()
912 // Make unfolding in the reconstruction point with several local maxima.
913 // Return the number of local maxima.
915 // if multiplicity less then 2 - nothing to unfold
916 if(GetMultiplicity()<2) return 1;
918 AliPHOSDigit * maxAt[1000];
919 Float_t maxAtEnergy[1000];
920 Float_t locMaxCut, logWeight;
925 AliPHOSClusterizer* clusterizer = GetClusterizer();
927 Fatal("UnfoldLocalMaxima", "Cannot get clusterizer") ;
931 locMaxCut = clusterizer->GetEmcLocalMaxCut();
932 logWeight = clusterizer->GetEmcLogWeight();
935 locMaxCut = clusterizer->GetCpvLocalMaxCut();
936 logWeight = clusterizer->GetCpvLogWeight();
939 AliPHOSLoader* fLoader = AliPHOSLoader::GetPHOSLoader(fEventFolderName);
940 const AliPHOSGeometry* fGeom = fLoader->PHOSGeometry();
941 TClonesArray* digits = fLoader->Digits();
943 // if number of local maxima less then 2 - nothing to unfold
944 Int_t nMax = GetNumberOfLocalMax(maxAt,maxAtEnergy,locMaxCut,digits);
945 if(nMax<2) return nMax;
947 Int_t* digitsList = GetDigitsList();
948 Int_t nDigits = GetMultiplicity();
949 Float_t* energies = GetEnergiesList();
950 Float_t* eFit = new Float_t[nDigits];
952 Int_t iDigit; //loop variable
954 for(iDigit=0; iDigit<nDigits; iDigit++)
959 for(iDigit=0; iDigit<nDigits; iDigit++)
962 AliPHOSDigit* digit = (AliPHOSDigit*)fLoader->Digits()->At( digitsList[iDigit] );
963 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
965 fGeom->RelPosInModule(relid, x, z);
967 for(Int_t iMax=0; iMax<nMax; iMax++)
969 AliPHOSDigit* digitMax = maxAt[iMax];
970 Float_t eMax = maxAtEnergy[iMax];
971 fGeom->AbsToRelNumbering(digitMax->GetId(), relid) ;
972 fGeom->RelPosInModule(relid, xMax, zMax);
973 Float_t dx = xMax - x;
974 Float_t dz = zMax - z;
976 GetReconstructionManager()->AG(eMax,dz,dx,amp,gx,gy);
983 for(Int_t iMax=0; iMax<nMax; iMax++)
985 AliPHOSDigit* digitMax = maxAt[iMax];
986 fGeom->AbsToRelNumbering(digitMax->GetId(), relid) ;
987 fGeom->RelPosInModule(relid, xMax, zMax);
988 Float_t eMax = maxAtEnergy[iMax];
990 AliPHOSEvalRecPoint* newRP = new AliPHOSEvalRecPoint(IsCPV(),this);
991 newRP->AddDigit(*digitMax,maxAtEnergy[iMax]);
993 //Neighbous ( matrix 3x3 around the local maximum)
994 for(Int_t iDigit=0; iDigit<nDigits; iDigit++)
996 AliPHOSDigit* digit = (AliPHOSDigit*)fLoader->Digits()->At( digitsList[iDigit] );
997 Float_t eDigit = energies[iDigit];
998 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
1000 fGeom->RelPosInModule(relid, ix, iz);
1002 if(AreNeighbours(digitMax,digit))
1004 Float_t dx = xMax - ix;
1005 Float_t dz = zMax - iz;
1006 Float_t singleShowerGain,gxMax,gyMax;
1007 GetReconstructionManager()->AG(eMax,dz,dx,singleShowerGain,gxMax,gyMax);
1008 Float_t totalGain = eFit[iDigit];
1009 Float_t ratio = singleShowerGain/totalGain;
1010 Info("UnfoldLocalMaxima", " ratio -> %f", ratio) ;
1011 eDigit = eDigit*ratio;
1012 newRP->AddDigit(*digit,eDigit);
1016 newRP->EvalLocalPosition(logWeight,digits);
1017 Info("UnfoldLocalMaxima", "======= Unfold: daughter rec point %d =================", iMax) ;
1021 // RemoveFromWorkingPool(this);
1029 void AliPHOSEvalRecPoint::PrintPoint()
1031 // print rec.point to stdout
1032 AliPHOSCpvRecPoint::Print();
1035 GetLocalPosition(lpos);
1038 message = " Chi2/dof = %f" ;
1039 message += " Local (x,z) = (%f, %f) in module %d" ;
1040 Info("Print", message.Data(), Chi2Dof(), lpos.X(), lpos.Z(), GetPHOSMod()) ;
1044 AliPHOSRecManager* AliPHOSEvalRecPoint::GetReconstructionManager() const
1046 // returns reconstruction manager
1047 TFolder* aliceF = (TFolder*)gROOT->FindObjectAny("YSAlice");
1048 TFolder* wPoolF = (TFolder*)aliceF->FindObject("WhiteBoard/RecPoints/PHOS/SmP");
1049 AliPHOSRecManager* recMng = (AliPHOSRecManager*)wPoolF->FindObject("AliPHOSRecManager");
1051 Fatal("GetReconstructionManager", "Couldn't find Reconstruction Manager") ;
1058 AliPHOSRecPoint* AliPHOSEvalRecPoint::Parent()
1061 if(fParent<0) return NULL;
1063 return (AliPHOSRecPoint*)GetFromWorkingPool(fParent);
1067 Float_t AliPHOSEvalRecPoint::Chi2Dof() const
1069 // returns a chi^2 per degree of freedom
1073 const TObject* AliPHOSEvalRecPoint::GetWorkingPool()
1075 // returns a pool of rec.points
1076 TFolder* aliceF = (TFolder*)gROOT->FindObjectAny("YSAlice");
1077 TFolder* wPoolF = (TFolder*)aliceF->FindObject("WhiteBoard/RecPoints/PHOS/SmP");
1078 TObject* wPool = wPoolF->FindObject("SmartPoints");
1080 Fatal("GetWorkingPool", "Couldn't find Working Pool") ;
1086 void AliPHOSEvalRecPoint::AddToWorkingPool(TObject* obj)
1088 // add a rec.point to a pool
1089 ((TObjArray*)GetWorkingPool())->Add(obj);
1092 TObject* AliPHOSEvalRecPoint::GetFromWorkingPool(Int_t i)
1094 // return a rec.point from a pool at an index i
1095 // return fWorkingPool->At(i);
1096 return ((TObjArray*)GetWorkingPool())->At(i);
1099 Int_t AliPHOSEvalRecPoint::InWorkingPool()
1101 // return the number of rec.points in a pool
1102 return ((TObjArray*)GetWorkingPool())->GetEntries();
1105 void AliPHOSEvalRecPoint::RemoveFromWorkingPool(TObject* obj)
1107 // remove a rec.point from a pool
1108 ((TObjArray*)GetWorkingPool())->Remove(obj);
1109 ((TObjArray*)GetWorkingPool())->Compress();
1112 void AliPHOSEvalRecPoint::PrintWorkingPool()
1114 // print pool of rec.points to stdout
1115 ((TObjArray*)GetWorkingPool())->Print();