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