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 | |
9a1398dd |
18 | //*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute) |
d15a28e7 |
19 | ////////////////////////////////////////////////////////////////////////////// |
9a1398dd |
20 | // Clusterization class. Performs clusterization (collects neighbouring active cells) and |
a4e98857 |
21 | // unfolds the clusters having several local maxima. |
22 | // Results are stored in TreeR#, branches PHOSEmcRP (EMC recPoints), |
9a1398dd |
23 | // PHOSCpvRP (CPV RecPoints) and AliPHOSClusterizer (Clusterizer with all |
f035f6ce |
24 | // parameters including input digits branch title, thresholds etc.) |
a4e98857 |
25 | // This TTask is normally called from Reconstructioner, but can as well be used in |
26 | // standalone mode. |
27 | // Use Case: |
fc12304f |
28 | // root [0] AliPHOSClusterizerv1 * cl = new AliPHOSClusterizerv1("galice.root", "recpointsname", "digitsname") |
a4e98857 |
29 | // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated |
fc12304f |
30 | // // reads gAlice from header file "galice.root", uses digits stored in the branch names "digitsname" (default = "Default") |
31 | // // and saves recpoints in branch named "recpointsname" (default = "digitsname") |
a4e98857 |
32 | // root [1] cl->ExecuteTask() |
9a1398dd |
33 | // //finds RecPoints in all events stored in galice.root |
a4e98857 |
34 | // root [2] cl->SetDigitsBranch("digits2") |
f035f6ce |
35 | // //sets another title for Digitis (input) branch |
a4e98857 |
36 | // root [3] cl->SetRecPointsBranch("recp2") |
f035f6ce |
37 | // //sets another title four output branches |
a4e98857 |
38 | // root [4] cl->SetEmcLocalMaxCut(0.03) |
9a1398dd |
39 | // //set clusterization parameters |
a4e98857 |
40 | // root [5] cl->ExecuteTask("deb all time") |
9a1398dd |
41 | // //once more finds RecPoints options are |
42 | // // deb - print number of found rec points |
43 | // // deb all - print number of found RecPoints and some their characteristics |
44 | // // time - print benchmarking results |
ed4205d8 |
45 | |
d15a28e7 |
46 | // --- ROOT system --- |
47 | |
48 | #include "TMath.h" |
9a1398dd |
49 | #include "TMinuit.h" |
50 | #include "TTree.h" |
9a1398dd |
51 | #include "TBenchmark.h" |
d15a28e7 |
52 | |
53 | // --- Standard library --- |
54 | |
d15a28e7 |
55 | // --- AliRoot header files --- |
88cb7938 |
56 | #include "AliPHOSGetter.h" |
e957fea8 |
57 | #include "AliPHOSGeometry.h" |
d15a28e7 |
58 | #include "AliPHOSClusterizerv1.h" |
88cb7938 |
59 | #include "AliPHOSEmcRecPoint.h" |
9a1398dd |
60 | #include "AliPHOSCpvRecPoint.h" |
d15a28e7 |
61 | #include "AliPHOSDigit.h" |
9a1398dd |
62 | #include "AliPHOSDigitizer.h" |
d15a28e7 |
63 | |
64 | ClassImp(AliPHOSClusterizerv1) |
f035f6ce |
65 | |
d15a28e7 |
66 | //____________________________________________________________________________ |
7b7c1533 |
67 | AliPHOSClusterizerv1::AliPHOSClusterizerv1() : AliPHOSClusterizer() |
d15a28e7 |
68 | { |
f035f6ce |
69 | // default ctor (to be used mainly by Streamer) |
f035f6ce |
70 | |
8d0f3f77 |
71 | InitParameters() ; |
92f521a9 |
72 | fDefaultInit = kTRUE ; |
9a1398dd |
73 | } |
7b7c1533 |
74 | |
9a1398dd |
75 | //____________________________________________________________________________ |
88cb7938 |
76 | AliPHOSClusterizerv1::AliPHOSClusterizerv1(const TString alirunFileName, const TString eventFolderName) |
77 | :AliPHOSClusterizer(alirunFileName, eventFolderName) |
9a1398dd |
78 | { |
a4e98857 |
79 | // ctor with the indication of the file where header Tree and digits Tree are stored |
9a1398dd |
80 | |
8d0f3f77 |
81 | InitParameters() ; |
2bd5457f |
82 | Init() ; |
92f521a9 |
83 | fDefaultInit = kFALSE ; |
9a1398dd |
84 | } |
fc12304f |
85 | |
9688c1dd |
86 | //____________________________________________________________________________ |
87 | AliPHOSClusterizerv1::~AliPHOSClusterizerv1() |
88 | { |
0bc3b8ed |
89 | // dtor |
90 | |
88cb7938 |
91 | // delete fPedestals ; |
92 | // delete fGains ; |
9688c1dd |
93 | } |
fc12304f |
94 | |
88cb7938 |
95 | |
fc12304f |
96 | //____________________________________________________________________________ |
97 | const TString AliPHOSClusterizerv1::BranchName() const |
98 | { |
88cb7938 |
99 | return GetName(); |
fc12304f |
100 | } |
101 | |
9a1398dd |
102 | //____________________________________________________________________________ |
3758d9fc |
103 | Float_t AliPHOSClusterizerv1::Calibrate(Int_t amp, Int_t absId) const |
88cb7938 |
104 | { |
105 | //To be replaced later by the method, reading individual parameters from the database |
106 | // if(fCalibrationDB) |
107 | // return fCalibrationDB->Calibrate(amp,absId) ; |
108 | // else{ //simulation |
109 | if(absId <= fEmcCrystals) //calibrate as EMC |
110 | return fADCpedestalEmc + amp*fADCchanelEmc ; |
111 | else //calibrate as CPV |
112 | return fADCpedestalCpv+ amp*fADCchanelCpv ; |
113 | // } |
3758d9fc |
114 | } |
548f0134 |
115 | |
3758d9fc |
116 | //____________________________________________________________________________ |
a4e98857 |
117 | void AliPHOSClusterizerv1::Exec(Option_t * option) |
118 | { |
7b7c1533 |
119 | // Steering method |
9a1398dd |
120 | |
121 | if(strstr(option,"tim")) |
7d493c2c |
122 | gBenchmark->Start("PHOSClusterizer"); |
9a1398dd |
123 | |
7b7c1533 |
124 | if(strstr(option,"print")) |
88cb7938 |
125 | Print() ; |
bed9e3fb |
126 | |
88cb7938 |
127 | AliPHOSGetter * gime = AliPHOSGetter::Instance() ; |
128 | |
fbf811ec |
129 | Int_t nevents = gime->MaxEvent() ; |
7b7c1533 |
130 | Int_t ievent ; |
fbf811ec |
131 | |
88cb7938 |
132 | for(ievent = 0; ievent < nevents; ievent++) |
133 | { |
134 | gime->Event(ievent, "D"); |
135 | |
136 | GetCalibrationParameters() ; |
bed9e3fb |
137 | |
548f0134 |
138 | fNumberOfEmcClusters = fNumberOfCpvClusters = 0 ; |
9a1398dd |
139 | |
88cb7938 |
140 | MakeClusters() ; |
141 | |
fc12304f |
142 | if(fToUnfold) |
143 | MakeUnfolding() ; |
7b7c1533 |
144 | |
88cb7938 |
145 | WriteRecPoints(); |
7b7c1533 |
146 | |
94de8339 |
147 | if(strstr(option,"deb")) |
148 | PrintRecPoints(option) ; |
149 | |
88cb7938 |
150 | //increment the total number of recpoints per run |
fbf811ec |
151 | fRecPointsInRun += gime->EmcRecPoints()->GetEntriesFast() ; |
152 | fRecPointsInRun += gime->CpvRecPoints()->GetEntriesFast() ; |
88cb7938 |
153 | } |
154 | |
155 | Unload(); |
9a1398dd |
156 | |
157 | if(strstr(option,"tim")){ |
158 | gBenchmark->Stop("PHOSClusterizer"); |
21cd0c07 |
159 | Info("Exec", " took %f seconds for Clusterizing %f seconds per event \n", |
160 | gBenchmark->GetCpuTime("PHOSClusterizer"), |
161 | gBenchmark->GetCpuTime("PHOSClusterizer")/nevents ) ; |
88cb7938 |
162 | } |
9a1398dd |
163 | } |
164 | |
165 | //____________________________________________________________________________ |
a0636361 |
166 | Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** maxAt, Float_t * maxAtEnergy, |
d1de15f5 |
167 | Int_t nPar, Float_t * fitparameters) const |
9a1398dd |
168 | { |
169 | // Calls TMinuit to fit the energy distribution of a cluster with several maxima |
a4e98857 |
170 | // The initial values for fitting procedure are set equal to the positions of local maxima. |
f035f6ce |
171 | // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers |
9a1398dd |
172 | |
88cb7938 |
173 | |
174 | AliPHOSGetter * gime = AliPHOSGetter::Instance(); |
175 | TClonesArray * digits = gime->Digits(); |
7b7c1533 |
176 | |
177 | |
9a1398dd |
178 | gMinuit->mncler(); // Reset Minuit's list of paramters |
179 | gMinuit->SetPrintLevel(-1) ; // No Printout |
180 | gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ; |
181 | // To set the address of the minimization function |
182 | |
183 | TList * toMinuit = new TList(); |
184 | toMinuit->AddAt(emcRP,0) ; |
7b7c1533 |
185 | toMinuit->AddAt(digits,1) ; |
9a1398dd |
186 | |
187 | gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare |
188 | |
189 | // filling initial values for fit parameters |
190 | AliPHOSDigit * digit ; |
191 | |
192 | Int_t ierflg = 0; |
193 | Int_t index = 0 ; |
194 | Int_t nDigits = (Int_t) nPar / 3 ; |
195 | |
196 | Int_t iDigit ; |
197 | |
7b7c1533 |
198 | const AliPHOSGeometry * geom = gime->PHOSGeometry() ; |
9a1398dd |
199 | |
200 | for(iDigit = 0; iDigit < nDigits; iDigit++){ |
a0636361 |
201 | digit = maxAt[iDigit]; |
9a1398dd |
202 | |
203 | Int_t relid[4] ; |
7b7c1533 |
204 | Float_t x = 0.; |
205 | Float_t z = 0.; |
206 | geom->AbsToRelNumbering(digit->GetId(), relid) ; |
207 | geom->RelPosInModule(relid, x, z) ; |
9a1398dd |
208 | |
209 | Float_t energy = maxAtEnergy[iDigit] ; |
210 | |
211 | gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ; |
212 | index++ ; |
213 | if(ierflg != 0){ |
21cd0c07 |
214 | Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ; |
9a1398dd |
215 | return kFALSE; |
216 | } |
217 | gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ; |
218 | index++ ; |
219 | if(ierflg != 0){ |
21cd0c07 |
220 | Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ; |
9a1398dd |
221 | return kFALSE; |
222 | } |
223 | gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ; |
224 | index++ ; |
225 | if(ierflg != 0){ |
21cd0c07 |
226 | Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ; |
9a1398dd |
227 | return kFALSE; |
228 | } |
229 | } |
230 | |
231 | Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly |
232 | // depends on it. |
233 | Double_t p1 = 1.0 ; |
234 | Double_t p2 = 0.0 ; |
235 | |
236 | gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls |
237 | gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient |
238 | gMinuit->SetMaxIterations(5); |
239 | gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings |
240 | |
241 | gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize |
242 | |
243 | if(ierflg == 4){ // Minimum not found |
21cd0c07 |
244 | Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" ); |
9a1398dd |
245 | return kFALSE ; |
246 | } |
247 | for(index = 0; index < nPar; index++){ |
248 | Double_t err ; |
249 | Double_t val ; |
250 | gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index |
251 | fitparameters[index] = val ; |
252 | } |
253 | |
254 | delete toMinuit ; |
255 | return kTRUE; |
c3c187e5 |
256 | |
d15a28e7 |
257 | } |
258 | |
3758d9fc |
259 | //____________________________________________________________________________ |
260 | void AliPHOSClusterizerv1::GetCalibrationParameters() |
261 | { |
fbf811ec |
262 | |
88cb7938 |
263 | AliPHOSGetter * gime = AliPHOSGetter::Instance(); |
c2cd5471 |
264 | |
88cb7938 |
265 | if ( !gime->Digitizer() ) |
266 | gime->LoadDigitizer(); |
267 | AliPHOSDigitizer * dig = gime->Digitizer(); |
268 | |
269 | fADCchanelEmc = dig->GetEMCchannel() ; |
270 | fADCpedestalEmc = dig->GetEMCpedestal(); |
271 | |
272 | fADCchanelCpv = dig->GetCPVchannel() ; |
273 | fADCpedestalCpv = dig->GetCPVpedestal() ; |
274 | |
275 | // fCalibrationDB = gime->CalibrationDB(); |
3758d9fc |
276 | } |
548f0134 |
277 | |
d15a28e7 |
278 | //____________________________________________________________________________ |
a4e98857 |
279 | void AliPHOSClusterizerv1::Init() |
280 | { |
2bd5457f |
281 | // Make all memory allocations which can not be done in default constructor. |
282 | // Attach the Clusterizer task to the list of PHOS tasks |
88cb7938 |
283 | |
284 | AliPHOSGetter* gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data()); |
fbf811ec |
285 | |
88cb7938 |
286 | AliPHOSGeometry * geom = gime->PHOSGeometry(); |
287 | |
3758d9fc |
288 | fEmcCrystals = geom->GetNModules() * geom->GetNCristalsInModule() ; |
289 | |
7b7c1533 |
290 | if(!gMinuit) |
88cb7938 |
291 | gMinuit = new TMinuit(100); |
fbf811ec |
292 | |
88cb7938 |
293 | if ( !gime->Clusterizer() ) |
294 | gime->PostClusterizer(this); |
9a1398dd |
295 | } |
7b7c1533 |
296 | |
8d0f3f77 |
297 | //____________________________________________________________________________ |
298 | void AliPHOSClusterizerv1::InitParameters() |
299 | { |
300 | |
301 | fNumberOfCpvClusters = 0 ; |
302 | fNumberOfEmcClusters = 0 ; |
303 | |
304 | fCpvClusteringThreshold = 0.0; |
305 | fEmcClusteringThreshold = 0.2; |
306 | |
307 | fEmcLocMaxCut = 0.03 ; |
308 | fCpvLocMaxCut = 0.03 ; |
309 | |
310 | fW0 = 4.5 ; |
311 | fW0CPV = 4.0 ; |
312 | |
313 | fEmcTimeGate = 1.e-8 ; |
314 | |
315 | fToUnfold = kTRUE ; |
88cb7938 |
316 | |
8d0f3f77 |
317 | fRecPointsInRun = 0 ; |
c2cd5471 |
318 | |
8d0f3f77 |
319 | } |
320 | |
9a1398dd |
321 | //____________________________________________________________________________ |
322 | Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const |
d15a28e7 |
323 | { |
b2a60966 |
324 | // Gives the neighbourness of two digits = 0 are not neighbour but continue searching |
325 | // = 1 are neighbour |
326 | // = 2 are not neighbour but do not continue searching |
a4e98857 |
327 | // neighbours are defined as digits having at least a common vertex |
b2a60966 |
328 | // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster |
329 | // which is compared to a digit (d2) not yet in a cluster |
330 | |
88cb7938 |
331 | AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ; |
7b7c1533 |
332 | |
d15a28e7 |
333 | Int_t rv = 0 ; |
334 | |
d15a28e7 |
335 | Int_t relid1[4] ; |
7b7c1533 |
336 | geom->AbsToRelNumbering(d1->GetId(), relid1) ; |
d15a28e7 |
337 | |
338 | Int_t relid2[4] ; |
7b7c1533 |
339 | geom->AbsToRelNumbering(d2->GetId(), relid2) ; |
d15a28e7 |
340 | |
9688c1dd |
341 | if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module |
92862013 |
342 | Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ; |
343 | Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ; |
d15a28e7 |
344 | |
92862013 |
345 | if (( coldiff <= 1 ) && ( rowdiff <= 1 )){ |
9688c1dd |
346 | if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fEmcTimeGate)) |
d15a28e7 |
347 | rv = 1 ; |
348 | } |
349 | else { |
350 | if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1)) |
88cb7938 |
351 | rv = 2; // Difference in row numbers is too large to look further |
d15a28e7 |
352 | } |
353 | |
354 | } |
355 | else { |
356 | |
9688c1dd |
357 | if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) ) |
d15a28e7 |
358 | rv=2 ; |
359 | |
360 | } |
d72dfbc3 |
361 | |
d15a28e7 |
362 | return rv ; |
363 | } |
364 | |
d15a28e7 |
365 | |
366 | //____________________________________________________________________________ |
9a1398dd |
367 | Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const |
d15a28e7 |
368 | { |
b2a60966 |
369 | // Tells if (true) or not (false) the digit is in a PHOS-EMC module |
370 | |
9f616d61 |
371 | Bool_t rv = kFALSE ; |
88cb7938 |
372 | AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ; |
d15a28e7 |
373 | |
2b629790 |
374 | Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ(); |
d15a28e7 |
375 | |
2b629790 |
376 | if(digit->GetId() <= nEMC ) rv = kTRUE; |
ed4205d8 |
377 | |
378 | return rv ; |
379 | } |
380 | |
ed4205d8 |
381 | //____________________________________________________________________________ |
9a1398dd |
382 | Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const |
ed4205d8 |
383 | { |
fad3e5b9 |
384 | // Tells if (true) or not (false) the digit is in a PHOS-CPV module |
ed4205d8 |
385 | |
386 | Bool_t rv = kFALSE ; |
88cb7938 |
387 | |
388 | AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ; |
ed4205d8 |
389 | |
2b629790 |
390 | Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ(); |
ed4205d8 |
391 | |
2b629790 |
392 | if(digit->GetId() > nEMC ) rv = kTRUE; |
d15a28e7 |
393 | |
394 | return rv ; |
395 | } |
9a1398dd |
396 | |
397 | //____________________________________________________________________________ |
88cb7938 |
398 | void AliPHOSClusterizerv1::Unload() |
399 | { |
400 | AliPHOSGetter * gime = AliPHOSGetter::Instance() ; |
401 | gime->PhosLoader()->UnloadDigits() ; |
402 | gime->PhosLoader()->UnloadRecPoints() ; |
403 | } |
404 | |
405 | //____________________________________________________________________________ |
406 | void AliPHOSClusterizerv1::WriteRecPoints() |
a4e98857 |
407 | { |
7b7c1533 |
408 | |
409 | // Creates new branches with given title |
a4e98857 |
410 | // fills and writes into TreeR. |
88cb7938 |
411 | |
412 | AliPHOSGetter * gime = AliPHOSGetter::Instance(); |
413 | |
fbf811ec |
414 | TObjArray * emcRecPoints = gime->EmcRecPoints() ; |
415 | TObjArray * cpvRecPoints = gime->CpvRecPoints() ; |
416 | TClonesArray * digits = gime->Digits() ; |
fbf811ec |
417 | |
88cb7938 |
418 | TTree * treeR = gime->TreeR(); |
fbf811ec |
419 | |
9a1398dd |
420 | Int_t index ; |
01a599c9 |
421 | //Evaluate position, dispersion and other RecPoint properties... |
88cb7938 |
422 | for(index = 0; index < emcRecPoints->GetEntries(); index++) |
423 | dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) )->EvalAll(fW0,digits) ; |
424 | |
7b7c1533 |
425 | emcRecPoints->Sort() ; |
7b7c1533 |
426 | for(index = 0; index < emcRecPoints->GetEntries(); index++) |
88cb7938 |
427 | dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) )->SetIndexInList(index) ; |
fbf811ec |
428 | |
7b7c1533 |
429 | emcRecPoints->Expand(emcRecPoints->GetEntriesFast()) ; |
fbf811ec |
430 | |
9a1398dd |
431 | //Now the same for CPV |
7b7c1533 |
432 | for(index = 0; index < cpvRecPoints->GetEntries(); index++) |
e957fea8 |
433 | dynamic_cast<AliPHOSRecPoint *>( cpvRecPoints->At(index) )->EvalAll(digits) ; |
fbf811ec |
434 | |
7b7c1533 |
435 | cpvRecPoints->Sort() ; |
fbf811ec |
436 | |
7b7c1533 |
437 | for(index = 0; index < cpvRecPoints->GetEntries(); index++) |
88cb7938 |
438 | dynamic_cast<AliPHOSRecPoint *>( cpvRecPoints->At(index) )->SetIndexInList(index) ; |
fbf811ec |
439 | |
7b7c1533 |
440 | cpvRecPoints->Expand(cpvRecPoints->GetEntriesFast()) ; |
9a1398dd |
441 | |
9a1398dd |
442 | Int_t bufferSize = 32000 ; |
443 | Int_t splitlevel = 0 ; |
8d0f3f77 |
444 | |
01a599c9 |
445 | //First EMC |
8d0f3f77 |
446 | TBranch * emcBranch = treeR->Branch("PHOSEmcRP","TObjArray",&emcRecPoints,bufferSize,splitlevel); |
fc12304f |
447 | emcBranch->SetTitle(BranchName()); |
9a1398dd |
448 | |
449 | //Now CPV branch |
8d0f3f77 |
450 | TBranch * cpvBranch = treeR->Branch("PHOSCpvRP","TObjArray",&cpvRecPoints,bufferSize,splitlevel); |
fc12304f |
451 | cpvBranch->SetTitle(BranchName()); |
9a1398dd |
452 | |
01a599c9 |
453 | emcBranch ->Fill() ; |
454 | cpvBranch ->Fill() ; |
88cb7938 |
455 | |
456 | gime->WriteRecPoints("OVERWRITE"); |
457 | gime->WriteClusterizer("OVERWRITE"); |
9a1398dd |
458 | } |
459 | |
460 | //____________________________________________________________________________ |
461 | void AliPHOSClusterizerv1::MakeClusters() |
462 | { |
463 | // Steering method to construct the clusters stored in a list of Reconstructed Points |
464 | // A cluster is defined as a list of neighbour digits |
d15a28e7 |
465 | |
88cb7938 |
466 | |
467 | AliPHOSGetter * gime = AliPHOSGetter::Instance(); |
468 | |
469 | TObjArray * emcRecPoints = gime->EmcRecPoints() ; |
470 | TObjArray * cpvRecPoints = gime->CpvRecPoints() ; |
7b7c1533 |
471 | |
7d493c2c |
472 | emcRecPoints->Delete() ; |
473 | cpvRecPoints->Delete() ; |
7b7c1533 |
474 | |
fbf811ec |
475 | TClonesArray * digits = gime->Digits() ; |
88cb7938 |
476 | |
477 | TClonesArray * digitsC = static_cast<TClonesArray*>( digits->Clone() ) ; |
7b7c1533 |
478 | |
479 | |
d15a28e7 |
480 | // Clusterization starts |
9a1398dd |
481 | |
7b7c1533 |
482 | TIter nextdigit(digitsC) ; |
d15a28e7 |
483 | AliPHOSDigit * digit ; |
92862013 |
484 | Bool_t notremoved = kTRUE ; |
6ad0bfa0 |
485 | |
88cb7938 |
486 | while ( (digit = dynamic_cast<AliPHOSDigit *>( nextdigit()) ) ) { // scan over the list of digitsC |
487 | |
488 | |
9a1398dd |
489 | AliPHOSRecPoint * clu = 0 ; |
c161df70 |
490 | |
d1de15f5 |
491 | TArrayI clusterdigitslist(1500) ; |
d15a28e7 |
492 | Int_t index ; |
f2bc1b87 |
493 | |
3758d9fc |
494 | if (( IsInEmc (digit) && Calibrate(digit->GetAmp(),digit->GetId()) > fEmcClusteringThreshold ) || |
495 | ( IsInCpv (digit) && Calibrate(digit->GetAmp(),digit->GetId()) > fCpvClusteringThreshold ) ) { |
d15a28e7 |
496 | Int_t iDigitInCluster = 0 ; |
7b7c1533 |
497 | |
6ad0bfa0 |
498 | if ( IsInEmc(digit) ) { |
88cb7938 |
499 | // start a new EMC RecPoint |
500 | if(fNumberOfEmcClusters >= emcRecPoints->GetSize()) |
501 | emcRecPoints->Expand(2*fNumberOfEmcClusters+1) ; |
502 | |
503 | emcRecPoints->AddAt(new AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ; |
504 | clu = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(fNumberOfEmcClusters) ) ; |
505 | fNumberOfEmcClusters++ ; |
506 | clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->GetId())) ; |
507 | clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ; |
508 | iDigitInCluster++ ; |
509 | digitsC->Remove(digit) ; |
ca77b032 |
510 | |
d72dfbc3 |
511 | } else { |
88cb7938 |
512 | |
513 | // start a new CPV cluster |
514 | if(fNumberOfCpvClusters >= cpvRecPoints->GetSize()) |
515 | cpvRecPoints->Expand(2*fNumberOfCpvClusters+1); |
516 | |
517 | cpvRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ; |
518 | |
519 | clu = dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ; |
520 | fNumberOfCpvClusters++ ; |
521 | clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->GetId()) ) ; |
522 | clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ; |
523 | iDigitInCluster++ ; |
524 | digitsC->Remove(digit) ; |
d15a28e7 |
525 | nextdigit.Reset() ; |
88cb7938 |
526 | |
527 | // Here we remove remaining EMC digits, which cannot make a cluster |
528 | |
92862013 |
529 | if( notremoved ) { |
88cb7938 |
530 | while( ( digit = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) { |
9f616d61 |
531 | if( IsInEmc(digit) ) |
88cb7938 |
532 | digitsC->Remove(digit) ; |
d15a28e7 |
533 | else |
88cb7938 |
534 | break ; |
535 | } |
536 | notremoved = kFALSE ; |
537 | } |
538 | |
d15a28e7 |
539 | } // else |
540 | |
541 | nextdigit.Reset() ; |
fad3e5b9 |
542 | |
d15a28e7 |
543 | AliPHOSDigit * digitN ; |
544 | index = 0 ; |
9f616d61 |
545 | while (index < iDigitInCluster){ // scan over digits already in cluster |
88cb7938 |
546 | digit = dynamic_cast<AliPHOSDigit*>( digits->At(clusterdigitslist[index]) ) ; |
547 | index++ ; |
548 | while ( (digitN = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) { // scan over the reduced list of digits |
549 | Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!! |
d15a28e7 |
550 | switch (ineb ) { |
9f616d61 |
551 | case 0 : // not a neighbour |
88cb7938 |
552 | break ; |
553 | case 1 : // are neighbours |
554 | clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), digitN->GetId() ) ) ; |
555 | clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ; |
556 | iDigitInCluster++ ; |
557 | digitsC->Remove(digitN) ; |
558 | break ; |
9f616d61 |
559 | case 2 : // too far from each other |
88cb7938 |
560 | goto endofloop; |
561 | } // switch |
562 | |
563 | } // while digitN |
564 | |
d15a28e7 |
565 | endofloop: ; |
88cb7938 |
566 | nextdigit.Reset() ; |
567 | |
d15a28e7 |
568 | } // loop over cluster |
ad8cfaf4 |
569 | |
ad8cfaf4 |
570 | } // energy theshold |
ad8cfaf4 |
571 | |
31aa6d6c |
572 | |
d15a28e7 |
573 | } // while digit |
574 | |
7b7c1533 |
575 | delete digitsC ; |
9688c1dd |
576 | |
9a1398dd |
577 | } |
578 | |
579 | //____________________________________________________________________________ |
a4e98857 |
580 | void AliPHOSClusterizerv1::MakeUnfolding() |
581 | { |
582 | // Unfolds clusters using the shape of an ElectroMagnetic shower |
9688c1dd |
583 | // Performs unfolding of all EMC/CPV clusters |
9a1398dd |
584 | |
88cb7938 |
585 | AliPHOSGetter * gime = AliPHOSGetter::Instance(); |
fc12304f |
586 | |
7b7c1533 |
587 | const AliPHOSGeometry * geom = gime->PHOSGeometry() ; |
88cb7938 |
588 | |
589 | TObjArray * emcRecPoints = gime->EmcRecPoints() ; |
590 | TObjArray * cpvRecPoints = gime->CpvRecPoints() ; |
fbf811ec |
591 | TClonesArray * digits = gime->Digits() ; |
7b7c1533 |
592 | |
a4e98857 |
593 | // Unfold first EMC clusters |
9a1398dd |
594 | if(fNumberOfEmcClusters > 0){ |
595 | |
7b7c1533 |
596 | Int_t nModulesToUnfold = geom->GetNModules() ; |
9a1398dd |
597 | |
598 | Int_t numberofNotUnfolded = fNumberOfEmcClusters ; |
599 | Int_t index ; |
600 | for(index = 0 ; index < numberofNotUnfolded ; index++){ |
601 | |
88cb7938 |
602 | AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) ) ; |
9a1398dd |
603 | if(emcRecPoint->GetPHOSMod()> nModulesToUnfold) |
88cb7938 |
604 | break ; |
9a1398dd |
605 | |
606 | Int_t nMultipl = emcRecPoint->GetMultiplicity() ; |
a0636361 |
607 | AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ; |
9a1398dd |
608 | Float_t * maxAtEnergy = new Float_t[nMultipl] ; |
7b7c1533 |
609 | Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,digits) ; |
9a1398dd |
610 | |
611 | if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0 |
88cb7938 |
612 | UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ; |
613 | emcRecPoints->Remove(emcRecPoint); |
614 | emcRecPoints->Compress() ; |
615 | index-- ; |
616 | fNumberOfEmcClusters -- ; |
617 | numberofNotUnfolded-- ; |
9a1398dd |
618 | } |
619 | |
620 | delete[] maxAt ; |
621 | delete[] maxAtEnergy ; |
622 | } |
623 | } |
a4e98857 |
624 | // Unfolding of EMC clusters finished |
9a1398dd |
625 | |
626 | |
a4e98857 |
627 | // Unfold now CPV clusters |
9a1398dd |
628 | if(fNumberOfCpvClusters > 0){ |
629 | |
9688c1dd |
630 | Int_t nModulesToUnfold = geom->GetNModules() ; |
9a1398dd |
631 | |
632 | Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ; |
633 | Int_t index ; |
634 | for(index = 0 ; index < numberofCpvNotUnfolded ; index++){ |
635 | |
88cb7938 |
636 | AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( cpvRecPoints->At(index) ) ; |
9a1398dd |
637 | |
638 | if(recPoint->GetPHOSMod()> nModulesToUnfold) |
88cb7938 |
639 | break ; |
9a1398dd |
640 | |
88cb7938 |
641 | AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ; |
9a1398dd |
642 | |
643 | Int_t nMultipl = emcRecPoint->GetMultiplicity() ; |
a0636361 |
644 | AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ; |
9a1398dd |
645 | Float_t * maxAtEnergy = new Float_t[nMultipl] ; |
7b7c1533 |
646 | Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,digits) ; |
9a1398dd |
647 | |
648 | if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0 |
88cb7938 |
649 | UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ; |
650 | cpvRecPoints->Remove(emcRecPoint); |
651 | cpvRecPoints->Compress() ; |
652 | index-- ; |
653 | numberofCpvNotUnfolded-- ; |
654 | fNumberOfCpvClusters-- ; |
9a1398dd |
655 | } |
656 | |
657 | delete[] maxAt ; |
658 | delete[] maxAtEnergy ; |
659 | } |
660 | } |
661 | //Unfolding of Cpv clusters finished |
662 | |
663 | } |
664 | |
9a1398dd |
665 | //____________________________________________________________________________ |
666 | Double_t AliPHOSClusterizerv1::ShowerShape(Double_t r) |
667 | { |
668 | // Shape of the shower (see PHOS TDR) |
a4e98857 |
669 | // If you change this function, change also the gradient evaluation in ChiSquare() |
9a1398dd |
670 | |
671 | Double_t r4 = r*r*r*r ; |
672 | Double_t r295 = TMath::Power(r, 2.95) ; |
673 | Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ; |
674 | return shape ; |
675 | } |
676 | |
677 | //____________________________________________________________________________ |
678 | void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc, |
88cb7938 |
679 | Int_t nMax, |
680 | AliPHOSDigit ** maxAt, |
681 | Float_t * maxAtEnergy) |
9a1398dd |
682 | { |
683 | // Performs the unfolding of a cluster with nMax overlapping showers |
684 | |
88cb7938 |
685 | AliPHOSGetter * gime = AliPHOSGetter::Instance(); |
fc12304f |
686 | |
7b7c1533 |
687 | const AliPHOSGeometry * geom = gime->PHOSGeometry() ; |
88cb7938 |
688 | |
fbf811ec |
689 | const TClonesArray * digits = gime->Digits() ; |
690 | TObjArray * emcRecPoints = gime->EmcRecPoints() ; |
691 | TObjArray * cpvRecPoints = gime->CpvRecPoints() ; |
7b7c1533 |
692 | |
9a1398dd |
693 | Int_t nPar = 3 * nMax ; |
694 | Float_t * fitparameters = new Float_t[nPar] ; |
695 | |
696 | Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ; |
697 | if( !rv ) { |
88cb7938 |
698 | // Fit failed, return and remove cluster |
9a1398dd |
699 | delete[] fitparameters ; |
700 | return ; |
701 | } |
702 | |
703 | // create ufolded rec points and fill them with new energy lists |
704 | // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[] |
705 | // and later correct this number in acordance with actual energy deposition |
706 | |
707 | Int_t nDigits = iniEmc->GetMultiplicity() ; |
708 | Float_t * efit = new Float_t[nDigits] ; |
7b7c1533 |
709 | Float_t xDigit=0.,zDigit=0.,distance=0. ; |
710 | Float_t xpar=0.,zpar=0.,epar=0. ; |
9a1398dd |
711 | Int_t relid[4] ; |
7b7c1533 |
712 | AliPHOSDigit * digit = 0 ; |
9a1398dd |
713 | Int_t * emcDigits = iniEmc->GetDigitsList() ; |
714 | |
715 | Int_t iparam ; |
716 | Int_t iDigit ; |
717 | for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){ |
88cb7938 |
718 | digit = dynamic_cast<AliPHOSDigit*>( digits->At(emcDigits[iDigit] ) ) ; |
7b7c1533 |
719 | geom->AbsToRelNumbering(digit->GetId(), relid) ; |
720 | geom->RelPosInModule(relid, xDigit, zDigit) ; |
9a1398dd |
721 | efit[iDigit] = 0; |
722 | |
723 | iparam = 0 ; |
724 | while(iparam < nPar ){ |
725 | xpar = fitparameters[iparam] ; |
726 | zpar = fitparameters[iparam+1] ; |
727 | epar = fitparameters[iparam+2] ; |
728 | iparam += 3 ; |
729 | distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ; |
730 | distance = TMath::Sqrt(distance) ; |
731 | efit[iDigit] += epar * ShowerShape(distance) ; |
732 | } |
733 | } |
734 | |
735 | |
736 | // Now create new RecPoints and fill energy lists with efit corrected to fluctuations |
737 | // so that energy deposited in each cell is distributed betwin new clusters proportionally |
738 | // to its contribution to efit |
739 | |
740 | Float_t * emcEnergies = iniEmc->GetEnergiesList() ; |
741 | Float_t ratio ; |
742 | |
743 | iparam = 0 ; |
744 | while(iparam < nPar ){ |
745 | xpar = fitparameters[iparam] ; |
746 | zpar = fitparameters[iparam+1] ; |
747 | epar = fitparameters[iparam+2] ; |
748 | iparam += 3 ; |
749 | |
7b7c1533 |
750 | AliPHOSEmcRecPoint * emcRP = 0 ; |
9a1398dd |
751 | |
752 | if(iniEmc->IsEmc()){ //create new entries in fEmcRecPoints... |
753 | |
7b7c1533 |
754 | if(fNumberOfEmcClusters >= emcRecPoints->GetSize()) |
88cb7938 |
755 | emcRecPoints->Expand(2*fNumberOfEmcClusters) ; |
9a1398dd |
756 | |
73a68ccb |
757 | (*emcRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ; |
88cb7938 |
758 | emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(fNumberOfEmcClusters) ) ; |
9a1398dd |
759 | fNumberOfEmcClusters++ ; |
760 | } |
761 | else{//create new entries in fCpvRecPoints |
7b7c1533 |
762 | if(fNumberOfCpvClusters >= cpvRecPoints->GetSize()) |
88cb7938 |
763 | cpvRecPoints->Expand(2*fNumberOfCpvClusters) ; |
9a1398dd |
764 | |
73a68ccb |
765 | (*cpvRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ; |
88cb7938 |
766 | emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ; |
9a1398dd |
767 | fNumberOfCpvClusters++ ; |
768 | } |
769 | |
770 | Float_t eDigit ; |
771 | for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){ |
88cb7938 |
772 | digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) ) ; |
7b7c1533 |
773 | geom->AbsToRelNumbering(digit->GetId(), relid) ; |
774 | geom->RelPosInModule(relid, xDigit, zDigit) ; |
9a1398dd |
775 | distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ; |
776 | distance = TMath::Sqrt(distance) ; |
777 | ratio = epar * ShowerShape(distance) / efit[iDigit] ; |
778 | eDigit = emcEnergies[iDigit] * ratio ; |
779 | emcRP->AddDigit( *digit, eDigit ) ; |
88cb7938 |
780 | } |
9a1398dd |
781 | } |
782 | |
783 | delete[] fitparameters ; |
784 | delete[] efit ; |
785 | |
786 | } |
787 | |
788 | //_____________________________________________________________________________ |
789 | void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag) |
790 | { |
a4e98857 |
791 | // Calculates the Chi square for the cluster unfolding minimization |
9a1398dd |
792 | // Number of parameters, Gradient, Chi squared, parameters, what to do |
793 | |
88cb7938 |
794 | TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ; |
9a1398dd |
795 | |
88cb7938 |
796 | AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) ) ; |
797 | TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) ) ; |
9a1398dd |
798 | |
799 | |
800 | |
88cb7938 |
801 | // AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit |
9a1398dd |
802 | |
803 | Int_t * emcDigits = emcRP->GetDigitsList() ; |
804 | |
7b7c1533 |
805 | Int_t nOdigits = emcRP->GetDigitsMultiplicity() ; |
9a1398dd |
806 | |
807 | Float_t * emcEnergies = emcRP->GetEnergiesList() ; |
88cb7938 |
808 | |
809 | const AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ; |
9a1398dd |
810 | fret = 0. ; |
811 | Int_t iparam ; |
812 | |
813 | if(iflag == 2) |
814 | for(iparam = 0 ; iparam < nPar ; iparam++) |
815 | Grad[iparam] = 0 ; // Will evaluate gradient |
816 | |
817 | Double_t efit ; |
818 | |
819 | AliPHOSDigit * digit ; |
820 | Int_t iDigit ; |
821 | |
7b7c1533 |
822 | for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) { |
9a1398dd |
823 | |
88cb7938 |
824 | digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) ); |
9a1398dd |
825 | |
826 | Int_t relid[4] ; |
827 | Float_t xDigit ; |
828 | Float_t zDigit ; |
829 | |
830 | geom->AbsToRelNumbering(digit->GetId(), relid) ; |
831 | |
832 | geom->RelPosInModule(relid, xDigit, zDigit) ; |
833 | |
834 | if(iflag == 2){ // calculate gradient |
835 | Int_t iParam = 0 ; |
836 | efit = 0 ; |
837 | while(iParam < nPar ){ |
88cb7938 |
838 | Double_t distance = (xDigit - x[iParam]) * (xDigit - x[iParam]) ; |
839 | iParam++ ; |
840 | distance += (zDigit - x[iParam]) * (zDigit - x[iParam]) ; |
841 | distance = TMath::Sqrt( distance ) ; |
842 | iParam++ ; |
843 | efit += x[iParam] * ShowerShape(distance) ; |
844 | iParam++ ; |
9a1398dd |
845 | } |
88cb7938 |
846 | Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E) |
9a1398dd |
847 | iParam = 0 ; |
848 | while(iParam < nPar ){ |
88cb7938 |
849 | Double_t xpar = x[iParam] ; |
850 | Double_t zpar = x[iParam+1] ; |
851 | Double_t epar = x[iParam+2] ; |
852 | Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ); |
853 | Double_t shape = sum * ShowerShape(dr) ; |
854 | Double_t r4 = dr*dr*dr*dr ; |
855 | Double_t r295 = TMath::Power(dr,2.95) ; |
856 | Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) + |
857 | 0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ; |
858 | |
859 | Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x |
860 | iParam++ ; |
861 | Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z |
862 | iParam++ ; |
863 | Grad[iParam] += shape ; // Derivative over energy |
864 | iParam++ ; |
9a1398dd |
865 | } |
866 | } |
867 | efit = 0; |
868 | iparam = 0 ; |
869 | |
870 | while(iparam < nPar ){ |
871 | Double_t xpar = x[iparam] ; |
872 | Double_t zpar = x[iparam+1] ; |
873 | Double_t epar = x[iparam+2] ; |
874 | iparam += 3 ; |
875 | Double_t distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ; |
876 | distance = TMath::Sqrt(distance) ; |
877 | efit += epar * ShowerShape(distance) ; |
878 | } |
879 | |
88cb7938 |
880 | fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ; |
9a1398dd |
881 | // Here we assume, that sigma = sqrt(E) |
882 | } |
83974468 |
883 | |
d15a28e7 |
884 | } |
885 | |
886 | //____________________________________________________________________________ |
88cb7938 |
887 | void AliPHOSClusterizerv1::Print()const |
d15a28e7 |
888 | { |
d1de15f5 |
889 | // Print clusterizer parameters |
890 | |
21cd0c07 |
891 | TString message ; |
892 | TString taskName(GetName()) ; |
893 | taskName.ReplaceAll(Version(), "") ; |
894 | |
895 | if( strcmp(GetName(), "") !=0 ) { |
9a1398dd |
896 | // Print parameters |
21cd0c07 |
897 | message = "\n--------------- %s %s -----------\n" ; |
898 | message += "Clusterizing digits from the file: %s\n" ; |
899 | message += " Branch: %s\n" ; |
900 | message += " EMC Clustering threshold = %f\n" ; |
901 | message += " EMC Local Maximum cut = %f\n" ; |
902 | message += " EMC Logarothmic weight = %f\n" ; |
903 | message += " CPV Clustering threshold = %f\n" ; |
904 | message += " CPV Local Maximum cut = %f\n" ; |
905 | message += " CPV Logarothmic weight = %f\n" ; |
9a1398dd |
906 | if(fToUnfold) |
21cd0c07 |
907 | message += " Unfolding on\n" ; |
9a1398dd |
908 | else |
21cd0c07 |
909 | message += " Unfolding off\n" ; |
9a1398dd |
910 | |
21cd0c07 |
911 | message += "------------------------------------------------------------------" ; |
9a1398dd |
912 | } |
21cd0c07 |
913 | else |
914 | message = " AliPHOSClusterizerv1 not initialized " ; |
915 | |
916 | Info("Print", message.Data(), |
917 | taskName.Data(), |
918 | GetTitle(), |
919 | taskName.Data(), |
920 | GetName(), |
921 | fEmcClusteringThreshold, |
922 | fEmcLocMaxCut, |
923 | fW0, |
924 | fCpvClusteringThreshold, |
925 | fCpvLocMaxCut, |
926 | fW0CPV ) ; |
9a1398dd |
927 | } |
21cd0c07 |
928 | |
929 | |
9a1398dd |
930 | //____________________________________________________________________________ |
a4e98857 |
931 | void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option) |
932 | { |
933 | // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer |
9a1398dd |
934 | |
88cb7938 |
935 | AliPHOSGetter * gime = AliPHOSGetter::Instance(); |
936 | |
937 | TObjArray * emcRecPoints = gime->EmcRecPoints() ; |
938 | TObjArray * cpvRecPoints = gime->CpvRecPoints() ; |
7b7c1533 |
939 | |
709e117a |
940 | printf("\nevent %d \n", gAlice->GetEvNumber()) ; |
941 | printf(" Found %d EMC RecPoints and %d CPV RecPoints \n", |
942 | emcRecPoints->GetEntriesFast(),cpvRecPoints->GetEntriesFast()) ; |
88cb7938 |
943 | |
944 | fRecPointsInRun += emcRecPoints->GetEntriesFast() ; |
945 | fRecPointsInRun += cpvRecPoints->GetEntriesFast() ; |
946 | |
11f9c5ff |
947 | |
9a1398dd |
948 | if(strstr(option,"all")) { |
709e117a |
949 | printf("\n EMC clusters \n") ; |
950 | printf("Index Ene(MeV) Multi Module X Y Z Lambdas_1 Lambda_2 # of prim Primaries list\n") ; |
9a1398dd |
951 | Int_t index ; |
7b7c1533 |
952 | for (index = 0 ; index < emcRecPoints->GetEntries() ; index++) { |
953 | AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )emcRecPoints->At(index) ; |
9a1398dd |
954 | TVector3 locpos; |
955 | rp->GetLocalPosition(locpos); |
9a1398dd |
956 | Float_t lambda[2]; |
957 | rp->GetElipsAxis(lambda); |
9a1398dd |
958 | Int_t * primaries; |
959 | Int_t nprimaries; |
960 | primaries = rp->GetPrimaries(nprimaries); |
709e117a |
961 | printf("\n%6d %8.2f %3d %2d %4.1f %4.1f %4.1f %4f %4f %2d : ", |
11f9c5ff |
962 | rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(), |
963 | locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ; |
3bf72d32 |
964 | |
21cd0c07 |
965 | for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) { |
709e117a |
966 | printf("%d ", primaries[iprimary] ) ; |
21cd0c07 |
967 | } |
709e117a |
968 | printf("\n") ; |
969 | } |
970 | |
971 | //Now plot CPV recPoints |
972 | printf("\n CPV clusters \n") ; |
973 | printf("Index Ene(MeV) Module X Y Z \n") ; |
974 | for (index = 0 ; index < cpvRecPoints->GetEntries() ; index++) { |
975 | AliPHOSRecPoint * rp = (AliPHOSRecPoint * )cpvRecPoints->At(index) ; |
11f9c5ff |
976 | |
709e117a |
977 | TVector3 locpos; |
978 | rp->GetLocalPosition(locpos); |
11f9c5ff |
979 | |
709e117a |
980 | printf("\n%6d %8.2f %2d %4.1f %4.1f %4.1f \n", |
11f9c5ff |
981 | rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(), |
709e117a |
982 | locpos.X(), locpos.Y(), locpos.Z()) ; |
88cb7938 |
983 | } |
9a1398dd |
984 | } |
d15a28e7 |
985 | } |
9a1398dd |
986 | |