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 : /* History of cvs commits:
19 : *
20 : * $Log: AliPHOSClusterizerv1.cxx,v $
21 : * Revision 1.118 2007/12/11 21:23:26 kharlov
22 : * Added possibility to swith off unfolding
23 : *
24 : * Revision 1.117 2007/10/18 08:42:05 kharlov
25 : * Bad channels cleaned before clusterization
26 : *
27 : * Revision 1.116 2007/10/01 20:24:08 kharlov
28 : * Memory leaks fixed
29 : *
30 : * Revision 1.115 2007/09/26 14:22:17 cvetan
31 : * Important changes to the reconstructor classes. Complete elimination of the run-loaders, which are now steered only from AliReconstruction. Removal of the corresponding Reconstruct() and FillESD() methods.
32 : *
33 : * Revision 1.114 2007/09/06 16:06:44 kharlov
34 : * Absence of sorting results in loose of all unfolded clusters
35 : *
36 : * Revision 1.113 2007/08/28 12:55:07 policheh
37 : * Loaders removed from the reconstruction code (C.Cheshkov)
38 : *
39 : * Revision 1.112 2007/08/22 09:20:50 hristov
40 : * Updated QA classes (Yves)
41 : *
42 : * Revision 1.111 2007/08/08 12:11:28 kharlov
43 : * Protection against uninitialized fQADM
44 : *
45 : * Revision 1.110 2007/08/07 14:16:00 kharlov
46 : * Quality assurance added (Yves Schutz)
47 : *
48 : * Revision 1.109 2007/07/24 17:20:35 policheh
49 : * Usage of RecoParam objects instead of hardcoded parameters in reconstruction.
50 : * (See $ALICE_ROOT/PHOS/macros/BeamTest2006/RawReconstruction.C).
51 : *
52 : * Revision 1.108 2007/06/18 07:00:51 kharlov
53 : * Bug fix for attempt to use AliPHOSEmcRecPoint after its deletion
54 : *
55 : * Revision 1.107 2007/05/25 14:12:26 policheh
56 : * Local to tracking CS transformation added for CPV rec. points
57 : *
58 : * Revision 1.106 2007/05/24 13:01:22 policheh
59 : * Local to tracking CS transformation invoked for each EMC rec.point
60 : *
61 : * Revision 1.105 2007/05/02 13:41:22 kharlov
62 : * Mode protection against absence of calib.data from AliPHOSCalibData to AliPHOSClusterizerv1::GetCalibrationParameters()
63 : *
64 : * Revision 1.104 2007/04/27 16:55:53 kharlov
65 : * Calibration stops if PHOS CDB objects do not exist
66 : *
67 : * Revision 1.103 2007/04/11 11:55:45 policheh
68 : * SetDistancesToBadChannels() added.
69 : *
70 : * Revision 1.102 2007/03/28 19:18:15 kharlov
71 : * RecPoints recalculation in TSM removed
72 : *
73 : * Revision 1.101 2007/03/06 06:51:27 kharlov
74 : * Calculation of cluster properties dep. on vertex posponed to TrackSegmentMaker
75 : *
76 : * Revision 1.100 2007/01/10 11:07:26 kharlov
77 : * Raw digits writing to file (B.Polichtchouk)
78 : *
79 : * Revision 1.99 2006/11/07 16:49:51 kharlov
80 : * Corrections for next event switch in case of raw data (B.Polichtchouk)
81 : *
82 : * Revision 1.98 2006/10/27 17:14:27 kharlov
83 : * Introduce AliDebug and AliLog (B.Polichtchouk)
84 : *
85 : * Revision 1.97 2006/08/29 11:41:19 kharlov
86 : * Missing implementation of ctors and = operator are added
87 : *
88 : * Revision 1.96 2006/08/25 16:56:30 kharlov
89 : * Compliance with Effective C++
90 : *
91 : * Revision 1.95 2006/08/11 12:36:26 cvetan
92 : * Update of the PHOS code needed in order to read and reconstruct the beam test raw data (i.e. without an existing galice.root)
93 : *
94 : * Revision 1.94 2006/08/07 12:27:49 hristov
95 : * Removing obsolete code which affected the event numbering scheme
96 : *
97 : * Revision 1.93 2006/08/01 12:20:17 cvetan
98 : * 1. Adding a possibility to read and reconstruct an old rcu formatted raw data. This is controlled by an option of AliReconstruction and AliPHOSReconstructor. 2. In case of raw data processing (without galice.root) create the default AliPHOSGeometry object. Most likely this should be moved to the CDB
99 : *
100 : * Revision 1.92 2006/04/29 20:26:46 hristov
101 : * Separate EMC and CPV calibration (Yu.Kharlov)
102 : *
103 : * Revision 1.91 2006/04/22 10:30:17 hristov
104 : * Add fEnergy to AliPHOSDigit and operate with EMC amplitude in energy units (Yu.Kharlov)
105 : *
106 : * Revision 1.90 2006/04/11 15:22:59 hristov
107 : * run number in query set to -1: forces AliCDBManager to use its run number (A.Colla)
108 : *
109 : * Revision 1.89 2006/03/13 14:05:42 kharlov
110 : * Calibration objects for EMC and CPV
111 : *
112 : * Revision 1.88 2006/01/11 08:54:52 hristov
113 : * Additional protection in case no calibration entry was found
114 : *
115 : * Revision 1.87 2005/11/22 08:46:43 kharlov
116 : * Updated to new CDB (Boris Polichtchouk)
117 : *
118 : * Revision 1.86 2005/11/14 21:52:43 hristov
119 : * Coding conventions
120 : *
121 : * Revision 1.85 2005/09/27 16:08:08 hristov
122 : * New version of CDB storage framework (A.Colla)
123 : *
124 : * Revision 1.84 2005/09/21 10:02:47 kharlov
125 : * Reading calibration from CDB (Boris Polichtchouk)
126 : *
127 : * Revision 1.82 2005/09/02 15:43:13 kharlov
128 : * Add comments in GetCalibrationParameters and Calibrate
129 : *
130 : * Revision 1.81 2005/09/02 14:32:07 kharlov
131 : * Calibration of raw data
132 : *
133 : * Revision 1.80 2005/08/24 15:31:36 kharlov
134 : * Setting raw digits flag
135 : *
136 : * Revision 1.79 2005/07/25 15:53:53 kharlov
137 : * Read raw data
138 : *
139 : * Revision 1.78 2005/05/28 14:19:04 schutz
140 : * Compilation warnings fixed by T.P.
141 : *
142 : */
143 :
144 : //*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
145 : //////////////////////////////////////////////////////////////////////////////
146 : // Clusterization class. Performs clusterization (collects neighbouring active cells) and
147 : // unfolds the clusters having several local maxima.
148 : // Results are stored in TreeR#, branches PHOSEmcRP (EMC recPoints),
149 : // PHOSCpvRP (CPV RecPoints) and AliPHOSClusterizer (Clusterizer with all
150 : // parameters including input digits branch title, thresholds etc.)
151 : // This TTask is normally called from Reconstructioner, but can as well be used in
152 : // standalone mode.
153 : // Use Case:
154 : // root [0] AliPHOSClusterizerv1 * cl = new AliPHOSClusterizerv1(<pointer_to_phos_geometry_onject>)
155 : // root [1] cl->Digits2Clusters(digitsTree,clusterTree)
156 : // //finds RecPoints in the current event
157 : // root [2] cl->SetDigitsBranch("digits2")
158 : // //sets another title for Digitis (input) branch
159 : // root [3] cl->SetRecPointsBranch("recp2")
160 : // //sets another title four output branches
161 : // root [4] cl->SetEmcLocalMaxCut(0.03)
162 : // //set clusterization parameters
163 :
164 : // --- ROOT system ---
165 :
166 : #include "TMath.h"
167 : #include "TMinuit.h"
168 : #include "TTree.h"
169 : #include "TBenchmark.h"
170 : #include "TClonesArray.h"
171 :
172 : // --- Standard library ---
173 :
174 : // --- AliRoot header files ---
175 : #include "AliLog.h"
176 : #include "AliConfig.h"
177 : #include "AliPHOSGeometry.h"
178 : #include "AliPHOSClusterizerv1.h"
179 : #include "AliPHOSEmcRecPoint.h"
180 : #include "AliPHOSCpvRecPoint.h"
181 : #include "AliPHOSDigit.h"
182 : #include "AliPHOSDigitizer.h"
183 : #include "AliCDBManager.h"
184 : #include "AliCDBStorage.h"
185 : #include "AliCDBEntry.h"
186 : #include "AliPHOSRecoParam.h"
187 : #include "AliPHOSReconstructor.h"
188 : #include "AliPHOSCalibData.h"
189 :
190 22 : ClassImp(AliPHOSClusterizerv1)
191 :
192 : //____________________________________________________________________________
193 : AliPHOSClusterizerv1::AliPHOSClusterizerv1() :
194 0 : AliPHOSClusterizer(),
195 0 : fDefaultInit(0), fEmcCrystals(0), fToUnfold(0), fToUnfoldCPV(0),
196 0 : fWrite(0),
197 0 : fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
198 0 : fEmcClusteringThreshold(0), fCpvClusteringThreshold(0),
199 0 : fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
200 0 : fW0CPV(0),
201 0 : fTimeGateLowAmp(0.), fTimeGateLow(0.), fTimeGateHigh(0.),
202 0 : fEcoreRadius(0)
203 0 : {
204 : // default ctor (to be used mainly by Streamer)
205 :
206 0 : fDefaultInit = kTRUE ;
207 :
208 0 : for(Int_t i=0; i<53760; i++){
209 0 : fDigitsUsed[i]=0 ;
210 : }
211 0 : }
212 :
213 : //____________________________________________________________________________
214 : AliPHOSClusterizerv1::AliPHOSClusterizerv1(AliPHOSGeometry *geom) :
215 2 : AliPHOSClusterizer(geom),
216 2 : fDefaultInit(0), fEmcCrystals(0), fToUnfold(0), fToUnfoldCPV(0),
217 2 : fWrite(0),
218 2 : fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
219 2 : fEmcClusteringThreshold(0), fCpvClusteringThreshold(0),
220 2 : fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
221 2 : fW0CPV(0),
222 2 : fTimeGateLowAmp(0.), fTimeGateLow(0.), fTimeGateHigh(0.),
223 2 : fEcoreRadius(0)
224 10 : {
225 : // ctor with the indication of the file where header Tree and digits Tree are stored
226 :
227 215044 : for(Int_t i=0; i<53760; i++){
228 107520 : fDigitsUsed[i]=0 ;
229 : }
230 :
231 2 : Init() ;
232 2 : fDefaultInit = kFALSE ;
233 4 : }
234 :
235 : //____________________________________________________________________________
236 : AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
237 8 : {
238 : // dtor
239 :
240 8 : }
241 : //____________________________________________________________________________
242 : void AliPHOSClusterizerv1::Digits2Clusters(Option_t *option)
243 : {
244 : // Steering method to perform clusterization for one event
245 : // The input is the tree with digits.
246 : // The output is the tree with clusters.
247 :
248 16 : if(strstr(option,"tim"))
249 0 : gBenchmark->Start("PHOSClusterizer");
250 :
251 8 : if(strstr(option,"print")) {
252 0 : Print() ;
253 0 : return ;
254 : }
255 :
256 8 : MakeClusters() ;
257 :
258 24 : AliDebug(2,Form("Number of EMC clusters: %d, CPV clusters: %d\n",
259 : fEMCRecPoints->GetEntriesFast(), fCPVRecPoints->GetEntriesFast()));
260 8 : if(AliLog::GetGlobalDebugLevel()>1) {
261 0 : fEMCRecPoints->Print();
262 0 : fCPVRecPoints->Print();
263 0 : }
264 :
265 8 : MakeUnfolding();
266 :
267 8 : WriteRecPoints();
268 :
269 8 : if(strstr(option,"deb"))
270 0 : PrintRecPoints(option) ;
271 :
272 8 : if(strstr(option,"tim")){
273 0 : gBenchmark->Stop("PHOSClusterizer");
274 0 : AliInfo(Form("took %f seconds for Clusterizing\n",
275 : gBenchmark->GetCpuTime("PHOSClusterizer")));
276 0 : }
277 8 : fEMCRecPoints->Delete();
278 8 : fCPVRecPoints->Delete();
279 16 : }
280 :
281 : //____________________________________________________________________________
282 : Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** maxAt, Float_t * maxAtEnergy,
283 : Int_t nPar, Float_t * fitparameters) const
284 : {
285 : // Calls TMinuit to fit the energy distribution of a cluster with several maxima
286 : // The initial values for fitting procedure are set equal to the positions of local maxima.
287 : // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
288 :
289 :
290 0 : if(!gMinuit) //it was deleted by someone else
291 0 : gMinuit = new TMinuit(100) ;
292 0 : gMinuit->mncler(); // Reset Minuit's list of paramters
293 0 : gMinuit->SetPrintLevel(-1) ; // No Printout
294 0 : gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;
295 : // To set the address of the minimization function
296 :
297 0 : TList * toMinuit = new TList();
298 0 : toMinuit->AddAt(emcRP,0) ;
299 0 : toMinuit->AddAt(fDigitsArr,1) ;
300 0 : toMinuit->AddAt(fGeom,2) ;
301 :
302 0 : gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare
303 :
304 : // filling initial values for fit parameters
305 : AliPHOSDigit * digit ;
306 :
307 0 : Int_t ierflg = 0;
308 : Int_t index = 0 ;
309 0 : Int_t nDigits = (Int_t) nPar / 3 ;
310 :
311 : Int_t iDigit ;
312 :
313 0 : for(iDigit = 0; iDigit < nDigits; iDigit++){
314 0 : digit = maxAt[iDigit];
315 :
316 0 : Int_t relid[4] ;
317 0 : Float_t x = 0.;
318 0 : Float_t z = 0.;
319 0 : fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
320 0 : fGeom->RelPosInModule(relid, x, z) ;
321 :
322 0 : Float_t energy = maxAtEnergy[iDigit] ;
323 :
324 0 : gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ;
325 0 : index++ ;
326 0 : if(ierflg != 0){
327 0 : Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ;
328 0 : return kFALSE;
329 : }
330 0 : gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ;
331 0 : index++ ;
332 0 : if(ierflg != 0){
333 0 : Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ;
334 0 : return kFALSE;
335 : }
336 0 : gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;
337 0 : index++ ;
338 0 : if(ierflg != 0){
339 0 : Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ;
340 0 : return kFALSE;
341 : }
342 0 : }
343 :
344 0 : Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
345 : // depends on it.
346 0 : Double_t p1 = 1.0 ;
347 0 : Double_t p2 = 0.0 ;
348 :
349 0 : gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
350 0 : gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
351 0 : gMinuit->SetMaxIterations(5);
352 0 : gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
353 :
354 0 : gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
355 :
356 0 : if(ierflg == 4){ // Minimum not found
357 0 : Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" );
358 0 : return kFALSE ;
359 : }
360 0 : for(index = 0; index < nPar; index++){
361 0 : Double_t err ;
362 0 : Double_t val ;
363 0 : gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index
364 0 : fitparameters[index] = val ;
365 0 : }
366 :
367 0 : delete toMinuit ;
368 0 : return kTRUE;
369 :
370 0 : }
371 :
372 :
373 : //____________________________________________________________________________
374 : void AliPHOSClusterizerv1::Init()
375 : {
376 : // Make all memory allocations which can not be done in default constructor.
377 : // Attach the Clusterizer task to the list of PHOS tasks
378 :
379 4 : fEmcCrystals = fGeom->GetNModules() * fGeom->GetNCristalsInModule() ;
380 :
381 2 : if(!gMinuit)
382 4 : gMinuit = new TMinuit(100);
383 :
384 2 : if (!fgCalibData)
385 4 : fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
386 2 : if (fgCalibData->GetCalibDataEmc() == 0)
387 0 : AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n");
388 2 : if (fgCalibData->GetCalibDataCpv() == 0)
389 0 : AliFatal("Calibration parameters for PHOS CPV not found. Stop reconstruction.\n");
390 :
391 2 : }
392 :
393 : //____________________________________________________________________________
394 : void AliPHOSClusterizerv1::InitParameters()
395 : {
396 :
397 16 : fNumberOfCpvClusters = 0 ;
398 8 : fNumberOfEmcClusters = 0 ;
399 :
400 8 : const AliPHOSRecoParam* recoParam = AliPHOSReconstructor::GetRecoParam();
401 8 : if(!recoParam) AliFatal("Reconstruction parameters are not set!");
402 :
403 8 : recoParam->Print();
404 :
405 8 : fEmcClusteringThreshold = recoParam->GetEMCClusteringThreshold();
406 8 : fCpvClusteringThreshold = recoParam->GetCPVClusteringThreshold();
407 :
408 8 : fEmcLocMaxCut = recoParam->GetEMCLocalMaxCut();
409 8 : fCpvLocMaxCut = recoParam->GetCPVLocalMaxCut();
410 :
411 8 : fW0 = recoParam->GetEMCLogWeight();
412 8 : fW0CPV = recoParam->GetCPVLogWeight();
413 :
414 8 : fTimeGateLowAmp = recoParam->GetTimeGateAmpThresh() ;
415 8 : fTimeGateLow = recoParam->GetTimeGateLow() ;
416 8 : fTimeGateHigh = recoParam->GetTimeGateHigh() ;
417 :
418 8 : fEcoreRadius = recoParam->GetEMCEcoreRadius();
419 :
420 8 : fToUnfold = recoParam->EMCToUnfold() ;
421 8 : fToUnfoldCPV = recoParam->CPVToUnfold() ;
422 :
423 8 : fWrite = kTRUE ;
424 8 : }
425 :
426 : //____________________________________________________________________________
427 : Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
428 : {
429 : // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
430 : // = 1 are neighbour
431 : // = 2 are not neighbour but do not continue searching
432 : // =-1 are not neighbour, continue searching, but do not look before d2 next time
433 : // neighbours are defined as digits having at least a common vertex
434 : // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
435 : // which is compared to a digit (d2) not yet in a cluster
436 :
437 2150 : Int_t relid1[4] ;
438 1075 : fGeom->AbsToRelNumbering(d1->GetId(), relid1) ;
439 :
440 1075 : Int_t relid2[4] ;
441 1075 : fGeom->AbsToRelNumbering(d2->GetId(), relid2) ;
442 :
443 1963 : if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
444 888 : Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
445 888 : Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
446 :
447 888 : if (( coldiff <= 1 ) && ( rowdiff <= 1 )){ //At least common vertex
448 : // if (( relid1[2]==relid2[2] && coldiff <= 1 ) || ( relid1[3]==relid2[3] && rowdiff <= 1 )){ //common side
449 332 : if((relid1[1] != 0) || CheckTimeGate(d1->GetTime(),d1->GetEnergy(),d2->GetTime(),d2->GetEnergy()))
450 166 : return 1 ;
451 : }
452 : else {
453 1104 : if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
454 119 : return 2; // Difference in row numbers is too large to look further
455 : }
456 603 : return 0 ;
457 :
458 : }
459 : else {
460 374 : if(relid1[0] > relid2[0] && relid1[1]==relid2[1] ) //we switched to the next module
461 187 : return -1 ;
462 0 : if(relid1[1] < relid2[1]) //we switched from EMC(0) to CPV(-1)
463 0 : return -1 ;
464 :
465 0 : return 2 ;
466 :
467 : }
468 :
469 : return 0 ;
470 1075 : }
471 : //____________________________________________________________________________
472 : Bool_t AliPHOSClusterizerv1::CheckTimeGate(Float_t t1, Float_t amp1, Float_t t2, Float_t amp2)const{
473 : //Check if two cells have reasonable time difference
474 : //Note that at low amplitude time is defined up to 1 tick == 100 ns.
475 478 : if(amp1<fTimeGateLowAmp || amp2<fTimeGateLowAmp){
476 88 : return (TMath::Abs(t1 - t2 ) < fTimeGateLow) ;
477 : }
478 : else{ //Time should be measured with good accuracy
479 78 : return (TMath::Abs(t1 - t2 ) < fTimeGateHigh) ;
480 : }
481 :
482 166 : }
483 : //____________________________________________________________________________
484 : Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
485 : {
486 : // Tells if (true) or not (false) the digit is in a PHOS-EMC module
487 :
488 : Bool_t rv = kFALSE ;
489 :
490 284 : Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
491 :
492 284 : if(digit->GetId() <= nEMC ) rv = kTRUE;
493 :
494 142 : return rv ;
495 : }
496 :
497 : //____________________________________________________________________________
498 : Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
499 : {
500 : // Tells if (true) or not (false) the digit is in a PHOS-CPV module
501 :
502 : Bool_t rv = kFALSE ;
503 :
504 244 : Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
505 :
506 122 : if(digit->GetId() > nEMC ) rv = kTRUE;
507 :
508 122 : return rv ;
509 : }
510 :
511 : //____________________________________________________________________________
512 : void AliPHOSClusterizerv1::WriteRecPoints()
513 : {
514 :
515 : // Creates new branches with given title
516 : // fills and writes into TreeR.
517 :
518 : Int_t index ;
519 : //Evaluate position, dispersion and other RecPoint properties..
520 16 : Int_t nEmc = fEMCRecPoints->GetEntriesFast();
521 8 : Float_t emcMinE= AliPHOSReconstructor::GetRecoParam()->GetEMCMinE(); //Minimal digit energy
522 8 : TVector3 fakeVtx(0.,0.,0.) ;
523 36 : for(index = 0; index < nEmc; index++){
524 : AliPHOSEmcRecPoint * rp =
525 20 : static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) );
526 10 : rp->Purify(emcMinE,fDigitsArr) ;
527 10 : if(rp->GetMultiplicity()==0){
528 0 : fEMCRecPoints->RemoveAt(index) ;
529 0 : delete rp ;
530 0 : continue;
531 : }
532 :
533 : // No vertex is available now, calculate corrections in PID
534 10 : rp->EvalAll(fDigitsArr) ;
535 10 : rp->EvalCoreEnergy(fW0,fEcoreRadius,fDigitsArr) ;
536 10 : rp->EvalAll(fW0,fakeVtx,fDigitsArr) ;
537 10 : rp->EvalLocal2TrackingCSTransform();
538 10 : }
539 8 : fEMCRecPoints->Compress() ;
540 8 : fEMCRecPoints->Sort() ;
541 : // fEMCRecPoints->Expand(fEMCRecPoints->GetEntriesFast()) ;
542 54 : for(index = 0; index < fEMCRecPoints->GetEntries(); index++){
543 20 : static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) )->SetIndexInList(index) ;
544 : }
545 :
546 : //For each rec.point set the distance to the nearest bad crystal (BVP)
547 8 : SetDistancesToBadChannels();
548 :
549 : //Now the same for CPV
550 16 : Float_t cpvMinE= AliPHOSReconstructor::GetRecoParam()->GetCPVMinE(); //Minimal digit energy
551 24 : for(index = 0; index < fCPVRecPoints->GetEntries(); index++){
552 0 : AliPHOSCpvRecPoint * rp = static_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) );
553 0 : rp->Purify(cpvMinE,fDigitsArr) ;
554 0 : if(rp->GetMultiplicity()==0){
555 0 : fCPVRecPoints->RemoveAt(index) ;
556 0 : delete rp ;
557 0 : continue;
558 : }
559 :
560 0 : rp->EvalAll(fDigitsArr) ;
561 0 : rp->EvalAll(fW0CPV,fakeVtx,fDigitsArr) ;
562 0 : rp->EvalLocal2TrackingCSTransform();
563 0 : }
564 8 : fCPVRecPoints->Sort() ;
565 :
566 32 : for(index = 0; index < fCPVRecPoints->GetEntries(); index++)
567 8 : static_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) )->SetIndexInList(index) ;
568 :
569 16 : fCPVRecPoints->Expand(fCPVRecPoints->GetEntriesFast()) ;
570 :
571 8 : if(fWrite){ //We write TreeR
572 8 : fTreeR->Fill();
573 : }
574 8 : }
575 :
576 : //____________________________________________________________________________
577 : void AliPHOSClusterizerv1::MakeClusters()
578 : {
579 : // Steering method to construct the clusters stored in a list of Reconstructed Points
580 : // A cluster is defined as a list of neighbour digits
581 :
582 16 : fNumberOfCpvClusters = 0 ;
583 8 : fNumberOfEmcClusters = 0 ;
584 :
585 : //Mark all digits as unused yet
586 : const Int_t maxNDigits = 3584; // There is no clusters larger than PHOS module ;)
587 8 : Int_t nDigits=fDigitsArr->GetEntriesFast() ;
588 :
589 480 : for(Int_t i=0; i<nDigits; i++){
590 232 : fDigitsUsed[i]=0 ;
591 : }
592 : Int_t iFirst = 0 ; //first index of digit which potentially can be a part of cluster
593 : //e.g. first digit in this module, first CPV digit etc.
594 : AliPHOSDigit * digit ;
595 8 : TArrayI clusterdigitslist(maxNDigits) ;
596 : AliPHOSRecPoint * clu = 0 ;
597 488 : for(Int_t i=0; i<nDigits; i++){
598 232 : if(fDigitsUsed[i])
599 : continue ;
600 :
601 264 : digit=static_cast<AliPHOSDigit*>(fDigitsArr->At(i)) ;
602 :
603 : clu=0 ;
604 :
605 : Int_t index ;
606 :
607 : //is this digit so energetic that start cluster?
608 660 : AliDebug(2,Form("Digit %d, energy=%f, ID=%d",i,digit->GetEnergy(),digit->GetId()));
609 528 : if (( IsInEmc(digit) && Calibrate(digit->GetEnergy(),digit->GetId()) > fEmcClusteringThreshold ) ||
610 244 : ( IsInCpv(digit) && Calibrate(digit->GetEnergy(),digit->GetId()) > fCpvClusteringThreshold ) ) {
611 : Int_t iDigitInCluster = 0 ;
612 20 : if ( IsInEmc(digit) ) {
613 : // start a new EMC RecPoint
614 20 : if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
615 0 : fEMCRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
616 :
617 30 : fEMCRecPoints->AddAt(new AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
618 20 : clu = static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
619 10 : fNumberOfEmcClusters++ ;
620 30 : clu->AddDigit(*digit, Calibrate(digit->GetEnergy(),digit->GetId()),CalibrateT(digit->GetTime(),digit->GetId(),digit->IsLG())) ;
621 20 : clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
622 : iDigitInCluster++ ;
623 10 : fDigitsUsed[i]=kTRUE ;
624 10 : } else {
625 : // start a new CPV cluster
626 0 : if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
627 0 : fCPVRecPoints->Expand(2*fNumberOfCpvClusters+1);
628 :
629 0 : fCPVRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
630 0 : clu = static_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
631 0 : fNumberOfCpvClusters++ ;
632 0 : clu->AddDigit(*digit, Calibrate(digit->GetEnergy(),digit->GetId()),0.) ; // no timing information in CPV
633 0 : clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
634 : iDigitInCluster++ ;
635 0 : fDigitsUsed[i]=kTRUE ;
636 : } // else
637 :
638 : //Now scan remaining digits in list to find neigbours of our seed
639 :
640 : AliPHOSDigit * digitN ;
641 : index = 0 ;
642 196 : while (index < iDigitInCluster){ // scan over digits already in cluster
643 528 : digit = static_cast<AliPHOSDigit*>( fDigitsArr->At(clusterdigitslist[index]) ) ;
644 176 : index++ ;
645 8089 : for(Int_t j=iFirst; j<nDigits; j++){
646 3959 : if (iDigitInCluster >= maxNDigits) {
647 0 : AliError(Form("The number of digits in cluster is more than %d, skip the rest of event",
648 : maxNDigits));
649 0 : return;
650 : }
651 3959 : if(fDigitsUsed[j])
652 : continue ; //look through remaining digits
653 2150 : digitN = static_cast<AliPHOSDigit*>( fDigitsArr->At(j) ) ;
654 1075 : Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
655 1262 : switch (ineb ) {
656 : case -1: //too early (e.g. previous module), do not look before j at subsequent passes
657 : iFirst=j ;
658 187 : break ;
659 : case 0 : // not a neighbour
660 : break ;
661 : case 1 : // are neighbours
662 498 : clu->AddDigit(*digitN, Calibrate(digitN->GetEnergy(),digitN->GetId()),CalibrateT(digitN->GetTime(),digitN->GetId(),digit->IsLG())) ;
663 332 : clusterdigitslist[iDigitInCluster] = j ;
664 166 : iDigitInCluster++ ;
665 166 : fDigitsUsed[j]=kTRUE ;
666 166 : break ;
667 : case 2 : // too far from each other
668 119 : goto endOfLoop;
669 : } // switch
670 :
671 956 : }
672 :
673 : endOfLoop: ; //scanned all possible neighbours for this digit
674 :
675 : } // loop over cluster
676 10 : } // energy theshold
677 132 : }
678 :
679 16 : }
680 :
681 : //____________________________________________________________________________
682 : void AliPHOSClusterizerv1::MakeUnfolding()
683 : {
684 : // Unfolds clusters using the shape of an ElectroMagnetic shower
685 : // Performs unfolding of all EMC/CPV clusters
686 :
687 24 : if(fToUnfold && (fNumberOfEmcClusters > 0)){ // Unfold first EMC clusters
688 :
689 8 : Int_t nModulesToUnfold = fGeom->GetNModules() ;
690 :
691 8 : Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
692 : Int_t index ;
693 36 : for(index = 0 ; index < numberofNotUnfolded ; index++){
694 :
695 10 : AliPHOSEmcRecPoint * emcRecPoint = static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) ) ;
696 10 : if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
697 0 : break ;
698 :
699 10 : Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
700 10 : AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
701 10 : Float_t * maxAtEnergy = new Float_t[nMultipl] ;
702 10 : Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,fDigitsArr) ;
703 :
704 10 : if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
705 0 : UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
706 :
707 0 : fEMCRecPoints->Remove(emcRecPoint);
708 0 : fEMCRecPoints->Compress() ;
709 0 : index-- ;
710 0 : fNumberOfEmcClusters -- ;
711 0 : numberofNotUnfolded-- ;
712 0 : }
713 : else{
714 10 : emcRecPoint->SetNExMax(1) ; //Only one local maximum
715 : }
716 :
717 20 : delete[] maxAt ;
718 20 : delete[] maxAtEnergy ;
719 10 : }
720 8 : }
721 : // Unfolding of EMC clusters finished
722 :
723 :
724 : // Unfold now CPV clusters
725 16 : if(fToUnfoldCPV && (fNumberOfCpvClusters > 0)){
726 :
727 0 : Int_t nModulesToUnfold = fGeom->GetNModules() ;
728 :
729 0 : Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
730 : Int_t index ;
731 0 : for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
732 :
733 0 : AliPHOSRecPoint * recPoint = static_cast<AliPHOSRecPoint *>( fCPVRecPoints->At(index) ) ;
734 :
735 0 : if(recPoint->GetPHOSMod()> nModulesToUnfold)
736 0 : break ;
737 :
738 0 : AliPHOSEmcRecPoint * emcRecPoint = static_cast<AliPHOSEmcRecPoint*>(recPoint) ;
739 :
740 0 : Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
741 0 : AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
742 0 : Float_t * maxAtEnergy = new Float_t[nMultipl] ;
743 0 : Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,fDigitsArr) ;
744 :
745 : //Number of points to fit should be larger than number of parameters
746 0 : if( nMax > 1 && emcRecPoint->GetMultiplicity()>3*nMax ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
747 0 : UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
748 0 : fCPVRecPoints->Remove(emcRecPoint);
749 0 : fCPVRecPoints->Compress() ;
750 0 : index-- ;
751 0 : numberofCpvNotUnfolded-- ;
752 0 : fNumberOfCpvClusters-- ;
753 0 : }
754 :
755 0 : delete[] maxAt ;
756 0 : delete[] maxAtEnergy ;
757 0 : }
758 0 : }
759 : //Unfolding of Cpv clusters finished
760 :
761 8 : }
762 :
763 : //____________________________________________________________________________
764 : Double_t AliPHOSClusterizerv1::ShowerShape(Double_t x, Double_t z)
765 : {
766 : // Shape of the shower (see PHOS TDR)
767 : // If you change this function, change also the gradient evaluation in ChiSquare()
768 :
769 : //for the moment we neglect dependence on the incident angle.
770 :
771 0 : Double_t r2 = x*x + z*z ;
772 0 : Double_t r4 = r2*r2 ;
773 0 : Double_t r295 = TMath::Power(r2, 2.95/2.) ;
774 0 : Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
775 0 : return shape ;
776 : }
777 :
778 : //____________________________________________________________________________
779 : void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
780 : Int_t nMax,
781 : AliPHOSDigit ** maxAt,
782 : Float_t * maxAtEnergy)
783 : {
784 : // Performs the unfolding of a cluster with nMax overlapping showers
785 :
786 0 : Int_t nPar = 3 * nMax ;
787 0 : Float_t * fitparameters = new Float_t[nPar] ;
788 :
789 0 : Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
790 :
791 0 : if( !rv ) {
792 : // Fit failed, return and remove cluster
793 0 : iniEmc->SetNExMax(-1) ;
794 0 : delete[] fitparameters ;
795 0 : return ;
796 : }
797 :
798 : // create ufolded rec points and fill them with new energy lists
799 : // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
800 : // and later correct this number in acordance with actual energy deposition
801 :
802 0 : Int_t nDigits = iniEmc->GetMultiplicity() ;
803 0 : Float_t * efit = new Float_t[nDigits] ;
804 0 : Float_t xDigit=0.,zDigit=0. ;
805 : Float_t xpar=0.,zpar=0.,epar=0. ;
806 0 : Int_t relid[4] ;
807 : AliPHOSDigit * digit = 0 ;
808 0 : Int_t * emcDigits = iniEmc->GetDigitsList() ;
809 :
810 0 : TVector3 vIncid ;
811 :
812 : Int_t iparam ;
813 : Int_t iDigit ;
814 0 : for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
815 0 : digit = static_cast<AliPHOSDigit*>( fDigitsArr->At(emcDigits[iDigit] ) ) ;
816 0 : fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
817 0 : fGeom->RelPosInModule(relid, xDigit, zDigit) ;
818 0 : efit[iDigit] = 0;
819 :
820 : iparam = 0 ;
821 0 : while(iparam < nPar ){
822 0 : xpar = fitparameters[iparam] ;
823 0 : zpar = fitparameters[iparam+1] ;
824 0 : epar = fitparameters[iparam+2] ;
825 0 : iparam += 3 ;
826 : // fGeom->GetIncidentVector(fVtx,relid[0],xpar,zpar,vIncid) ;
827 : // efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) ;
828 0 : efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
829 : }
830 : }
831 :
832 : // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
833 : // so that energy deposited in each cell is distributed betwin new clusters proportionally
834 : // to its contribution to efit
835 :
836 0 : Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
837 : Float_t ratio ;
838 :
839 : iparam = 0 ;
840 0 : while(iparam < nPar ){
841 0 : xpar = fitparameters[iparam] ;
842 0 : zpar = fitparameters[iparam+1] ;
843 0 : epar = fitparameters[iparam+2] ;
844 0 : iparam += 3 ;
845 : // fGeom->GetIncidentVector(fVtx,iniEmc->GetPHOSMod(),xpar,zpar,vIncid) ;
846 :
847 : AliPHOSEmcRecPoint * emcRP = 0 ;
848 :
849 0 : if(iniEmc->IsEmc()){ //create new entries in fEMCRecPoints...
850 :
851 0 : if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
852 0 : fEMCRecPoints->Expand(2*fNumberOfEmcClusters) ;
853 :
854 0 : (*fEMCRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
855 0 : emcRP = static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
856 0 : fNumberOfEmcClusters++ ;
857 0 : emcRP->SetNExMax((Int_t)nPar/3) ;
858 0 : }
859 : else{//create new entries in fCPVRecPoints
860 0 : if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
861 0 : fCPVRecPoints->Expand(2*fNumberOfCpvClusters) ;
862 :
863 0 : (*fCPVRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
864 0 : emcRP = static_cast<AliPHOSEmcRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
865 0 : fNumberOfCpvClusters++ ;
866 : }
867 :
868 : Float_t eDigit ;
869 0 : for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
870 0 : digit = static_cast<AliPHOSDigit*>( fDigitsArr->At( emcDigits[iDigit] ) ) ;
871 0 : fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
872 0 : fGeom->RelPosInModule(relid, xDigit, zDigit) ;
873 0 : if(efit[iDigit]>0){
874 : // ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) / efit[iDigit] ;
875 0 : ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
876 0 : eDigit = emcEnergies[iDigit] * ratio ;
877 0 : emcRP->AddDigit( *digit, eDigit,CalibrateT(digit->GetTime(),digit->GetId(),digit->IsLG()) ) ;
878 : }
879 : }
880 : }
881 :
882 0 : delete[] fitparameters ;
883 0 : delete[] efit ;
884 :
885 0 : }
886 :
887 : //_____________________________________________________________________________
888 : void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
889 : {
890 : // Calculates the Chi square for the cluster unfolding minimization
891 : // Number of parameters, Gradient, Chi squared, parameters, what to do
892 :
893 0 : TList * toMinuit = static_cast<TList*>( gMinuit->GetObjectFit() ) ;
894 :
895 0 : AliPHOSEmcRecPoint * emcRP = static_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) ) ;
896 0 : TClonesArray * digits = static_cast<TClonesArray*>( toMinuit->At(1) ) ;
897 : // A bit buggy way to get an access to the geometry
898 : // To be revised!
899 0 : AliPHOSGeometry *geom = static_cast<AliPHOSGeometry *>(toMinuit->At(2));
900 :
901 : // TVector3 * vtx = static_cast<TVector3*>(toMinuit->At(3)) ; //Vertex position
902 :
903 : // AliPHOSEmcRecPoint * emcRP = static_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
904 :
905 0 : Int_t * emcDigits = emcRP->GetDigitsList() ;
906 :
907 0 : Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
908 :
909 0 : Float_t * emcEnergies = emcRP->GetEnergiesList() ;
910 :
911 : // TVector3 vInc ;
912 0 : fret = 0. ;
913 : Int_t iparam ;
914 :
915 0 : if(iflag == 2)
916 0 : for(iparam = 0 ; iparam < nPar ; iparam++)
917 0 : Grad[iparam] = 0 ; // Will evaluate gradient
918 :
919 : Double_t efit ;
920 :
921 : AliPHOSDigit * digit ;
922 : Int_t iDigit ;
923 :
924 0 : for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
925 :
926 0 : digit = static_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) );
927 :
928 0 : Int_t relid[4] ;
929 0 : Float_t xDigit ;
930 0 : Float_t zDigit ;
931 :
932 0 : geom->AbsToRelNumbering(digit->GetId(), relid) ;
933 :
934 0 : geom->RelPosInModule(relid, xDigit, zDigit) ;
935 :
936 0 : if(iflag == 2){ // calculate gradient
937 : Int_t iParam = 0 ;
938 : efit = 0 ;
939 0 : while(iParam < nPar ){
940 0 : Double_t dx = (xDigit - x[iParam]) ;
941 0 : iParam++ ;
942 0 : Double_t dz = (zDigit - x[iParam]) ;
943 0 : iParam++ ;
944 : // fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),x[iParam-2],x[iParam-1],vInc) ;
945 : // efit += x[iParam] * ShowerShape(dx,dz,vInc) ;
946 0 : efit += x[iParam] * ShowerShape(dx,dz) ;
947 0 : iParam++ ;
948 : }
949 0 : Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
950 : iParam = 0 ;
951 0 : while(iParam < nPar ){
952 0 : Double_t xpar = x[iParam] ;
953 0 : Double_t zpar = x[iParam+1] ;
954 0 : Double_t epar = x[iParam+2] ;
955 : // fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
956 0 : Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
957 : // Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
958 0 : Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ;
959 : //DP: No incident angle dependence in gradient yet!!!!!!
960 0 : Double_t r4 = dr*dr*dr*dr ;
961 0 : Double_t r295 = TMath::Power(dr,2.95) ;
962 0 : Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
963 0 : 0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
964 :
965 0 : Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x
966 : iParam++ ;
967 0 : Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z
968 0 : iParam++ ;
969 0 : Grad[iParam] += shape ; // Derivative over energy
970 0 : iParam++ ;
971 : }
972 0 : }
973 : efit = 0;
974 : iparam = 0 ;
975 :
976 0 : while(iparam < nPar ){
977 0 : Double_t xpar = x[iparam] ;
978 0 : Double_t zpar = x[iparam+1] ;
979 0 : Double_t epar = x[iparam+2] ;
980 0 : iparam += 3 ;
981 : // fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
982 : // efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
983 0 : efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
984 : }
985 :
986 0 : fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
987 : // Here we assume, that sigma = sqrt(E)
988 0 : }
989 :
990 0 : }
991 :
992 : //____________________________________________________________________________
993 : void AliPHOSClusterizerv1::Print(const Option_t *)const
994 : {
995 : // Print clusterizer parameters
996 :
997 0 : TString message ;
998 0 : TString taskName(GetName()) ;
999 0 : taskName.ReplaceAll(Version(), "") ;
1000 :
1001 0 : if( strcmp(GetName(), "") !=0 ) {
1002 : // Print parameters
1003 0 : message = "\n--------------- %s %s -----------\n" ;
1004 0 : message += "Clusterizing digits from the file: %s\n" ;
1005 0 : message += " Branch: %s\n" ;
1006 0 : message += " EMC Clustering threshold = %f\n" ;
1007 0 : message += " EMC Local Maximum cut = %f\n" ;
1008 0 : message += " EMC Logarothmic weight = %f\n" ;
1009 0 : message += " CPV Clustering threshold = %f\n" ;
1010 0 : message += " CPV Local Maximum cut = %f\n" ;
1011 0 : message += " CPV Logarothmic weight = %f\n" ;
1012 0 : if(fToUnfold)
1013 0 : message += " Unfolding on\n" ;
1014 : else
1015 0 : message += " Unfolding off\n" ;
1016 :
1017 0 : message += "------------------------------------------------------------------" ;
1018 : }
1019 : else
1020 0 : message = " AliPHOSClusterizerv1 not initialized " ;
1021 :
1022 0 : AliInfo(Form("%s, %s %s %s %s %f %f %f %f %f %f", message.Data(),
1023 : taskName.Data(),
1024 : GetTitle(),
1025 : taskName.Data(),
1026 : GetName(),
1027 : fEmcClusteringThreshold,
1028 : fEmcLocMaxCut,
1029 : fW0,
1030 : fCpvClusteringThreshold,
1031 : fCpvLocMaxCut,
1032 : fW0CPV )) ;
1033 0 : }
1034 : //____________________________________________________________________________
1035 : void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
1036 : {
1037 : // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
1038 :
1039 0 : AliInfo(Form("\nFound %d EMC RecPoints and %d CPV RecPoints",
1040 : fEMCRecPoints->GetEntriesFast(),
1041 : fCPVRecPoints->GetEntriesFast() )) ;
1042 :
1043 0 : if(strstr(option,"all")) {
1044 0 : printf("\n EMC clusters \n") ;
1045 0 : printf("Index Ene(MeV) Multi Module X Y Z Lambdas_1 Lambda_2 # of prim Primaries list\n") ;
1046 : Int_t index ;
1047 0 : for (index = 0 ; index < fEMCRecPoints->GetEntries() ; index++) {
1048 0 : AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )fEMCRecPoints->At(index) ;
1049 0 : TVector3 locpos;
1050 0 : rp->GetLocalPosition(locpos);
1051 0 : Float_t lambda[2];
1052 0 : rp->GetElipsAxis(lambda);
1053 : Int_t * primaries;
1054 0 : Int_t nprimaries;
1055 0 : primaries = rp->GetPrimaries(nprimaries);
1056 0 : printf("\n%6d %8.2f %3d %2d %4.1f %4.1f %4.1f %4f %4f %2d : ",
1057 0 : rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
1058 0 : locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
1059 :
1060 0 : for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
1061 0 : printf("%d ", primaries[iprimary] ) ;
1062 : }
1063 0 : printf("\n") ;
1064 0 : }
1065 :
1066 : //Now plot CPV recPoints
1067 0 : printf("\n CPV clusters \n") ;
1068 0 : printf("Index Ene(MeV) Module X Y Z \n") ;
1069 0 : for (index = 0 ; index < fCPVRecPoints->GetEntries() ; index++) {
1070 0 : AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )fCPVRecPoints->At(index) ;
1071 :
1072 0 : TVector3 locpos;
1073 0 : rp->GetLocalPosition(locpos);
1074 :
1075 0 : printf("\n%6d %8.2f %2d %4.1f %4.1f %4.1f \n",
1076 0 : rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
1077 0 : locpos.X(), locpos.Y(), locpos.Z()) ;
1078 0 : }
1079 0 : }
1080 0 : }
1081 :
1082 :
1083 : //____________________________________________________________________________
1084 : void AliPHOSClusterizerv1::SetDistancesToBadChannels()
1085 : {
1086 : //For each EMC rec. point set the distance to the nearest bad crystal.
1087 : //Author: Boris Polichtchouk
1088 :
1089 16 : if(!fgCalibData->GetNumOfEmcBadChannels()) return;
1090 :
1091 0 : Int_t badIds[8000];
1092 0 : memset(badIds,0,8000*sizeof(Int_t));
1093 0 : fgCalibData->EmcBadChannelIds(badIds);
1094 :
1095 : AliPHOSEmcRecPoint* rp;
1096 :
1097 0 : TMatrixF gmat;
1098 0 : TVector3 gposRecPoint; // global (in ALICE frame) position of rec. point
1099 0 : TVector3 gposBadChannel; // global position of bad crystal
1100 0 : TVector3 dR;
1101 :
1102 : Float_t dist=0.,minDist=0.;
1103 0 : Int_t relid[4]={0,0,0,0} ;
1104 0 : TVector3 lpos ;
1105 0 : for(Int_t iRP=0; iRP<fEMCRecPoints->GetEntries(); iRP++){
1106 0 : rp = (AliPHOSEmcRecPoint*)fEMCRecPoints->At(iRP);
1107 : //evaluate distance to border
1108 0 : relid[0]=rp->GetPHOSMod() ;
1109 : Float_t xcorner=0,zcorner=0.;
1110 0 : Float_t xmin,xmax,zmin,zmax ;
1111 0 : fGeom->GetCrystalsEdges(rp->GetPHOSMod(),xmin, zmin, xmax, zmax) ;
1112 0 : rp->GetLocalPosition(lpos) ;
1113 0 : Float_t minDistX = 2.2+TMath::Min(lpos.X()-xmin,xmax-lpos.X()); //2.2 - crystal size
1114 0 : Float_t minDistZ = 2.2+TMath::Min(lpos.Z()-zmin,zmax-lpos.Z()); //2.2 - crystal size
1115 0 : minDist=TMath::Min(minDistX,minDistZ) ;
1116 0 : Int_t ixmin=32+Int_t(xmin/2.2) ;
1117 0 : for(Int_t iBad=0; iBad<fgCalibData->GetNumOfEmcBadChannels(); iBad++) {
1118 0 : fGeom->AbsToRelNumbering(badIds[iBad],relid) ;
1119 0 : if(relid[0]!=rp->GetPHOSMod()) //We can not evaluate global position directly since
1120 : continue ; //bad channels can be in the module which does not exist in simulations.
1121 0 : if(relid[2]<=ixmin) //non-existing part of mod 4
1122 : continue ;
1123 0 : rp->GetGlobalPosition(gposRecPoint,gmat);
1124 0 : fGeom->RelPosInAlice(badIds[iBad],gposBadChannel);
1125 0 : AliDebug(2,Form("BC position:[%.3f,%.3f,%.3f], RP position:[%.3f,%.3f,%.3f]. E=%.3f\n",
1126 : gposBadChannel.X(),gposBadChannel.Y(),gposBadChannel.Z(),
1127 : gposRecPoint.X(),gposRecPoint.Y(),gposRecPoint.Z(),rp->GetEnergy()));
1128 0 : dR = gposBadChannel-gposRecPoint;
1129 0 : dist = dR.Mag();
1130 0 : if(dist<minDist) minDist = dist;
1131 : }
1132 :
1133 0 : rp->SetDistanceToBadCrystal(minDist);
1134 0 : }
1135 :
1136 8 : }
1137 : //==================================================================================
1138 : Float_t AliPHOSClusterizerv1::Calibrate(Float_t amp, Int_t absId) const{
1139 : // Calibrate EMC digit, i.e. multiply its Amp by a factor read from CDB
1140 :
1141 616 : const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
1142 :
1143 : //Determine rel.position of the cell absolute ID
1144 308 : Int_t relId[4];
1145 308 : geom->AbsToRelNumbering(absId,relId);
1146 308 : Int_t module=relId[0];
1147 308 : Int_t row =relId[2];
1148 308 : Int_t column=relId[3];
1149 308 : if(relId[1]){ //CPV
1150 0 : Float_t calibration = fgCalibData->GetADCchannelCpv(module,column,row);
1151 0 : return amp*calibration ;
1152 : }
1153 : else{ //EMC
1154 308 : Float_t calibration = fgCalibData->GetADCchannelEmc(module,column,row);
1155 308 : return amp*calibration ;
1156 : }
1157 308 : }
1158 : //==================================================================================
1159 : Float_t AliPHOSClusterizerv1::CalibrateT(Float_t time, Int_t absId,Bool_t isLG)const{
1160 : // Calibrate time in EMC digit
1161 :
1162 176 : const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
1163 :
1164 : //Determine rel.position of the cell absolute ID
1165 176 : Int_t relId[4];
1166 176 : geom->AbsToRelNumbering(absId,relId);
1167 176 : Int_t module=relId[0];
1168 176 : Int_t row =relId[2];
1169 176 : Int_t column=relId[3];
1170 176 : if(relId[1]){ //CPV
1171 0 : return 0. ;
1172 : }
1173 : else{ //EMC
1174 352 : if(isLG)
1175 194 : time += fgCalibData->GetLGTimeShiftEmc(module,column,row);
1176 : else
1177 158 : time += fgCalibData->GetTimeShiftEmc(module,column,row);
1178 176 : return time ;
1179 : }
1180 176 : }
1181 :
|
