]> 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
Some additional changes related to the previous changes. AliL3Transform
[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 ---
17#include <iostream.h>
18
19
20ClassImp(AliPHOSEvalRecPoint)
21
22 AliPHOSEvalRecPoint::AliPHOSEvalRecPoint(): AliPHOSCpvRecPoint()
23{
24 fParent=-333;
25 fChi2Dof=-111;
26 fIsCpv = kTRUE;
27 fIsEmc = kFALSE;
28}
29
30AliPHOSEvalRecPoint::AliPHOSEvalRecPoint(Bool_t cpv, AliPHOSEvalRecPoint* parent) : AliPHOSCpvRecPoint()
31{
32
33 fParent=-333;
34 fChi2Dof=-111;
35
36// fParent=parent;
37 TObjArray* wPool = (TObjArray*)GetWorkingPool();
38 if(!wPool) {
39 cout<<" Couldn't find working pool. Exit."<<endl;
40 exit(1);
41 }
42
43 fParent = wPool->IndexOf((TObject*)parent);
44 fChi2Dof = parent->Chi2Dof();
45
46 if(cpv) {
47 fIsEmc = kFALSE;
48 fIsCpv = kTRUE;
49 }
50 else {
51 fIsEmc = kTRUE;
52 fIsCpv = kFALSE;
53 }
54
55 // Add itself to working pool
56 AddToWorkingPool((TObject*)this);
57
58}
59
60AliPHOSEvalRecPoint::AliPHOSEvalRecPoint(Int_t i, Bool_t cpv) : AliPHOSCpvRecPoint()
61{
62
63 fChi2Dof=-111;
64 fParent=-333;
65
66 AliPHOSEmcRecPoint* rp=0;
67
68 AliPHOSGetter* fGetter = AliPHOSGetter::GetInstance();
69
70 if(cpv) {
71 rp = (AliPHOSCpvRecPoint*)fGetter->CpvRecPoints()->At(i);
72 fIsEmc = kFALSE;
73 fIsCpv = kTRUE;
74 }
75 else {
76 rp = (AliPHOSEmcRecPoint*)fGetter->EmcRecPoints()->At(i);
77 fIsEmc = kTRUE;
78 fIsCpv = kFALSE;
79 }
80
81 Int_t* Digits = rp->GetDigitsList();
82 Float_t* Energies = rp->GetEnergiesList();
83 Int_t nDigits = rp->GetMultiplicity();
84
85 for(Int_t iDigit=0; iDigit<nDigits; iDigit++) {
86 AliPHOSDigit* digit = (AliPHOSDigit*)fGetter->Digits()->At( Digits[iDigit] );
87 Float_t eDigit = Energies[iDigit];
88 this->AddDigit(*digit,eDigit);
89 }
90
91 TVector3 locpos;
92 rp->GetLocalPosition(locpos);
93 fLocPos = locpos;
94
95 // Add itself to working pool
96 AddToWorkingPool((TObject*)this);
97
98}
99
100AliPHOSClusterizer* AliPHOSEvalRecPoint::GetClusterizer()
101{
102 TFolder* aliceF = (TFolder*)gROOT->FindObjectAny("YSAlice");
103 TFolder* wPoolF = (TFolder*)aliceF->FindObject("WhiteBoard/RecPoints/PHOS/SmP");
104 AliPHOSClusterizer* clu = (AliPHOSClusterizer*)wPoolF->FindObject("PHOS:clu-v1");
105 if(!clu) {
106 cout<<" Couldn't find Clusterizer. Exit."<<endl;
107 exit(1);
108 }
109
110 return clu;
111}
112
113Bool_t AliPHOSEvalRecPoint::TooClose(AliPHOSRecPoint* pt) const
114{
115 TVector3 her_pos;
116 TVector3 my_pos;
117 pt->GetLocalPosition(her_pos);
118 this->GetLocalPosition(my_pos);
119 Float_t dx = her_pos.X() - my_pos.X();
120 Float_t dz = her_pos.Z() - my_pos.Z();
121 Float_t dr = TMath::Sqrt(dx*dx + dz*dz);
122
123 if(dr<GetReconstructionManager()->MergeGammasMinDistanceCut())
124 return kTRUE;
125 else
126 return kFALSE;
127
128}
129
130Bool_t AliPHOSEvalRecPoint::NeedToSplit() const
131{
132 return kFALSE;
133}
134
135void AliPHOSEvalRecPoint::DeleteParent()
136{
137 fParent=-333;
138}
139
140void AliPHOSEvalRecPoint::UpdateWorkingPool()
141{
c9eeb00e 142
143 Int_t i; //loop variable
cbd576a6 144
c9eeb00e 145 for(i=0; i<this->InWorkingPool(); i++) {
cbd576a6 146 AliPHOSEvalRecPoint* parent = (AliPHOSEvalRecPoint*)GetFromWorkingPool(i);
147 TObjArray children;
148 Int_t nChild = parent->HasChild(children);
149 for(Int_t iChild=0; iChild<nChild; iChild++)
150 {
151 ((AliPHOSEvalRecPoint*)children.At(iChild))->DeleteParent();
152 }
153
154 if(nChild) {
155 RemoveFromWorkingPool(parent);
156 delete parent;
157 }
158
159 }
160
c9eeb00e 161 for(i=0; i<this->InWorkingPool(); i++) {
cbd576a6 162 AliPHOSEvalRecPoint* weak = (AliPHOSEvalRecPoint*)GetFromWorkingPool(i);
163 if (weak->KillWeakPoint()) delete weak;
164 }
165
c9eeb00e 166 for(i=0; i<this->InWorkingPool(); i++) {
cbd576a6 167 AliPHOSEvalRecPoint* close = (AliPHOSEvalRecPoint*)GetFromWorkingPool(i);
168 close->MergeClosePoint();
169 }
170
c9eeb00e 171 for(i=0; i<this->InWorkingPool(); i++)
cbd576a6 172 ((AliPHOSEvalRecPoint*)AliPHOSEvalRecPoint::GetFromWorkingPool(i))->SetIndexInList(i);
173
174}
175
176Int_t AliPHOSEvalRecPoint::HasChild(TObjArray& children)
177{
178 for(Int_t iChild=0; iChild<InWorkingPool(); iChild++)
179 {
180 AliPHOSEvalRecPoint* child = (AliPHOSEvalRecPoint*)GetFromWorkingPool(iChild);
181 TObject* parent = (TObject*)child->Parent();
182 TObject* me = (TObject*)this;
183 if(parent==me) children.Add(child);
184 }
185
186 return children.GetEntries();
187}
188
189void AliPHOSEvalRecPoint::Init()
190{
191
192 AliPHOSClusterizer* clusterizer = GetClusterizer();
193 if(!clusterizer) {
194 cout<<" Cannot get clusterizer. Exit."<<endl;
195 exit(1);
196 }
197
198 TClonesArray* digits = AliPHOSGetter::GetInstance()->Digits();
199 Float_t LogWeight=0;
200
201 if(this->IsEmc()) {
202 LogWeight = clusterizer->GetEmcLogWeight();
203 }
204 else {
205 LogWeight = clusterizer->GetCpvLogWeight();
206 }
207
208 EvalLocalPosition(LogWeight,digits); // evaluate initial position
209}
210
211
212void AliPHOSEvalRecPoint::MakeJob()
213{
214 // Reconstruction algoritm implementation.
215
216 AliPHOSRecManager* recMng = GetReconstructionManager();
217
218 Init();
219
220 UnfoldLocalMaxima();
221
222 TObjArray children;
223 Int_t nChild = HasChild(children);
224
225 if(!nChild) {
226 EvaluatePosition();
227 if(Chi2Dof()>recMng->OneGamChisqCut()) SplitMergedShowers();
228 }
229
230 for(Int_t iChild=0; iChild<nChild; iChild++) {
231 AliPHOSEvalRecPoint* child = (AliPHOSEvalRecPoint*)children.At(iChild);
232 child->EvaluatePosition();
233
234 if(child->Chi2Dof()>recMng->OneGamChisqCut())
235 child->SplitMergedShowers();
236 }
237
238}
239
240void AliPHOSEvalRecPoint::InitTwoGam(Float_t* gamma1, Float_t* gamma2)
241{
242 //Compute start values for two gamma fit algorithm.
243 // gamma1[0], gamma1[1], gamma1[2] are Energy,x,z of the reconstructed "gammas".
244
245
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();
251
252 cout<<" (x,z,E)[old] = ("<<yy<<","<<xx<<","<<E<<")"<<endl;
253
254// xx = XY(xx/E);
255// yy = XY(yy/E);
256
257
258 Float_t eDigit ;
259 AliPHOSDigit * digit ;
260 Int_t nDigits = GetMultiplicity();
261 Int_t * Digits = GetDigitsList();
262 Float_t * Energies = GetEnergiesList();
263
264 Float_t ix ;
265 Float_t iy ;
266 Int_t relid[4] ;
267
268 Float_t exx = 0;
269 Float_t eyy = 0;
270 Float_t exy = 0;
271 Float_t sqr;
272 Float_t cos2fi = 1.;
273
274 AliPHOSGetter* fGetter = AliPHOSGetter::GetInstance();
275 const AliPHOSGeometry* fGeom = fGetter->PHOSGeometry();
276
c9eeb00e 277 Int_t iDigit; //loop variable
278
279 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++)
cbd576a6 280 {
281 digit = (AliPHOSDigit*)fGetter->Digits()->At( Digits[iDigit] );
282 eDigit = Energies[iDigit];
283 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
284 fGeom->RelPosInModule(relid, iy, ix);
285
286 Float_t dx = ix - xx;
287 Float_t dy = iy - yy;
288 exx += eDigit*dx*dx;
289 eyy += eDigit*dy*dy;
290 exy += eDigit*dx*dy;
291
292 }
293
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;
300
301 Float_t eu3 = 0;
c9eeb00e 302 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++)
cbd576a6 303 {
304 digit = (AliPHOSDigit*)fGetter->Digits()->At( Digits[iDigit] );
305 eDigit = Energies[iDigit];
306 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
307 fGeom->RelPosInModule(relid, iy, ix);
308
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;
313 }
314
315 Float_t c = E*eu3*eu3/(euu*euu*euu)/2.;
316 Float_t sign = 1.;
317 if(eu3<0) sign = -1.;
318 Float_t r = 1.+ c + sign*TMath::Sqrt(2.*c + c*c);
319 Float_t u = 0;
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);
322
323 Float_t e2c = E/(1.+r);
324 Float_t e1c = E-e2c;
325 Float_t u1 = -u/(1.+r);
326 Float_t u2 = u+u1;
327
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;
332
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);
339
340 gamma1[0] = e1c;
341 gamma1[1] = y1c;
342 gamma1[2] = x1c;
343
344 gamma2[0] = e2c;
345 gamma2[1] = y2c;
346 gamma2[2] = x2c;
347
348}
349
350void AliPHOSEvalRecPoint::TwoGam(Float_t* gamma1, Float_t* gamma2)
351{
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.
356
357
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.
363
364 Float_t chisq = 100.; //Initial chisquare.
365
366 Int_t nadc = GetMultiplicity();
367 if(nadc<3)
368 fChi2Dof= -111.;
369
370 Int_t dof = nadc - 5;
371 if(dof<1) dof=1;
372 Float_t chstop = chmin*dof;
373
374 Float_t ch = 1.e+20;
375 Float_t st = st0;
376 Float_t grx1 = 0.;
377 Float_t gry1 = 0.;
378 Float_t grx2 = 0.;
379 Float_t gry2 = 0.;
380 Float_t gre = 0.;
381 Float_t gr = 1.;
382 Float_t ee1=0,xx1=0,yy1=0,ee2=0,xx2=0,yy2=0,cosi;
383
384 Float_t e1c = gamma1[0];
385 Float_t y1c = gamma1[1];
386 Float_t x1c = gamma1[2];
387
388 Float_t e2c = gamma2[0];
389 Float_t y2c = gamma2[1];
390 Float_t x2c = gamma2[2];
391
392 Float_t E = GetEnergy();
393
394 Float_t eDigit ;
395 AliPHOSDigit * digit ;
396 Int_t nDigits = GetMultiplicity();
397 Int_t * Digits = GetDigitsList();
398 Float_t * Energies = GetEnergiesList();
399
400 Float_t ix ;
401 Float_t iy ;
402 Int_t relid[4] ;
403
404 AliPHOSGetter* fGetter = AliPHOSGetter::GetInstance();
405 const AliPHOSGeometry* fGeom = fGetter->PHOSGeometry();
406
407 for(Int_t Iter=0; Iter<Niter; Iter++)
408 {
409 Float_t chisqc = 0.;
410 Float_t grx1c = 0.;
411 Float_t gry1c = 0.;
412 Float_t grx2c = 0.;
413 Float_t gry2c = 0.;
414 Float_t grec = 0.;
415
416 for(Int_t iDigit = 0 ; iDigit < nDigits ; iDigit ++)
417 {
418 digit = (AliPHOSDigit*)fGetter->Digits()->At( Digits[iDigit] );
419 eDigit = Energies[iDigit];
420 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
421 fGeom->RelPosInModule(relid, iy, ix);
422
423 Float_t a1,gx1,gy1;
424 Float_t a2,gx2,gy2;
425
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);
431
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);
437
438 Float_t A = a1+a2;
439// Float_t D = Const*A*(1. - A/E);
440// if(D<0) D=0;
441
442// // D = 9.; // ????
443
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));
448
449 Float_t dd;
450 chisqc += GetReconstructionManager()->TwoGamChi2(A,eDigit,E,dd);
451
452 grx1c += gx1*dd;
453 gry1c += gy1*dd;
454 grx2c += gx2*dd;
455 gry2c += gy2*dd;
456 grec += (a1/e1c - a2/e2c)*E*dd;
457
458 }
459
460 Float_t grc = TMath::Sqrt(grx1c*grx1c + gry1c*gry1c + grx2c*grx2c + gry2c*gry2c);
461 if(grc<1.e-10) grc=1.e-10;
462
463 Float_t sc = 1. + chisqc/ch;
464 st = st/sc;
465
466 if(chisqc>ch)
467 goto loop20;
468 cosi = (grx1*grx1c + gry1*gry1c + grx2*grx2c + gry2*gry2c + gre*grec)/gr/grc;
469 st = st*sc/(1.4 - cosi);
470
471 ee1 = e1c;
472 xx1 = x1c;
473 yy1 = y1c;
474 ee2 = e2c;
475 xx2 = x2c;
476 yy2 = y2c;
477
478 ch = chisqc;
479
480 if(ch<chstop)
481 goto loop101;
482
483 grx1 = grx1c;
484 gry1 = gry1c;
485 grx2 = grx2c;
486 gry2 = gry2c;
487 gre = grec;
488 gr = grc;
489
490 loop20: ;
491 Float_t step = st*gr;
492
493
494 cout<<" Iteration "<<Iter<<" dof "<<dof<<" chisq/dof "<<ch/dof<<" chstop/dof "<<chstop/dof<<" step " <<step<<" stpmin "<<stpmin<<endl;
495
496
497 if(step<stpmin)
498 goto loop101;
499
500 Float_t de = st*gre*E;
501 e1c = ee1 - de;
502 e2c = ee2 + de;
503
504 if( (e1c>emin) && (e2c>emin) )
505 goto loop25;
506 st = st/2.;
507 goto loop20;
508
509 loop25: ;
510 x1c = xx1 - st*grx1;
511 y1c = yy1 - st*gry1;
512 x2c = xx2 - st*grx2;
513 y2c = yy2 - st*gry2;
514
515
516 }
517
518 loop101:
519
520// if(ch>chisq*(nadc-2)-delch)
521// return ch/dof;
522
523 chisq = ch/dof;
524 gamma1[0] = ee1;
525 gamma1[1] = yy1;
526 gamma1[2] = xx1;
527
528 gamma2[0] = ee2;
529 gamma2[1] = yy2;
530 gamma2[2] = xx2;
531
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;
538
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;
541
542 fChi2Dof = chisq;
543
544}
545
546void AliPHOSEvalRecPoint::UnfoldTwoMergedPoints(Float_t* gamma1, Float_t* gamma2)
547{
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.
551
552
553 Int_t Nmax = 2;
554 Float_t* gamma;
555
556 Int_t* Digits = GetDigitsList();
557 Int_t nDigits = GetMultiplicity();
558 Float_t* Energies = GetEnergiesList();
559 Float_t* eFit = new Float_t[nDigits];
560
561 AliPHOSGetter* fGetter = AliPHOSGetter::GetInstance();
562 const AliPHOSGeometry* fGeom = fGetter->PHOSGeometry();
563
564 for(Int_t iDigit=0; iDigit<nDigits; iDigit++)
565 {
566 AliPHOSDigit* digit = (AliPHOSDigit*)fGetter->Digits()->At( Digits[iDigit] );
567 Int_t relid[4] ;
568 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
569 Float_t x,z;
570 fGeom->RelPosInModule(relid, x, z);
571
572 Float_t gain = 0.;
573 for(Int_t iMax=0; iMax<Nmax; iMax++)
574 {
575 if(iMax==0)
576 gamma = gamma1;
577 else
578 gamma = gamma2;
579
580 Float_t eMax = gamma[0];
581 Float_t xMax = gamma[1];
582 Float_t zMax = gamma[2];
583
584 Float_t dx = xMax - x;
585 Float_t dz = zMax - z;
586 Float_t Amp,gx,gy;
587 GetReconstructionManager()->AG(eMax,dz,dx,Amp,gx,gy);
588 gain += Amp;
589 }
590
591 eFit[iDigit] = gain;
592 }
593
594
595 for( Int_t iMax=0; iMax<Nmax; iMax++)
596 {
597 if(iMax==0)
598 gamma = gamma1;
599 else
600 gamma = gamma2;
601
602 Float_t eMax = gamma[0];
603 Float_t xMax = gamma[1];
604 Float_t zMax = gamma[2];
605
606 AliPHOSEvalRecPoint* newRP = new AliPHOSEvalRecPoint(IsCPV(),this);
607 TVector3 newpos(xMax,0,zMax);
608 newRP->SetLocalPosition(newpos);
609
610 for( Int_t iDigit = 0 ; iDigit < nDigits ; iDigit ++)
611 {
612 AliPHOSDigit* digit = (AliPHOSDigit*)fGetter->Digits()->At( Digits[iDigit] );
613 Float_t eDigit = Energies[iDigit];
614 Int_t relid[4] ;
615 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
616 Float_t ix,iz;
617 fGeom->RelPosInModule(relid, ix, iz);
618
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);
628 }
629
630 cout<<"======= Split: daughter rec point "<<iMax<<" ================="<<endl;
631 newRP->Print("");
632
633 }
634
635 delete[] eFit;
636 cout<<endl;
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;
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);
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;
724 cout<<" chisq: "<<chisq<<endl;
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
743
744 cout<<" Iteration "<<Iter<<" dof "<<dof<<" chisq/dof "<<chisq/dof<<" chstop/dof "<<chstop/dof<<" step " <<step<<" Stpmin "<<Stpmin<<endl;
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);
777// cout<<" (x,z)[old] = ("<<x_old<<","<<z_old<<")"<<endl;
778// cout<<" (x,z)[new] = ("<<pos.X()<<","<<pos.Z()<<")"<<endl;
779
780
781// cout<<" gamma[1]: "<<gamma1[1]<<" gamma1[2]: "<<gamma1[2]<<endl;
782
783 fChi2Dof = chisq;
784
785}
786
787
788Bool_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 cout<<"+++++++ Killing this rec point ++++++++++"<<endl;
797 RemoveFromWorkingPool(this);
798 return kTRUE;
799 }
800
801 return kFALSE;
802}
803
804
805void 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
822void AliPHOSEvalRecPoint::MergeClosePoint()
823{
824
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();
829
830 for (Int_t i=0;i<InWorkingPool(); i++)
831 {
832 TObject* obj = GetFromWorkingPool(i);
366ed113 833 if(!((TObject*)this)->IsEqual(obj))
cbd576a6 834 {
835 AliPHOSRecPoint* rp = (AliPHOSRecPoint*)obj;
836 if(GetPHOSMod() == rp->GetPHOSMod())
837 {
838 if(TooClose(rp))
839 {
840 cout<<"+++++++ Merging point 1: ++++++"<<endl;
841 this->Print("");
842 cout<<"+++++++ and point 2: ++++++++++"<<endl;
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*)fGetter->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 cout<<"++++++ Resulting point: ++++++++"<<endl;
872 this->Print("");
873
874 break;
875 }
876 }
877 }
878 }
879
880}
881
882Int_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
a0636361 890 AliPHOSDigit * maxAt[1000];
cbd576a6 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 = fGetter->Clusterizer("PHOSclu-v1");
898 AliPHOSClusterizer* clusterizer = GetClusterizer();
899 if(!clusterizer) {
900 cout<<" Cannot get clusterizer. Exit."<<endl;
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 AliPHOSGetter* fGetter = AliPHOSGetter::GetInstance();
914 const AliPHOSGeometry* fGeom = fGetter->PHOSGeometry();
915 TClonesArray* digits = fGetter->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
c9eeb00e 926 Int_t iDigit; //loop variable
927
928 for(iDigit=0; iDigit<nDigits; iDigit++)
cbd576a6 929 {
930 eFit[iDigit] = 0.;
931 }
932
c9eeb00e 933 for(iDigit=0; iDigit<nDigits; iDigit++)
cbd576a6 934 {
935
936 AliPHOSDigit* digit = (AliPHOSDigit*)fGetter->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 {
a0636361 943 AliPHOSDigit* digitMax = maxAt[iMax];
cbd576a6 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 {
a0636361 959 AliPHOSDigit* digitMax = maxAt[iMax];
cbd576a6 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*)fGetter->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 cout<<" ratio -> "<<ratio<<endl;
985 eDigit = eDigit*ratio;
986 newRP->AddDigit(*digit,eDigit);
987 }
988 }
989
990 newRP->EvalLocalPosition(LogWeight,digits);
991 cout<<"======= Unfold: daughter rec point "<<iMax<<" ================="<<endl;
992 newRP->Print("");
993 }
994
995// RemoveFromWorkingPool(this);
996
997 delete[] eFit;
998 cout<<endl;
999
1000 return Nmax;
1001
1002}
1003
1004void AliPHOSEvalRecPoint::PrintPoint(Option_t* opt)
1005{
1006 AliPHOSCpvRecPoint::Print(opt);
1007
1008 TVector3 lpos;
1009 GetLocalPosition(lpos);
1010
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;
1015}
1016
1017AliPHOSRecManager* 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 cout<<" Couldn't find Reconstruction Manager. Exit."<<endl;
1025 exit(1);
1026 }
1027
1028 return recMng;
1029}
1030
1031
1032AliPHOSRecPoint* AliPHOSEvalRecPoint::Parent()
1033{
1034 if(fParent<0) return NULL;
1035 else
1036 return (AliPHOSRecPoint*)GetFromWorkingPool(fParent);
1037// return fParent;
1038}
1039
1040Float_t AliPHOSEvalRecPoint::Chi2Dof() const
1041{
1042 return fChi2Dof;
1043}
1044
1045const 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 cout<<" Couldn't find Working Pool. Exit."<<endl;
1053 exit(1);
1054 }
1055
1056 return wPool;
1057}
1058
1059void AliPHOSEvalRecPoint::AddToWorkingPool(TObject* obj)
1060{
1061 ((TObjArray*)GetWorkingPool())->Add(obj);
1062}
1063
1064TObject* AliPHOSEvalRecPoint::GetFromWorkingPool(Int_t i)
1065{
1066// return fWorkingPool->At(i);
1067 return ((TObjArray*)GetWorkingPool())->At(i);
1068}
1069
1070Int_t AliPHOSEvalRecPoint::InWorkingPool()
1071{
1072 return ((TObjArray*)GetWorkingPool())->GetEntries();
1073}
1074
1075void AliPHOSEvalRecPoint::RemoveFromWorkingPool(TObject* obj)
1076{
1077 ((TObjArray*)GetWorkingPool())->Remove(obj);
1078 ((TObjArray*)GetWorkingPool())->Compress();
1079}
1080
1081void AliPHOSEvalRecPoint::PrintWorkingPool()
1082{
1083 ((TObjArray*)GetWorkingPool())->Print("");
1084}
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099