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