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