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