]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PHOS/AliPHOSClusterizerv1.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
events quality macro
[u/mrichter/AliRoot.git] / PHOS / AliPHOSClusterizerv1.cxx
... / ...
CommitLineData
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$
21 * Revision 1.95 2006/08/11 12:36:26 cvetan
22 * Update of the PHOS code needed in order to read and reconstruct the beam test raw data (i.e. without an existing galice.root)
23 *
24 * Revision 1.94 2006/08/07 12:27:49 hristov
25 * Removing obsolete code which affected the event numbering scheme
26 *
27 * Revision 1.93 2006/08/01 12:20:17 cvetan
28 * 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
29 *
30 * Revision 1.92 2006/04/29 20:26:46 hristov
31 * Separate EMC and CPV calibration (Yu.Kharlov)
32 *
33 * Revision 1.91 2006/04/22 10:30:17 hristov
34 * Add fEnergy to AliPHOSDigit and operate with EMC amplitude in energy units (Yu.Kharlov)
35 *
36 * Revision 1.90 2006/04/11 15:22:59 hristov
37 * run number in query set to -1: forces AliCDBManager to use its run number (A.Colla)
38 *
39 * Revision 1.89 2006/03/13 14:05:42 kharlov
40 * Calibration objects for EMC and CPV
41 *
42 * Revision 1.88 2006/01/11 08:54:52 hristov
43 * Additional protection in case no calibration entry was found
44 *
45 * Revision 1.87 2005/11/22 08:46:43 kharlov
46 * Updated to new CDB (Boris Polichtchouk)
47 *
48 * Revision 1.86 2005/11/14 21:52:43 hristov
49 * Coding conventions
50 *
51 * Revision 1.85 2005/09/27 16:08:08 hristov
52 * New version of CDB storage framework (A.Colla)
53 *
54 * Revision 1.84 2005/09/21 10:02:47 kharlov
55 * Reading calibration from CDB (Boris Polichtchouk)
56 *
57 * Revision 1.82 2005/09/02 15:43:13 kharlov
58 * Add comments in GetCalibrationParameters and Calibrate
59 *
60 * Revision 1.81 2005/09/02 14:32:07 kharlov
61 * Calibration of raw data
62 *
63 * Revision 1.80 2005/08/24 15:31:36 kharlov
64 * Setting raw digits flag
65 *
66 * Revision 1.79 2005/07/25 15:53:53 kharlov
67 * Read raw data
68 *
69 * Revision 1.78 2005/05/28 14:19:04 schutz
70 * Compilation warnings fixed by T.P.
71 *
72 */
73
74//*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
75//////////////////////////////////////////////////////////////////////////////
76// Clusterization class. Performs clusterization (collects neighbouring active cells) and
77// unfolds the clusters having several local maxima.
78// Results are stored in TreeR#, branches PHOSEmcRP (EMC recPoints),
79// PHOSCpvRP (CPV RecPoints) and AliPHOSClusterizer (Clusterizer with all
80// parameters including input digits branch title, thresholds etc.)
81// This TTask is normally called from Reconstructioner, but can as well be used in
82// standalone mode.
83// Use Case:
84// root [0] AliPHOSClusterizerv1 * cl = new AliPHOSClusterizerv1("galice.root", "recpointsname", "digitsname")
85// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
86// // reads gAlice from header file "galice.root", uses digits stored in the branch names "digitsname" (default = "Default")
87// // and saves recpoints in branch named "recpointsname" (default = "digitsname")
88// root [1] cl->ExecuteTask()
89// //finds RecPoints in all events stored in galice.root
90// root [2] cl->SetDigitsBranch("digits2")
91// //sets another title for Digitis (input) branch
92// root [3] cl->SetRecPointsBranch("recp2")
93// //sets another title four output branches
94// root [4] cl->SetEmcLocalMaxCut(0.03)
95// //set clusterization parameters
96// root [5] cl->ExecuteTask("deb all time")
97// //once more finds RecPoints options are
98// // deb - print number of found rec points
99// // deb all - print number of found RecPoints and some their characteristics
100// // time - print benchmarking results
101
102// --- ROOT system ---
103
104#include "TMath.h"
105#include "TMinuit.h"
106#include "TTree.h"
107#include "TBenchmark.h"
108
109// --- Standard library ---
110
111// --- AliRoot header files ---
112#include "AliLog.h"
113#include "AliPHOSGetter.h"
114#include "AliPHOSGeometry.h"
115#include "AliPHOSClusterizerv1.h"
116#include "AliPHOSEmcRecPoint.h"
117#include "AliPHOSCpvRecPoint.h"
118#include "AliPHOSDigit.h"
119#include "AliPHOSDigitizer.h"
120#include "AliPHOSCalibrationDB.h"
121#include "AliCDBManager.h"
122#include "AliCDBStorage.h"
123#include "AliCDBEntry.h"
124
125ClassImp(AliPHOSClusterizerv1)
126
127//____________________________________________________________________________
128AliPHOSClusterizerv1::AliPHOSClusterizerv1() :
129 AliPHOSClusterizer(),
130 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
131 fWrite(0), fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
132 fCalibData(0), fADCchanelEmc(0), fADCpedestalEmc(0),
133 fADCchanelCpv(0), fADCpedestalCpv(0), fEmcClusteringThreshold(0),
134 fCpvClusteringThreshold(0), fEmcMinE(0), fCpvMinE(0),
135 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
136 fW0CPV(0), fRecPointsInRun(0), fEmcTimeGate(0),
137 fIsOldRCUFormat(0)
138{
139 // default ctor (to be used mainly by Streamer)
140
141 InitParameters() ;
142 fDefaultInit = kTRUE ;
143}
144
145//____________________________________________________________________________
146AliPHOSClusterizerv1::AliPHOSClusterizerv1(const TString alirunFileName, const TString eventFolderName) :
147 AliPHOSClusterizer(alirunFileName, eventFolderName),
148 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
149 fWrite(0), fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
150 fCalibData(0), fADCchanelEmc(0), fADCpedestalEmc(0),
151 fADCchanelCpv(0), fADCpedestalCpv(0), fEmcClusteringThreshold(0),
152 fCpvClusteringThreshold(0), fEmcMinE(0), fCpvMinE(0),
153 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
154 fW0CPV(0), fRecPointsInRun(0), fEmcTimeGate(0),
155 fIsOldRCUFormat(0)
156{
157 // ctor with the indication of the file where header Tree and digits Tree are stored
158
159 InitParameters() ;
160 Init() ;
161 fDefaultInit = kFALSE ;
162}
163
164//____________________________________________________________________________
165 AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
166{
167 // dtor
168
169}
170//____________________________________________________________________________
171const TString AliPHOSClusterizerv1::BranchName() const
172{
173 return GetName();
174}
175
176//____________________________________________________________________________
177Float_t AliPHOSClusterizerv1::CalibrateEMC(Float_t amp, Int_t absId)
178{
179 // Convert EMC measured amplitude into real energy.
180 // Calibration parameters are taken from calibration data base for raw data,
181 // or from digitizer parameters for simulated data.
182
183 if(fCalibData){
184 Int_t relId[4];
185 AliPHOSGetter *gime = AliPHOSGetter::Instance();
186 gime->PHOSGeometry()->AbsToRelNumbering(absId,relId) ;
187 Int_t module = relId[0];
188 Int_t column = relId[3];
189 Int_t row = relId[2];
190 if(absId <= fEmcCrystals) { // this is EMC
191 fADCchanelEmc = fCalibData->GetADCchannelEmc (module,column,row);
192 return amp*fADCchanelEmc ;
193 }
194 }
195 else{ //simulation
196 if(absId <= fEmcCrystals) // this is EMC
197 return fADCpedestalEmc + amp*fADCchanelEmc ;
198 }
199 return 0;
200}
201
202//____________________________________________________________________________
203Float_t AliPHOSClusterizerv1::CalibrateCPV(Int_t amp, Int_t absId)
204{
205 // Convert digitized CPV amplitude into charge.
206 // Calibration parameters are taken from calibration data base for raw data,
207 // or from digitizer parameters for simulated data.
208
209 if(fCalibData){
210 Int_t relId[4];
211 AliPHOSGetter *gime = AliPHOSGetter::Instance();
212 gime->PHOSGeometry()->AbsToRelNumbering(absId,relId) ;
213 Int_t module = relId[0];
214 Int_t column = relId[3];
215 Int_t row = relId[2];
216 if(absId > fEmcCrystals) { // this is CPV
217 fADCchanelCpv = fCalibData->GetADCchannelCpv (module,column,row);
218 fADCpedestalCpv = fCalibData->GetADCpedestalCpv(module,column,row);
219 return fADCpedestalCpv + amp*fADCchanelCpv ;
220 }
221 }
222 else{ //simulation
223 if(absId > fEmcCrystals) // this is CPV
224 return fADCpedestalCpv+ amp*fADCchanelCpv ;
225 }
226 return 0;
227}
228
229//____________________________________________________________________________
230void AliPHOSClusterizerv1::Exec(Option_t *option)
231{
232 // Steering method to perform clusterization for events
233 // in the range from fFirstEvent to fLastEvent.
234 // This range is optionally set by SetEventRange().
235 // if fLastEvent=-1 (by default), then process events until the end.
236
237 if(strstr(option,"tim"))
238 gBenchmark->Start("PHOSClusterizer");
239
240 if(strstr(option,"print")) {
241 Print() ;
242 return ;
243 }
244
245 GetCalibrationParameters() ;
246
247 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
248 if (fRawReader == 0)
249 gime->SetRawDigits(kFALSE);
250 else
251 gime->SetRawDigits(kTRUE);
252
253 if (fLastEvent == -1)
254 fLastEvent = gime->MaxEvent() - 1 ;
255 else
256 fLastEvent = TMath::Min(fFirstEvent, gime->MaxEvent()); // one event at the time
257 Int_t nEvents = fLastEvent - fFirstEvent + 1;
258
259 Int_t ievent ;
260 for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
261 if (fRawReader == 0)
262 gime->Event(ievent ,"D"); // Read digits from simulated data
263 else {
264 gime->Event(fRawReader,"W",fIsOldRCUFormat); // Read digits from raw data
265 }
266 fNumberOfEmcClusters = fNumberOfCpvClusters = 0 ;
267
268 MakeClusters() ;
269
270 if(fToUnfold)
271 MakeUnfolding() ;
272
273 WriteRecPoints();
274
275 if(strstr(option,"deb"))
276 PrintRecPoints(option) ;
277
278 //increment the total number of recpoints per run
279 fRecPointsInRun += gime->EmcRecPoints()->GetEntriesFast() ;
280 fRecPointsInRun += gime->CpvRecPoints()->GetEntriesFast() ;
281 }
282
283 if(fWrite) //do not unload in "on flight" mode
284 Unload();
285
286 if(strstr(option,"tim")){
287 gBenchmark->Stop("PHOSClusterizer");
288 AliInfo(Form("took %f seconds for Clusterizing %f seconds per event \n",
289 gBenchmark->GetCpuTime("PHOSClusterizer"),
290 gBenchmark->GetCpuTime("PHOSClusterizer")/nEvents )) ;
291 }
292}
293
294//____________________________________________________________________________
295Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** maxAt, Float_t * maxAtEnergy,
296 Int_t nPar, Float_t * fitparameters) const
297{
298 // Calls TMinuit to fit the energy distribution of a cluster with several maxima
299 // The initial values for fitting procedure are set equal to the positions of local maxima.
300 // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
301
302
303 AliPHOSGetter * gime = AliPHOSGetter::Instance();
304 TClonesArray * digits = gime->Digits();
305
306
307 gMinuit->mncler(); // Reset Minuit's list of paramters
308 gMinuit->SetPrintLevel(-1) ; // No Printout
309 gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;
310 // To set the address of the minimization function
311
312 TList * toMinuit = new TList();
313 toMinuit->AddAt(emcRP,0) ;
314 toMinuit->AddAt(digits,1) ;
315
316 gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare
317
318 // filling initial values for fit parameters
319 AliPHOSDigit * digit ;
320
321 Int_t ierflg = 0;
322 Int_t index = 0 ;
323 Int_t nDigits = (Int_t) nPar / 3 ;
324
325 Int_t iDigit ;
326
327 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
328
329 for(iDigit = 0; iDigit < nDigits; iDigit++){
330 digit = maxAt[iDigit];
331
332 Int_t relid[4] ;
333 Float_t x = 0.;
334 Float_t z = 0.;
335 geom->AbsToRelNumbering(digit->GetId(), relid) ;
336 geom->RelPosInModule(relid, x, z) ;
337
338 Float_t energy = maxAtEnergy[iDigit] ;
339
340 gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ;
341 index++ ;
342 if(ierflg != 0){
343 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ;
344 return kFALSE;
345 }
346 gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ;
347 index++ ;
348 if(ierflg != 0){
349 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ;
350 return kFALSE;
351 }
352 gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;
353 index++ ;
354 if(ierflg != 0){
355 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ;
356 return kFALSE;
357 }
358 }
359
360 Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
361 // depends on it.
362 Double_t p1 = 1.0 ;
363 Double_t p2 = 0.0 ;
364
365 gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
366 gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
367 gMinuit->SetMaxIterations(5);
368 gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
369
370 gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
371
372 if(ierflg == 4){ // Minimum not found
373 Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" );
374 return kFALSE ;
375 }
376 for(index = 0; index < nPar; index++){
377 Double_t err ;
378 Double_t val ;
379 gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index
380 fitparameters[index] = val ;
381 }
382
383 delete toMinuit ;
384 return kTRUE;
385
386}
387
388//____________________________________________________________________________
389void AliPHOSClusterizerv1::GetCalibrationParameters()
390{
391 // Set calibration parameters:
392 // if calibration database exists, they are read from database,
393 // otherwise, they are taken from digitizer.
394 //
395 // It is a user responsilibity to open CDB before reconstruction, for example:
396 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
397
398 AliPHOSGetter * gime = AliPHOSGetter::Instance();
399 // fCalibData = new AliPHOSCalibData(gAlice->GetRunNumber()); //original
400
401 fCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
402
403 if(!fCalibData)
404 {
405 if ( !gime->Digitizer() )
406 gime->LoadDigitizer();
407 AliPHOSDigitizer * dig = gime->Digitizer();
408 fADCchanelEmc = dig->GetEMCchannel() ;
409 fADCpedestalEmc = dig->GetEMCpedestal();
410
411 fADCchanelCpv = dig->GetCPVchannel() ;
412 fADCpedestalCpv = dig->GetCPVpedestal() ;
413 }
414}
415
416//____________________________________________________________________________
417void AliPHOSClusterizerv1::Init()
418{
419 // Make all memory allocations which can not be done in default constructor.
420 // Attach the Clusterizer task to the list of PHOS tasks
421
422 AliPHOSGetter* gime = AliPHOSGetter::Instance() ;
423 if(!gime)
424 gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data());
425
426 AliPHOSGeometry * geom = gime->PHOSGeometry();
427
428 fEmcCrystals = geom->GetNModules() * geom->GetNCristalsInModule() ;
429
430 if(!gMinuit)
431 gMinuit = new TMinuit(100);
432
433 if ( !gime->Clusterizer() ) {
434 gime->PostClusterizer(this);
435 }
436}
437
438//____________________________________________________________________________
439void AliPHOSClusterizerv1::InitParameters()
440{
441
442 fNumberOfCpvClusters = 0 ;
443 fNumberOfEmcClusters = 0 ;
444
445 fCpvClusteringThreshold = 0.0;
446 fEmcClusteringThreshold = 0.2;
447
448 fEmcLocMaxCut = 0.03 ;
449 fCpvLocMaxCut = 0.03 ;
450
451 fEmcMinE = 0.01 ;
452 fCpvMinE = 0.0 ;
453
454 fW0 = 4.5 ;
455 fW0CPV = 4.0 ;
456
457 fEmcTimeGate = 1.e-8 ;
458
459 fToUnfold = kTRUE ;
460
461 fRecPointsInRun = 0 ;
462
463 fWrite = kTRUE ;
464
465 fCalibData = 0 ;
466
467 SetEventRange(0,-1) ;
468
469 fIsOldRCUFormat = kFALSE;
470}
471
472//____________________________________________________________________________
473Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
474{
475 // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
476 // = 1 are neighbour
477 // = 2 are not neighbour but do not continue searching
478 // neighbours are defined as digits having at least a common vertex
479 // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
480 // which is compared to a digit (d2) not yet in a cluster
481
482 AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
483
484 Int_t rv = 0 ;
485
486 Int_t relid1[4] ;
487 geom->AbsToRelNumbering(d1->GetId(), relid1) ;
488
489 Int_t relid2[4] ;
490 geom->AbsToRelNumbering(d2->GetId(), relid2) ;
491
492 if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
493 Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
494 Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
495
496 if (( coldiff <= 1 ) && ( rowdiff <= 1 )){
497 if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fEmcTimeGate))
498 rv = 1 ;
499 }
500 else {
501 if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
502 rv = 2; // Difference in row numbers is too large to look further
503 }
504
505 }
506 else {
507
508 if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) )
509 rv=2 ;
510
511 }
512
513 return rv ;
514}
515//____________________________________________________________________________
516void AliPHOSClusterizerv1::CleanDigits(TClonesArray * digits)
517{
518 // Remove digits with amplitudes below threshold
519
520 for(Int_t i=0; i<digits->GetEntriesFast(); i++){
521 AliPHOSDigit * digit = static_cast<AliPHOSDigit*>(digits->At(i)) ;
522 if ( (IsInEmc(digit) && CalibrateEMC(digit->GetEnergy(),digit->GetId()) < fEmcMinE) ||
523 (IsInCpv(digit) && CalibrateCPV(digit->GetAmp() ,digit->GetId()) < fCpvMinE) )
524 digits->RemoveAt(i) ;
525 }
526 digits->Compress() ;
527 for (Int_t i = 0 ; i < digits->GetEntriesFast() ; i++) {
528 AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
529 digit->SetIndexInList(i) ;
530 }
531
532}
533//____________________________________________________________________________
534Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
535{
536 // Tells if (true) or not (false) the digit is in a PHOS-EMC module
537
538 Bool_t rv = kFALSE ;
539 AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
540
541 Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
542
543 if(digit->GetId() <= nEMC ) rv = kTRUE;
544
545 return rv ;
546}
547
548//____________________________________________________________________________
549Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
550{
551 // Tells if (true) or not (false) the digit is in a PHOS-CPV module
552
553 Bool_t rv = kFALSE ;
554
555 AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
556
557 Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
558
559 if(digit->GetId() > nEMC ) rv = kTRUE;
560
561 return rv ;
562}
563
564//____________________________________________________________________________
565void AliPHOSClusterizerv1::Unload()
566{
567 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
568 gime->PhosLoader()->UnloadDigits() ;
569 gime->PhosLoader()->UnloadRecPoints() ;
570}
571
572//____________________________________________________________________________
573void AliPHOSClusterizerv1::WriteRecPoints()
574{
575
576 // Creates new branches with given title
577 // fills and writes into TreeR.
578
579 AliPHOSGetter * gime = AliPHOSGetter::Instance();
580
581 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
582 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
583 TClonesArray * digits = gime->Digits() ;
584
585
586 Int_t index ;
587 //Evaluate position, dispersion and other RecPoint properties..
588 Int_t nEmc = emcRecPoints->GetEntriesFast();
589 for(index = 0; index < nEmc; index++){
590 AliPHOSEmcRecPoint * rp = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) );
591 rp->Purify(fEmcMinE) ;
592 if(rp->GetMultiplicity()>0.) //If this RP is not empty
593 rp->EvalAll(fW0,digits) ;
594 else{
595 emcRecPoints->RemoveAt(index) ;
596 delete rp ;
597 }
598 }
599 emcRecPoints->Compress() ;
600 emcRecPoints->Sort() ;
601 // emcRecPoints->Expand(emcRecPoints->GetEntriesFast()) ;
602 for(index = 0; index < emcRecPoints->GetEntries(); index++){
603 dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) )->SetIndexInList(index) ;
604 }
605
606 //Now the same for CPV
607 for(index = 0; index < cpvRecPoints->GetEntries(); index++){
608 AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) );
609 rp->EvalAll(fW0CPV,digits) ;
610 }
611 cpvRecPoints->Sort() ;
612
613 for(index = 0; index < cpvRecPoints->GetEntries(); index++)
614 dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) )->SetIndexInList(index) ;
615
616 cpvRecPoints->Expand(cpvRecPoints->GetEntriesFast()) ;
617
618 if(fWrite){ //We write TreeR
619 TTree * treeR = gime->TreeR();
620
621 Int_t bufferSize = 32000 ;
622 Int_t splitlevel = 0 ;
623
624 //First EMC
625 TBranch * emcBranch = treeR->Branch("PHOSEmcRP","TObjArray",&emcRecPoints,bufferSize,splitlevel);
626 emcBranch->SetTitle(BranchName());
627
628 //Now CPV branch
629 TBranch * cpvBranch = treeR->Branch("PHOSCpvRP","TObjArray",&cpvRecPoints,bufferSize,splitlevel);
630 cpvBranch->SetTitle(BranchName());
631
632 emcBranch ->Fill() ;
633 cpvBranch ->Fill() ;
634
635 gime->WriteRecPoints("OVERWRITE");
636 gime->WriteClusterizer("OVERWRITE");
637 }
638}
639
640//____________________________________________________________________________
641void AliPHOSClusterizerv1::MakeClusters()
642{
643 // Steering method to construct the clusters stored in a list of Reconstructed Points
644 // A cluster is defined as a list of neighbour digits
645
646
647 AliPHOSGetter * gime = AliPHOSGetter::Instance();
648
649 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
650 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
651
652 emcRecPoints->Delete() ;
653 cpvRecPoints->Delete() ;
654
655 TClonesArray * digits = gime->Digits() ;
656
657 //Remove digits below threshold
658 CleanDigits(digits) ;
659
660 TClonesArray * digitsC = static_cast<TClonesArray*>( digits->Clone() ) ;
661
662 // Clusterization starts
663
664 TIter nextdigit(digitsC) ;
665 AliPHOSDigit * digit ;
666 Bool_t notremoved = kTRUE ;
667
668 while ( (digit = dynamic_cast<AliPHOSDigit *>( nextdigit()) ) ) { // scan over the list of digitsC
669
670
671 AliPHOSRecPoint * clu = 0 ;
672
673 TArrayI clusterdigitslist(1500) ;
674 Int_t index ;
675
676 if (( IsInEmc (digit) &&
677 CalibrateEMC(digit->GetEnergy(),digit->GetId()) > fEmcClusteringThreshold ) ||
678 ( IsInCpv (digit) &&
679 CalibrateCPV(digit->GetAmp() ,digit->GetId()) > fCpvClusteringThreshold ) ) {
680 Int_t iDigitInCluster = 0 ;
681
682 if ( IsInEmc(digit) ) {
683 // start a new EMC RecPoint
684 if(fNumberOfEmcClusters >= emcRecPoints->GetSize())
685 emcRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
686
687 emcRecPoints->AddAt(new AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
688 clu = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(fNumberOfEmcClusters) ) ;
689 fNumberOfEmcClusters++ ;
690 clu->AddDigit(*digit, CalibrateEMC(digit->GetEnergy(),digit->GetId())) ;
691 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
692 iDigitInCluster++ ;
693 digitsC->Remove(digit) ;
694
695 } else {
696
697 // start a new CPV cluster
698 if(fNumberOfCpvClusters >= cpvRecPoints->GetSize())
699 cpvRecPoints->Expand(2*fNumberOfCpvClusters+1);
700
701 cpvRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
702
703 clu = dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ;
704 fNumberOfCpvClusters++ ;
705 clu->AddDigit(*digit, CalibrateCPV(digit->GetAmp(),digit->GetId()) ) ;
706 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
707 iDigitInCluster++ ;
708 digitsC->Remove(digit) ;
709 nextdigit.Reset() ;
710
711 // Here we remove remaining EMC digits, which cannot make a cluster
712
713 if( notremoved ) {
714 while( ( digit = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) {
715 if( IsInEmc(digit) )
716 digitsC->Remove(digit) ;
717 else
718 break ;
719 }
720 notremoved = kFALSE ;
721 }
722
723 } // else
724
725 nextdigit.Reset() ;
726
727 AliPHOSDigit * digitN ;
728 index = 0 ;
729 while (index < iDigitInCluster){ // scan over digits already in cluster
730 digit = dynamic_cast<AliPHOSDigit*>( digits->At(clusterdigitslist[index]) ) ;
731 index++ ;
732 while ( (digitN = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) { // scan over the reduced list of digits
733 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
734 switch (ineb ) {
735 case 0 : // not a neighbour
736 break ;
737 case 1 : // are neighbours
738 if (IsInEmc (digitN))
739 clu->AddDigit(*digitN, CalibrateEMC( digitN->GetEnergy(), digitN->GetId() ) );
740 else
741 clu->AddDigit(*digitN, CalibrateCPV( digitN->GetAmp() , digitN->GetId() ) );
742
743 clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ;
744 iDigitInCluster++ ;
745 digitsC->Remove(digitN) ;
746 break ;
747 case 2 : // too far from each other
748 goto endofloop;
749 } // switch
750
751 } // while digitN
752
753 endofloop: ;
754 nextdigit.Reset() ;
755
756 } // loop over cluster
757
758 } // energy theshold
759
760
761 } // while digit
762
763 delete digitsC ;
764
765}
766
767//____________________________________________________________________________
768void AliPHOSClusterizerv1::MakeUnfolding()
769{
770 // Unfolds clusters using the shape of an ElectroMagnetic shower
771 // Performs unfolding of all EMC/CPV clusters
772
773 AliPHOSGetter * gime = AliPHOSGetter::Instance();
774
775 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
776
777 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
778 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
779 TClonesArray * digits = gime->Digits() ;
780
781 // Unfold first EMC clusters
782 if(fNumberOfEmcClusters > 0){
783
784 Int_t nModulesToUnfold = geom->GetNModules() ;
785
786 Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
787 Int_t index ;
788 for(index = 0 ; index < numberofNotUnfolded ; index++){
789
790 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) ) ;
791 if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
792 break ;
793
794 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
795 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
796 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
797 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,digits) ;
798
799 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
800 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
801 emcRecPoints->Remove(emcRecPoint);
802 emcRecPoints->Compress() ;
803 index-- ;
804 fNumberOfEmcClusters -- ;
805 numberofNotUnfolded-- ;
806 }
807 else{
808 emcRecPoint->SetNExMax(1) ; //Only one local maximum
809 }
810
811 delete[] maxAt ;
812 delete[] maxAtEnergy ;
813 }
814 }
815 // Unfolding of EMC clusters finished
816
817
818 // Unfold now CPV clusters
819 if(fNumberOfCpvClusters > 0){
820
821 Int_t nModulesToUnfold = geom->GetNModules() ;
822
823 Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
824 Int_t index ;
825 for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
826
827 AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( cpvRecPoints->At(index) ) ;
828
829 if(recPoint->GetPHOSMod()> nModulesToUnfold)
830 break ;
831
832 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ;
833
834 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
835 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
836 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
837 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,digits) ;
838
839 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
840 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
841 cpvRecPoints->Remove(emcRecPoint);
842 cpvRecPoints->Compress() ;
843 index-- ;
844 numberofCpvNotUnfolded-- ;
845 fNumberOfCpvClusters-- ;
846 }
847
848 delete[] maxAt ;
849 delete[] maxAtEnergy ;
850 }
851 }
852 //Unfolding of Cpv clusters finished
853
854}
855
856//____________________________________________________________________________
857Double_t AliPHOSClusterizerv1::ShowerShape(Double_t r)
858{
859 // Shape of the shower (see PHOS TDR)
860 // If you change this function, change also the gradient evaluation in ChiSquare()
861
862 Double_t r4 = r*r*r*r ;
863 Double_t r295 = TMath::Power(r, 2.95) ;
864 Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
865 return shape ;
866}
867
868//____________________________________________________________________________
869void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
870 Int_t nMax,
871 AliPHOSDigit ** maxAt,
872 Float_t * maxAtEnergy)
873{
874 // Performs the unfolding of a cluster with nMax overlapping showers
875
876 AliPHOSGetter * gime = AliPHOSGetter::Instance();
877
878 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
879
880 const TClonesArray * digits = gime->Digits() ;
881 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
882 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
883
884 Int_t nPar = 3 * nMax ;
885 Float_t * fitparameters = new Float_t[nPar] ;
886
887 Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
888 if( !rv ) {
889 // Fit failed, return and remove cluster
890 iniEmc->SetNExMax(-1) ;
891 delete[] fitparameters ;
892 return ;
893 }
894
895 // create ufolded rec points and fill them with new energy lists
896 // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
897 // and later correct this number in acordance with actual energy deposition
898
899 Int_t nDigits = iniEmc->GetMultiplicity() ;
900 Float_t * efit = new Float_t[nDigits] ;
901 Float_t xDigit=0.,zDigit=0.,distance=0. ;
902 Float_t xpar=0.,zpar=0.,epar=0. ;
903 Int_t relid[4] ;
904 AliPHOSDigit * digit = 0 ;
905 Int_t * emcDigits = iniEmc->GetDigitsList() ;
906
907 Int_t iparam ;
908 Int_t iDigit ;
909 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
910 digit = dynamic_cast<AliPHOSDigit*>( digits->At(emcDigits[iDigit] ) ) ;
911 geom->AbsToRelNumbering(digit->GetId(), relid) ;
912 geom->RelPosInModule(relid, xDigit, zDigit) ;
913 efit[iDigit] = 0;
914
915 iparam = 0 ;
916 while(iparam < nPar ){
917 xpar = fitparameters[iparam] ;
918 zpar = fitparameters[iparam+1] ;
919 epar = fitparameters[iparam+2] ;
920 iparam += 3 ;
921 distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ;
922 distance = TMath::Sqrt(distance) ;
923 efit[iDigit] += epar * ShowerShape(distance) ;
924 }
925 }
926
927
928 // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
929 // so that energy deposited in each cell is distributed betwin new clusters proportionally
930 // to its contribution to efit
931
932 Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
933 Float_t ratio ;
934
935 iparam = 0 ;
936 while(iparam < nPar ){
937 xpar = fitparameters[iparam] ;
938 zpar = fitparameters[iparam+1] ;
939 epar = fitparameters[iparam+2] ;
940 iparam += 3 ;
941
942 AliPHOSEmcRecPoint * emcRP = 0 ;
943
944 if(iniEmc->IsEmc()){ //create new entries in fEmcRecPoints...
945
946 if(fNumberOfEmcClusters >= emcRecPoints->GetSize())
947 emcRecPoints->Expand(2*fNumberOfEmcClusters) ;
948
949 (*emcRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
950 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(fNumberOfEmcClusters) ) ;
951 fNumberOfEmcClusters++ ;
952 emcRP->SetNExMax((Int_t)nPar/3) ;
953 }
954 else{//create new entries in fCpvRecPoints
955 if(fNumberOfCpvClusters >= cpvRecPoints->GetSize())
956 cpvRecPoints->Expand(2*fNumberOfCpvClusters) ;
957
958 (*cpvRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
959 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ;
960 fNumberOfCpvClusters++ ;
961 }
962
963 Float_t eDigit ;
964 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
965 digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) ) ;
966 geom->AbsToRelNumbering(digit->GetId(), relid) ;
967 geom->RelPosInModule(relid, xDigit, zDigit) ;
968 distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ;
969 distance = TMath::Sqrt(distance) ;
970 ratio = epar * ShowerShape(distance) / efit[iDigit] ;
971 eDigit = emcEnergies[iDigit] * ratio ;
972 emcRP->AddDigit( *digit, eDigit ) ;
973 }
974 }
975
976 delete[] fitparameters ;
977 delete[] efit ;
978
979}
980
981//_____________________________________________________________________________
982void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
983{
984 // Calculates the Chi square for the cluster unfolding minimization
985 // Number of parameters, Gradient, Chi squared, parameters, what to do
986
987 TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
988
989 AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) ) ;
990 TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) ) ;
991
992
993
994 // AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
995
996 Int_t * emcDigits = emcRP->GetDigitsList() ;
997
998 Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
999
1000 Float_t * emcEnergies = emcRP->GetEnergiesList() ;
1001
1002 const AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
1003 fret = 0. ;
1004 Int_t iparam ;
1005
1006 if(iflag == 2)
1007 for(iparam = 0 ; iparam < nPar ; iparam++)
1008 Grad[iparam] = 0 ; // Will evaluate gradient
1009
1010 Double_t efit ;
1011
1012 AliPHOSDigit * digit ;
1013 Int_t iDigit ;
1014
1015 for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
1016
1017 digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) );
1018
1019 Int_t relid[4] ;
1020 Float_t xDigit ;
1021 Float_t zDigit ;
1022
1023 geom->AbsToRelNumbering(digit->GetId(), relid) ;
1024
1025 geom->RelPosInModule(relid, xDigit, zDigit) ;
1026
1027 if(iflag == 2){ // calculate gradient
1028 Int_t iParam = 0 ;
1029 efit = 0 ;
1030 while(iParam < nPar ){
1031 Double_t distance = (xDigit - x[iParam]) * (xDigit - x[iParam]) ;
1032 iParam++ ;
1033 distance += (zDigit - x[iParam]) * (zDigit - x[iParam]) ;
1034 distance = TMath::Sqrt( distance ) ;
1035 iParam++ ;
1036 efit += x[iParam] * ShowerShape(distance) ;
1037 iParam++ ;
1038 }
1039 Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
1040 iParam = 0 ;
1041 while(iParam < nPar ){
1042 Double_t xpar = x[iParam] ;
1043 Double_t zpar = x[iParam+1] ;
1044 Double_t epar = x[iParam+2] ;
1045 Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
1046 Double_t shape = sum * ShowerShape(dr) ;
1047 Double_t r4 = dr*dr*dr*dr ;
1048 Double_t r295 = TMath::Power(dr,2.95) ;
1049 Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
1050 0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
1051
1052 Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x
1053 iParam++ ;
1054 Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z
1055 iParam++ ;
1056 Grad[iParam] += shape ; // Derivative over energy
1057 iParam++ ;
1058 }
1059 }
1060 efit = 0;
1061 iparam = 0 ;
1062
1063 while(iparam < nPar ){
1064 Double_t xpar = x[iparam] ;
1065 Double_t zpar = x[iparam+1] ;
1066 Double_t epar = x[iparam+2] ;
1067 iparam += 3 ;
1068 Double_t distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ;
1069 distance = TMath::Sqrt(distance) ;
1070 efit += epar * ShowerShape(distance) ;
1071 }
1072
1073 fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
1074 // Here we assume, that sigma = sqrt(E)
1075 }
1076
1077}
1078
1079//____________________________________________________________________________
1080void AliPHOSClusterizerv1::Print(const Option_t *)const
1081{
1082 // Print clusterizer parameters
1083
1084 TString message ;
1085 TString taskName(GetName()) ;
1086 taskName.ReplaceAll(Version(), "") ;
1087
1088 if( strcmp(GetName(), "") !=0 ) {
1089 // Print parameters
1090 message = "\n--------------- %s %s -----------\n" ;
1091 message += "Clusterizing digits from the file: %s\n" ;
1092 message += " Branch: %s\n" ;
1093 message += " EMC Clustering threshold = %f\n" ;
1094 message += " EMC Local Maximum cut = %f\n" ;
1095 message += " EMC Logarothmic weight = %f\n" ;
1096 message += " CPV Clustering threshold = %f\n" ;
1097 message += " CPV Local Maximum cut = %f\n" ;
1098 message += " CPV Logarothmic weight = %f\n" ;
1099 if(fToUnfold)
1100 message += " Unfolding on\n" ;
1101 else
1102 message += " Unfolding off\n" ;
1103
1104 message += "------------------------------------------------------------------" ;
1105 }
1106 else
1107 message = " AliPHOSClusterizerv1 not initialized " ;
1108
1109 AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(),
1110 taskName.Data(),
1111 GetTitle(),
1112 taskName.Data(),
1113 GetName(),
1114 fEmcClusteringThreshold,
1115 fEmcLocMaxCut,
1116 fW0,
1117 fCpvClusteringThreshold,
1118 fCpvLocMaxCut,
1119 fW0CPV )) ;
1120}
1121
1122
1123//____________________________________________________________________________
1124void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
1125{
1126 // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
1127
1128 AliPHOSGetter * gime = AliPHOSGetter::Instance();
1129
1130 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
1131 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
1132
1133 AliInfo(Form("\nevent %d \n Found %d EMC RecPoints and %d CPV RecPoints",
1134 gAlice->GetEvNumber(),
1135 emcRecPoints->GetEntriesFast(),
1136 cpvRecPoints->GetEntriesFast() )) ;
1137
1138 fRecPointsInRun += emcRecPoints->GetEntriesFast() ;
1139 fRecPointsInRun += cpvRecPoints->GetEntriesFast() ;
1140
1141
1142 if(strstr(option,"all")) {
1143 printf("\n EMC clusters \n") ;
1144 printf("Index Ene(MeV) Multi Module X Y Z Lambdas_1 Lambda_2 # of prim Primaries list\n") ;
1145 Int_t index ;
1146 for (index = 0 ; index < emcRecPoints->GetEntries() ; index++) {
1147 AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )emcRecPoints->At(index) ;
1148 TVector3 locpos;
1149 rp->GetLocalPosition(locpos);
1150 Float_t lambda[2];
1151 rp->GetElipsAxis(lambda);
1152 Int_t * primaries;
1153 Int_t nprimaries;
1154 primaries = rp->GetPrimaries(nprimaries);
1155 printf("\n%6d %8.2f %3d %2d %4.1f %4.1f %4.1f %4f %4f %2d : ",
1156 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
1157 locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
1158
1159 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
1160 printf("%d ", primaries[iprimary] ) ;
1161 }
1162 printf("\n") ;
1163 }
1164
1165 //Now plot CPV recPoints
1166 printf("\n CPV clusters \n") ;
1167 printf("Index Ene(MeV) Module X Y Z \n") ;
1168 for (index = 0 ; index < cpvRecPoints->GetEntries() ; index++) {
1169 AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )cpvRecPoints->At(index) ;
1170
1171 TVector3 locpos;
1172 rp->GetLocalPosition(locpos);
1173
1174 printf("\n%6d %8.2f %2d %4.1f %4.1f %4.1f \n",
1175 rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
1176 locpos.X(), locpos.Y(), locpos.Z()) ;
1177 }
1178 }
1179}
1180