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