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