Correct path of images and src
[u/mrichter/AliRoot.git] / PHOS / AliPHOSTrackSegmentMakerv1.cxx
CommitLineData
d15a28e7 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
b2a60966 16/* $Id$ */
17
d15a28e7 18//_________________________________________________________________________
b2a60966 19// Implementation version 1 of algorithm class to construct PHOS track segments
20// Associates EMC and PPSD clusters
21// Unfolds the EMC cluster
22//
23//*-- Author: Dmitri Peressounko (RRC Ki & SUBATECH)
24//
d15a28e7 25
26// --- ROOT system ---
27
28#include "TObjArray.h"
29#include "TClonesArray.h"
9f616d61 30#include "TObjectTable.h"
d15a28e7 31
32// --- Standard library ---
33
9f616d61 34#include <iostream>
35#include <cassert>
d15a28e7 36
37// --- AliRoot header files ---
38
39#include "AliPHOSTrackSegmentMakerv1.h"
40#include "AliPHOSTrackSegment.h"
41#include "AliPHOSLink.h"
42#include "AliPHOSv0.h"
43#include "AliRun.h"
44
92862013 45extern void UnfoldingChiSquare(Int_t &nPar, Double_t *Grad, Double_t & fret, Double_t *x, Int_t iflag) ;
9f616d61 46
d15a28e7 47ClassImp( AliPHOSTrackSegmentMakerv1)
48
49
50//____________________________________________________________________________
b2a60966 51 AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1() : AliPHOSTrackSegmentMaker()
d15a28e7 52{
9f616d61 53 // ctor
b2a60966 54
d15a28e7 55 fR0 = 4. ;
56 AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
57 //clusters are sorted in "rows" and "columns" of width geom->GetCrystalSize(0),
58 fDelta = fR0 + geom->GetCrystalSize(0) ;
6727be7e 59 fMinuit = new TMinuit(100) ;
60 fUnfoldFlag = kTRUE ;
d15a28e7 61}
62
9f616d61 63//____________________________________________________________________________
64 AliPHOSTrackSegmentMakerv1::~AliPHOSTrackSegmentMakerv1()
65{
66 // dtor
b2a60966 67
68 delete fMinuit ;
9f616d61 69}
b2a60966 70
d15a28e7 71//____________________________________________________________________________
72Bool_t AliPHOSTrackSegmentMakerv1::FindFit(AliPHOSEmcRecPoint * emcRP, int * maxAt, Float_t * maxAtEnergy,
92862013 73 Int_t nPar, Float_t * fitparameters)
9f616d61 74{
b2a60966 75 // Calls TMinuit to fit the energy distribution of a cluster with several maxima
6ad0bfa0 76
d15a28e7 77 AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
d15a28e7 78
9f616d61 79 gMinuit->SetPrintLevel(-1) ; // No Printout
80 gMinuit->SetFCN(UnfoldingChiSquare) ; // To set the address of the minimization function
81 gMinuit->SetObjectFit(emcRP) ; // To tranfer pointer to UnfoldingChiSquare
82
83 // filling initial values for fit parameters
d15a28e7 84 AliPHOSDigit * digit ;
9f616d61 85
86 Int_t ierflg = 0;
87 Int_t index = 0 ;
92862013 88 Int_t nDigits = (Int_t) nPar / 3 ;
9f616d61 89
d15a28e7 90 Int_t iDigit ;
9f616d61 91
92
92862013 93 for(iDigit = 0; iDigit < nDigits; iDigit++){
d15a28e7 94 digit = (AliPHOSDigit *) maxAt[iDigit];
95
92862013 96 Int_t relid[4] ;
d15a28e7 97 Float_t x ;
98 Float_t z ;
92862013 99 geom->AbsToRelNumbering(digit->GetId(), relid) ;
100 geom->RelPosInModule(relid, x, z) ;
d15a28e7 101
92862013 102 Float_t energy = maxAtEnergy[iDigit] ;
d15a28e7 103
9f616d61 104 gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ;
105 index++ ;
d15a28e7 106 if(ierflg != 0){
107 cout << "PHOS Unfolding> Unable to set initial value for fit procedure : x = " << x << endl ;
108 return kFALSE;
109 }
9f616d61 110 gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ;
111 index++ ;
d15a28e7 112 if(ierflg != 0){
113 cout << "PHOS Unfolding> Unable to set initial value for fit procedure : z = " << z << endl ;
114 return kFALSE;
115 }
92862013 116 gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;
9f616d61 117 index++ ;
d15a28e7 118 if(ierflg != 0){
92862013 119 cout << "PHOS Unfolding> Unable to set initial value for fit procedure : energy = " << energy << endl ;
d15a28e7 120 return kFALSE;
6ad0bfa0 121 }
d15a28e7 122 }
123
9f616d61 124 Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
125 // depends on it.
126 Double_t p1 = 1.0 ;
127 Double_t p2 = 0.0 ;
128
129 gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TgMinuit to reduce function calls
130 gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
d15a28e7 131 gMinuit->SetMaxIterations(5);
9f616d61 132 gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
133 gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
134 if(ierflg == 4){ // Minimum not found
135 cout << "PHOS Unfolding> Fit not converged, cluster abandoned "<< endl ;
d15a28e7 136 return kFALSE ;
137 }
92862013 138 for(index = 0; index < nPar; index++){
d15a28e7 139 Double_t err ;
140 Double_t val ;
9f616d61 141 gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index
92862013 142 fitparameters[index] = val ;
d15a28e7 143 }
d15a28e7 144 return kTRUE;
9f616d61 145
d15a28e7 146}
9f616d61 147
d15a28e7 148//____________________________________________________________________________
b2a60966 149void AliPHOSTrackSegmentMakerv1::FillOneModule(DigitsList * dl, RecPointsList * emcIn, TObjArray * emcOut,
d15a28e7 150 RecPointsList * ppsdIn, TObjArray * ppsdOutUp,
92862013 151 TObjArray * ppsdOutLow, Int_t & phosmod, Int_t & emcStopedAt,
d15a28e7 152 Int_t & ppsdStopedAt)
9f616d61 153{
154 // Unfold clusters and fill xxxOut arrays with clusters from one PHOS module
155
156 AliPHOSEmcRecPoint * emcRecPoint ;
d15a28e7 157 AliPHOSPpsdRecPoint * ppsdRecPoint ;
158 Int_t index ;
6ad0bfa0 159
92862013 160 Int_t nEmcUnfolded = emcIn->GetEntries() ;
161 for(index = emcStopedAt; index < nEmcUnfolded; index++){
d15a28e7 162 emcRecPoint = (AliPHOSEmcRecPoint *) (*emcIn)[index] ;
6ad0bfa0 163
92862013 164 if(emcRecPoint->GetPHOSMod() != phosmod )
9f616d61 165 break ;
166
92862013 167 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
168 Int_t maxAt[nMultipl] ;
169 Float_t maxAtEnergy[nMultipl] ;
170 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy) ;
6a3f1304 171
6727be7e 172 if(nMax <= 1 ) // if cluster is very flat (no pronounced maximum) then nMax = 0
d15a28e7 173 emcOut->Add(emcRecPoint) ;
6727be7e 174 else if (fUnfoldFlag) {
b2a60966 175 UnfoldClusters(dl, emcIn, emcRecPoint, nMax, maxAt, maxAtEnergy, emcOut) ;
d15a28e7 176 emcIn->Remove(emcRecPoint);
177 emcIn->Compress() ;
92862013 178 nEmcUnfolded-- ;
d15a28e7 179 index-- ;
180 }
181 }
182 emcStopedAt = index ;
183
184 for(index = ppsdStopedAt; index < ppsdIn->GetEntries(); index++){
185 ppsdRecPoint = (AliPHOSPpsdRecPoint *) (*ppsdIn)[index] ;
92862013 186 if(ppsdRecPoint->GetPHOSMod() != phosmod )
9f616d61 187 break ;
d15a28e7 188 if(ppsdRecPoint->GetUp() )
189 ppsdOutUp->Add(ppsdRecPoint) ;
190 else
191 ppsdOutLow->Add(ppsdRecPoint) ;
192 }
193 ppsdStopedAt = index ;
6ad0bfa0 194
d15a28e7 195 emcOut->Sort() ;
196 ppsdOutUp->Sort() ;
6ad0bfa0 197 ppsdOutLow->Sort() ;
d15a28e7 198}
199//____________________________________________________________________________
92862013 200Float_t AliPHOSTrackSegmentMakerv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcclu,AliPHOSPpsdRecPoint * PpsdClu, Bool_t &toofar)
d15a28e7 201{
b2a60966 202 // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
203
92862013 204 Float_t r = fR0 ;
d15a28e7 205
206 TVector3 vecEmc ;
207 TVector3 vecPpsd ;
208
92862013 209 emcclu->GetLocalPosition(vecEmc) ;
d15a28e7 210 PpsdClu->GetLocalPosition(vecPpsd) ;
92862013 211 if(emcclu->GetPHOSMod() == PpsdClu->GetPHOSMod()){
d15a28e7 212 if(vecPpsd.X() >= vecEmc.X() - fDelta ){
213 if(vecPpsd.Z() >= vecEmc.Z() - fDelta ){
214 AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
6ad0bfa0 215 // Correct to difference in CPV and EMC position due to different distance to center.
216 // we assume, that particle moves from center
92862013 217 Float_t dCPV = geom->GetIPtoOuterCoverDistance();
218 Float_t dEMC = geom->GetIPtoCrystalSurface() ;
219 dEMC = dEMC / dCPV ;
220 vecPpsd = dEMC * vecPpsd - vecEmc ;
221 r = vecPpsd.Mag() ;
d15a28e7 222 } // if zPpsd >= zEmc - fDelta
92862013 223 toofar = kFALSE ;
d15a28e7 224 } // if xPpsd >= xEmc - fDelta
225 else
92862013 226 toofar = kTRUE ;
d15a28e7 227 }
228 else
92862013 229 toofar = kTRUE ;
d15a28e7 230
92862013 231 return r ;
d15a28e7 232}
233
234//____________________________________________________________________________
92862013 235void AliPHOSTrackSegmentMakerv1::MakeLinks(TObjArray * emcRecPoints, TObjArray * ppsdRecPointsUp,
236 TObjArray * ppsdRecPointsLow, TClonesArray * linklowArray,
237 TClonesArray *linkupArray)
d15a28e7 238{
b2a60966 239 // Finds distances (links) between all EMC and PPSD clusters, which are not further apart from each other than fR0
240
92862013 241 TIter nextEmc(emcRecPoints) ;
d15a28e7 242 Int_t iEmcClu = 0 ;
243
92862013 244 AliPHOSPpsdRecPoint * ppsdlow ;
245 AliPHOSPpsdRecPoint * ppsdup ;
246 AliPHOSEmcRecPoint * emcclu ;
d15a28e7 247
248 Int_t iLinkLow = 0 ;
249 Int_t iLinkUp = 0 ;
250
92862013 251 while( (emcclu = (AliPHOSEmcRecPoint*)nextEmc() ) ) {
252 Bool_t toofar ;
253 TIter nextPpsdLow(ppsdRecPointsLow ) ;
d15a28e7 254 Int_t iPpsdLow = 0 ;
255
92862013 256 while( (ppsdlow = (AliPHOSPpsdRecPoint*)nextPpsdLow() ) ) {
257 Float_t r = GetDistanceInPHOSPlane(emcclu, ppsdlow, toofar) ;
d15a28e7 258
92862013 259 if(toofar)
d15a28e7 260 break ;
92862013 261 if(r < fR0)
262 new( (*linklowArray)[iLinkLow++]) AliPHOSLink(r, iEmcClu, iPpsdLow) ;
d15a28e7 263
264 iPpsdLow++ ;
265
266 }
267
92862013 268 TIter nextPpsdUp(ppsdRecPointsUp ) ;
d15a28e7 269 Int_t iPpsdUp = 0 ;
270
92862013 271 while( (ppsdup = (AliPHOSPpsdRecPoint*)nextPpsdUp() ) ) {
272 Float_t r = GetDistanceInPHOSPlane(emcclu, ppsdup, toofar) ;
d15a28e7 273
92862013 274 if(toofar)
d15a28e7 275 break ;
92862013 276 if(r < fR0)
277 new( (*linkupArray)[iLinkUp++]) AliPHOSLink(r, iEmcClu, iPpsdUp) ;
d15a28e7 278
279 iPpsdUp++ ;
280
281 }
282
283 iEmcClu++ ;
284
285 } // while nextEmC
286
92862013 287 linklowArray->Sort() ; //first links with smallest distances
288 linkupArray->Sort() ;
d15a28e7 289}
290
291//____________________________________________________________________________
92862013 292void AliPHOSTrackSegmentMakerv1::MakePairs(TObjArray * emcRecPoints, TObjArray * ppsdRecPointsUp,
293 TObjArray * ppsdRecPointsLow, TClonesArray * linklowArray,
294 TClonesArray * linkupArray, TrackSegmentsList * trsl)
6ad0bfa0 295{
296
297 // Finds the smallest links and makes pairs of PPSD and EMC clusters with smallest distance
298
92862013 299 TIter nextLow(linklowArray) ;
300 TIter nextUp(linkupArray) ;
d15a28e7 301
302 AliPHOSLink * linkLow ;
303 AliPHOSLink * linkUp ;
304
305 AliPHOSEmcRecPoint * emc ;
306 AliPHOSPpsdRecPoint * ppsdLow ;
307 AliPHOSPpsdRecPoint * ppsdUp ;
308
92862013 309 AliPHOSRecPoint * nullpointer = 0 ;
9f616d61 310
d15a28e7 311 while ( (linkLow = (AliPHOSLink *)nextLow() ) ){
92862013 312 emc = (AliPHOSEmcRecPoint *) emcRecPoints->At(linkLow->GetEmc()) ;
313 ppsdLow = (AliPHOSPpsdRecPoint *) ppsdRecPointsLow->At(linkLow->GetPpsd()) ;
6ad0bfa0 314 if( (emc) && (ppsdLow) ){ // RecPoints not removed yet
d15a28e7 315 ppsdUp = 0 ;
316
317 while ( (linkUp = (AliPHOSLink *)nextUp() ) ){
318 if(linkLow->GetEmc() == linkUp->GetEmc() ){
92862013 319 ppsdUp = (AliPHOSPpsdRecPoint *) ppsdRecPointsUp->At(linkUp->GetPpsd()) ;
d15a28e7 320 break ;
321 }
322
323 } // while nextUp
324
325 nextUp.Reset();
b2a60966 326// AliPHOSTrackSegment * subtr = new AliPHOSTrackSegment(emc, ppsdUp, ppsdLow ) ;
327// trsl->Add(subtr) ;
328 new( (*trsl)[fNTrackSegments] ) AliPHOSTrackSegment(emc, ppsdUp, ppsdLow ) ;
329 fNTrackSegments++ ;
92862013 330 emcRecPoints->AddAt(nullpointer,linkLow->GetEmc()) ;
331 ppsdRecPointsLow->AddAt(nullpointer,linkLow->GetPpsd()) ;
d15a28e7 332
333 if(ppsdUp)
92862013 334 ppsdRecPointsUp->AddAt(nullpointer,linkUp->GetPpsd()) ;
d15a28e7 335
6ad0bfa0 336 }
d15a28e7 337 }
338
92862013 339 TIter nextEmc(emcRecPoints) ;
d15a28e7 340 nextEmc.Reset() ;
341
92862013 342 while( (emc = (AliPHOSEmcRecPoint*)nextEmc()) ){ //to create pairs if no ppsdlow
9f616d61 343 ppsdLow = 0 ;
344 ppsdUp = 0 ;
d15a28e7 345
346 while ( (linkUp = (AliPHOSLink *)nextUp() ) ){
347
92862013 348 if(emcRecPoints->IndexOf(emc) == linkUp->GetEmc() ){
349 ppsdUp = (AliPHOSPpsdRecPoint *) ppsdRecPointsUp->At(linkUp->GetPpsd()) ;
d15a28e7 350 break ;
351 }
352
353 }
b2a60966 354// AliPHOSTrackSegment * subtr = new AliPHOSTrackSegment(emc, ppsdUp, ppsdLow ) ;
355// trsl->Add(subtr) ;
356 new( (*trsl)[fNTrackSegments] ) AliPHOSTrackSegment(emc, ppsdUp, ppsdLow ) ;
357 fNTrackSegments++ ;
358
359
d15a28e7 360 if(ppsdUp)
92862013 361 ppsdRecPointsUp->AddAt(nullpointer,linkUp->GetPpsd()) ;
d15a28e7 362 }
363
364}
365
366//____________________________________________________________________________
b2a60966 367void AliPHOSTrackSegmentMakerv1::MakeTrackSegments(DigitsList * dl, RecPointsList * emcl,
d15a28e7 368 RecPointsList * ppsdl, TrackSegmentsList * trsl)
369{
b2a60966 370 // Makes the track segments out of the list of EMC and PPSD Recpoints and stores them in a list
d15a28e7 371
92862013 372 Int_t phosmod = 1 ;
9f616d61 373 Int_t emcStopedAt = 0 ;
d15a28e7 374 Int_t ppsdStopedAt = 0 ;
375
92862013 376 TObjArray * emcRecPoints = new TObjArray(100) ; // these arrays keep pointers
377 TObjArray * ppsdRecPointsUp = new TObjArray(100) ; // to RecPoints, which are
378 TObjArray * ppsdRecPointsLow = new TObjArray(100) ; // kept in TClonesArray's emcl and ppsdl
d15a28e7 379
380
92862013 381 TClonesArray * linklowArray = new TClonesArray("AliPHOSLink", 100);
382 TClonesArray * linkupArray = new TClonesArray("AliPHOSLink", 100);
d15a28e7 383
384 AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
385
92862013 386 while(phosmod <= geom->GetNModules() ){
d15a28e7 387
b2a60966 388 FillOneModule(dl, emcl, emcRecPoints, ppsdl, ppsdRecPointsUp, ppsdRecPointsLow, phosmod, emcStopedAt, ppsdStopedAt) ;
6ad0bfa0 389
92862013 390 MakeLinks(emcRecPoints, ppsdRecPointsUp, ppsdRecPointsLow, linklowArray, linkupArray) ;
d15a28e7 391
92862013 392 MakePairs(emcRecPoints, ppsdRecPointsUp, ppsdRecPointsLow, linklowArray, linkupArray, trsl) ;
d15a28e7 393
92862013 394 emcRecPoints->Clear() ;
9f616d61 395
92862013 396 ppsdRecPointsUp->Clear() ;
9f616d61 397
92862013 398 ppsdRecPointsLow->Clear() ;
9f616d61 399
92862013 400 linkupArray->Clear() ;
9f616d61 401
92862013 402 linklowArray->Clear() ;
9f616d61 403
92862013 404 phosmod++ ;
d15a28e7 405 }
92862013 406 delete emcRecPoints ;
407 emcRecPoints = 0 ;
9f616d61 408
92862013 409 delete ppsdRecPointsUp ;
410 ppsdRecPointsUp = 0 ;
9f616d61 411
92862013 412 delete ppsdRecPointsLow ;
413 ppsdRecPointsLow = 0 ;
9f616d61 414
92862013 415 delete linkupArray ;
416 linkupArray = 0 ;
9f616d61 417
92862013 418 delete linklowArray ;
419 linklowArray = 0 ;
d15a28e7 420}
421
422//____________________________________________________________________________
9f616d61 423Double_t AliPHOSTrackSegmentMakerv1::ShowerShape(Double_t r)
424{
b2a60966 425 // Shape of the shower (see PHOS TDR)
426 // If you change this function, change also the gradien evaluation in ChiSquare()
427
6ad0bfa0 428 Double_t r4 = r*r*r*r ;
429 Double_t r295 = TMath::Power(r, 2.95) ;
d15a28e7 430 Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
431 return shape ;
432}
433
434//____________________________________________________________________________
b2a60966 435void AliPHOSTrackSegmentMakerv1::UnfoldClusters(DigitsList * dl, RecPointsList * emcIn, AliPHOSEmcRecPoint * iniEmc,
92862013 436 Int_t nMax, int * maxAt, Float_t * maxAtEnergy, TObjArray * emcList)
d15a28e7 437{
b2a60966 438 // Performs the unfolding of a cluster with nMax overlapping showers
439 // This is time consuming (use the (Un)SetUnfolFlag() )
440
92862013 441 Int_t nPar = 3 * nMax ;
442 Float_t fitparameters[nPar] ;
d15a28e7 443 AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
9f616d61 444
92862013 445 Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
9f616d61 446 if( !rv ) // Fit failed, return and remove cluster
d15a28e7 447 return ;
448
449 Float_t xDigit ;
450 Float_t zDigit ;
92862013 451 Int_t relid[4] ;
d15a28e7 452
92862013 453 Int_t nDigits = iniEmc->GetMultiplicity() ;
d15a28e7 454 Float_t xpar ;
455 Float_t zpar ;
92862013 456 Float_t epar ;
457 Float_t distance ;
458 Float_t ratio ;
459 Float_t efit[nDigits] ;
d15a28e7 460 Int_t iparam ;
461 Int_t iDigit ;
462
463 AliPHOSDigit * digit ;
464 AliPHOSEmcRecPoint * emcRP ;
6ad0bfa0 465 Int_t * emcDigits = iniEmc->GetDigitsList() ;
d15a28e7 466 Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
467
468 Int_t iRecPoint = emcIn->GetEntries() ;
469
92862013 470 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
d15a28e7 471 digit = (AliPHOSDigit *) emcDigits[iDigit];
92862013 472 geom->AbsToRelNumbering(digit->GetId(), relid) ;
473 geom->RelPosInModule(relid, xDigit, zDigit) ;
474 efit[iDigit] = 0;
d15a28e7 475 iparam = 0 ;
476
92862013 477 while(iparam < nPar ){
478 xpar = fitparameters[iparam] ;
479 zpar = fitparameters[iparam+1] ;
480 epar = fitparameters[iparam+2] ;
d15a28e7 481 iparam += 3 ;
92862013 482 distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ;
483 distance = TMath::Sqrt(distance) ;
b2a60966 484 efit[iDigit] += epar * AliPHOSTrackSegmentMakerv1::ShowerShape(distance) ;
d15a28e7 485 }
d15a28e7 486 }
487
488 iparam = 0 ;
489 Float_t eDigit ;
490
92862013 491 while(iparam < nPar ){
492 xpar = fitparameters[iparam] ;
493 zpar = fitparameters[iparam+1] ;
494 epar = fitparameters[iparam+2] ;
d15a28e7 495 iparam += 3 ;
496 new ((*emcIn)[iRecPoint]) AliPHOSEmcRecPoint( iniEmc->GetLogWeightCut(), iniEmc->GetLocMaxCut() ) ;
497 emcRP = (AliPHOSEmcRecPoint *) emcIn->At(iRecPoint++);
498
92862013 499 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
d15a28e7 500 digit = (AliPHOSDigit *) emcDigits[iDigit];
92862013 501 geom->AbsToRelNumbering(digit->GetId(), relid) ;
502 geom->RelPosInModule(relid, xDigit, zDigit) ;
503 distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ;
504 distance = TMath::Sqrt(distance) ;
b2a60966 505 ratio = epar * AliPHOSTrackSegmentMakerv1::ShowerShape(distance) / efit[iDigit] ;
92862013 506 eDigit = emcEnergies[iDigit] * ratio ;
9f616d61 507 emcRP->AddDigit( *digit, eDigit ) ;
d15a28e7 508 }
509
510 emcList->Add(emcRP) ;
511
512 }
513}
6ad0bfa0 514
d15a28e7 515//______________________________________________________________________________
92862013 516void UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
d15a28e7 517{
b2a60966 518 // Calculates th Chi square for the cluster unfolding minimization
519 // Number of parameters, Gradient, Chi squared, parameters, what to do
9f616d61 520
d15a28e7 521 AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
d15a28e7 522
9f616d61 523 AliPHOSEmcRecPoint * emcRP = (AliPHOSEmcRecPoint *) gMinuit->GetObjectFit() ; // EmcRecPoint to fit
6ad0bfa0 524 Int_t * emcDigits = emcRP->GetDigitsList() ;
d15a28e7 525 Float_t * emcEnergies = emcRP->GetEnergiesList() ;
526 fret = 0. ;
527 Int_t iparam ;
528
9f616d61 529 if(iflag == 2)
92862013 530 for(iparam = 0 ; iparam < nPar ; iparam++)
9f616d61 531 Grad[iparam] = 0 ; // Will evaluate gradient
532
92862013 533 Double_t efit ;
9f616d61 534
d15a28e7 535 AliPHOSDigit * digit ;
536 Int_t iDigit = 0 ;
9f616d61 537
d15a28e7 538 while ( (digit = (AliPHOSDigit *)emcDigits[iDigit] )){
92862013 539 Int_t relid[4] ;
d15a28e7 540 Float_t xDigit ;
541 Float_t zDigit ;
92862013 542 geom->AbsToRelNumbering(digit->GetId(), relid) ;
543 geom->RelPosInModule(relid, xDigit, zDigit) ;
d15a28e7 544
9f616d61 545 if(iflag == 2){ // calculate gradient
546 Int_t iParam = 0 ;
92862013 547 efit = 0 ;
548 while(iParam < nPar ){
549 Double_t distance = (xDigit - x[iParam]) * (xDigit - x[iParam]) ;
9f616d61 550 iParam++ ;
92862013 551 distance += (zDigit - x[iParam]) * (zDigit - x[iParam]) ;
552 distance = TMath::Sqrt( distance ) ;
9f616d61 553 iParam++ ;
92862013 554 efit += x[iParam] * AliPHOSTrackSegmentMakerv1::ShowerShape(distance) ;
9f616d61 555 iParam++ ;
556 }
92862013 557 Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
9f616d61 558 iParam = 0 ;
92862013 559 while(iParam < nPar ){
9f616d61 560 Double_t xpar = x[iParam] ;
561 Double_t zpar = x[iParam+1] ;
92862013 562 Double_t epar = x[iParam+2] ;
9f616d61 563 Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
564 Double_t shape = sum * AliPHOSTrackSegmentMakerv1::ShowerShape(dr) ;
565 Double_t r4 = dr*dr*dr*dr ;
566 Double_t r295 = TMath::Power(dr,2.95) ;
567 Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
568 0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
569
92862013 570 Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x
9f616d61 571 iParam++ ;
92862013 572 Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z
9f616d61 573 iParam++ ;
574 Grad[iParam] += shape ; // Derivative over energy
575 iParam++ ;
576 }
577 }
92862013 578 efit = 0;
9f616d61 579 iparam = 0 ;
92862013 580 while(iparam < nPar ){
9f616d61 581 Double_t xpar = x[iparam] ;
582 Double_t zpar = x[iparam+1] ;
92862013 583 Double_t epar = x[iparam+2] ;
9f616d61 584 iparam += 3 ;
92862013 585 Double_t distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ;
586 distance = TMath::Sqrt(distance) ;
587 efit += epar * AliPHOSTrackSegmentMakerv1::ShowerShape(distance) ;
9f616d61 588 }
92862013 589 fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
9f616d61 590 // Here we assume, that sigma = sqrt(E)
591 iDigit++ ;
d15a28e7 592 }
593}