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