git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Updated version with the Bari/Salerno model as default
[u/mrichter/AliRoot.git] / PHOS / AliPHOSClusterizerv1.cxx
CommitLineData
d15a28e7 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
b2a60966 16/* $Id$ */
17
d72dfbc3 18/* $Log:
19 1 October 2000. Yuri Kharlov:
20 AreNeighbours()
21 PPSD upper layer is considered if number of layers>1
22
23 18 October 2000. Yuri Kharlov:
24 AliPHOSClusterizerv1()
25 CPV clusterizing parameters added
26
27 MakeClusters()
28 After first PPSD digit remove EMC digits only once
29*/
9a1398dd 30//*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
d15a28e7 31//////////////////////////////////////////////////////////////////////////////
9a1398dd 32// Clusterization class. Performs clusterization (collects neighbouring active cells) and
33// unfolding of the clusters with several local maxima.
34// results are stored in TreeR#, branches PHOSEmcRP (EMC recPoints),
35// PHOSCpvRP (CPV RecPoints) and AliPHOSClusterizer (Clusterizer with all
f035f6ce 36// parameters including input digits branch title, thresholds etc.)
9a1398dd 37// This TTask normally called from Reconstructioner, but as well can be used it in
38// standalone mode:
39// root [0] AliPHOSClusterizerv1 * cl = new AliPHOSClusterizerv1("galice.root")
40// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
41// //reads gAlice from header file "..."
42// root [1] cl->ExecuteTask()
43// //finds RecPoints in all events stored in galice.root
44// root [2] cl->SetDigitsBranch("digits2")
f035f6ce 45// //sets another title for Digitis (input) branch
9a1398dd 46// root [3] cl->SetRecPointsBranch("recp2")
f035f6ce 47// //sets another title four output branches
9a1398dd 48// root [4] cl->SetEmcLocalMaxCut(0.03)
49// //set clusterization parameters
50// root [5] cl->ExecuteTask("deb all time")
51// //once more finds RecPoints options are
52// // deb - print number of found rec points
53// // deb all - print number of found RecPoints and some their characteristics
54// // time - print benchmarking results
ed4205d8 55
d15a28e7 56// --- ROOT system ---
57
9a1398dd 58#include "TROOT.h"
59#include "TFile.h"
60#include "TFolder.h"
d15a28e7 61#include "TMath.h"
9a1398dd 62#include "TMinuit.h"
63#include "TTree.h"
64#include "TSystem.h"
65#include "TBenchmark.h"
d15a28e7 66
67// --- Standard library ---
68
de9ec31b 69#include <iostream.h>
9a1398dd 70#include <iomanip.h>
d15a28e7 71
72// --- AliRoot header files ---
73
74#include "AliPHOSClusterizerv1.h"
9a1398dd 75#include "AliPHOSCpvRecPoint.h"
d15a28e7 76#include "AliPHOSDigit.h"
9a1398dd 77#include "AliPHOSDigitizer.h"
d15a28e7 78#include "AliPHOSEmcRecPoint.h"
9a1398dd 79#include "AliPHOS.h"
d15a28e7 80#include "AliPHOSPpsdRecPoint.h"
9a1398dd 81#include "AliRun.h"
d15a28e7 82
83ClassImp(AliPHOSClusterizerv1)
f035f6ce 84
d15a28e7 85//____________________________________________________________________________
9a1398dd 86 AliPHOSClusterizerv1::AliPHOSClusterizerv1():AliPHOSClusterizer()
d15a28e7 87{
f035f6ce 88 // default ctor (to be used mainly by Streamer)
89 SetName("AliPHOSClusterizer");
9a1398dd 90 SetTitle("Version 1") ;
f035f6ce 91
9a1398dd 92 fNumberOfCpvClusters = 0 ;
d15a28e7 93 fNumberOfEmcClusters = 0 ;
f035f6ce 94
9a1398dd 95 fCpvClusteringThreshold = 0.0;
c3c187e5 96 fEmcClusteringThreshold = 0.2;
c3c187e5 97 fPpsdClusteringThreshold = 0.0000002 ;
9a1398dd 98
99 fEmcLocMaxCut = 0.03 ;
100 fCpvLocMaxCut = 0.03 ;
101
102 fW0 = 4.5 ;
103 fW0CPV = 4.0 ;
f035f6ce 104
9a1398dd 105 fGeom = 0 ;
106
107 fDigits = 0 ;
108 fDigitizer = 0 ;
109 fEmcRecPoints = 0 ;
110 fCpvRecPoints = 0 ;
111
112 fIsInitialized = kFALSE ;
113
114}
115//____________________________________________________________________________
f035f6ce 116AliPHOSClusterizerv1::AliPHOSClusterizerv1(const char* headerFile,const char* digitsFile):AliPHOSClusterizer()
9a1398dd 117{
118 SetName("AliPHOSClusterizer");
119 SetTitle("Version 1") ;
120
121 fNumberOfCpvClusters = 0 ;
122 fNumberOfEmcClusters = 0 ;
f035f6ce 123
ed4205d8 124 fCpvClusteringThreshold = 0.0;
9a1398dd 125 fEmcClusteringThreshold = 0.2;
126 fPpsdClusteringThreshold = 0.0000002 ;
127
128 fEmcLocMaxCut = 0.03 ;
129 fCpvLocMaxCut = 0.03 ;
130
d15a28e7 131 fW0 = 4.5 ;
c3c187e5 132 fW0CPV = 4.0 ;
9a1398dd 133
134 fToUnfold = kTRUE ;
135
136 fHeaderFileName = headerFile ;
137 fDigitsBranchTitle = digitsFile ;
138
139 TFile * file = (TFile*) gROOT->GetFile(fHeaderFileName.Data() ) ;
140
141 if(file == 0){
f035f6ce 142 if(fHeaderFileName.Contains("rfio")) // if we read file using HPSS
143 file = TFile::Open(fHeaderFileName.Data(),"update") ;
144 else
145 file = new TFile(fHeaderFileName.Data(),"update") ;
9a1398dd 146 gAlice = (AliRun *) file->Get("gAlice") ;
147 }
148
149 AliPHOS * phos = (AliPHOS *) gAlice->GetDetector("PHOS") ;
150 fGeom = AliPHOSGeometry::GetInstance(phos->GetGeometry()->GetName(),phos->GetGeometry()->GetTitle() );
151
152 fDigits = new TClonesArray("AliPHOSDigit",10) ;
153 fDigitizer = new AliPHOSDigitizer() ;
154 fEmcRecPoints = new TObjArray(200) ;
155 fCpvRecPoints = new TObjArray(200) ;
7f5d510c 156 fEmcRecPoints->SetOwner(); // This lets Clear() really detete rec.points in array
157 fCpvRecPoints->SetOwner();
9a1398dd 158
159 if(!gMinuit) gMinuit = new TMinuit(100) ;
160
161 // add Task to //root/Tasks folder
162 TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ;
163 roottasks->Add(this) ;
164
165
166 fIsInitialized = kTRUE ;
167
168}
169//____________________________________________________________________________
170void AliPHOSClusterizerv1::Exec(Option_t * option){
171 // Steerign function
172
173 if(!fIsInitialized) Init() ;
174
175 if(strstr(option,"tim"))
176 gBenchmark->Start("PHOSClusterizer");
177
178 Int_t nEvents = (Int_t) gAlice->TreeE()->GetEntries() ;
179
180 for(fEvent = 0;fEvent< nEvents; fEvent++){
181 if(!ReadDigits()) //reads digits for event fEvent
182 return;
183 MakeClusters() ;
184
185 if(fToUnfold) MakeUnfolding() ;
186 WriteRecPoints() ;
187 if(strstr(option,"deb"))
188 PrintRecPoints(option) ;
189 }
190
191 if(strstr(option,"tim")){
192 gBenchmark->Stop("PHOSClusterizer");
193 cout << "AliPHOSClusterizer:" << endl ;
194 cout << " took " << gBenchmark->GetCpuTime("PHOSClusterizer") << " seconds for Clusterizing "
195 << gBenchmark->GetCpuTime("PHOSClusterizer")/nEvents << " seconds per event " << endl ;
196 cout << endl ;
197 }
198
199
200}
201
202//____________________________________________________________________________
203Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, int * maxAt, Float_t * maxAtEnergy,
204 Int_t nPar, Float_t * fitparameters)
205{
206 // Calls TMinuit to fit the energy distribution of a cluster with several maxima
f035f6ce 207 // the initial values for fitting procedure are set in the positions of local maxima.
208 // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
9a1398dd 209
210 gMinuit->mncler(); // Reset Minuit's list of paramters
211 gMinuit->SetPrintLevel(-1) ; // No Printout
212 gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;
213 // To set the address of the minimization function
214
215 TList * toMinuit = new TList();
216 toMinuit->AddAt(emcRP,0) ;
217 toMinuit->AddAt(fDigits,1) ;
218
219 gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare
220
221 // filling initial values for fit parameters
222 AliPHOSDigit * digit ;
223
224 Int_t ierflg = 0;
225 Int_t index = 0 ;
226 Int_t nDigits = (Int_t) nPar / 3 ;
227
228 Int_t iDigit ;
229
230
231 for(iDigit = 0; iDigit < nDigits; iDigit++){
232 digit = (AliPHOSDigit *) maxAt[iDigit];
233
234 Int_t relid[4] ;
235 Float_t x ;
236 Float_t z ;
237 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
238 fGeom->RelPosInModule(relid, x, z) ;
239
240 Float_t energy = maxAtEnergy[iDigit] ;
241
242 gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ;
243 index++ ;
244 if(ierflg != 0){
245 cout << "PHOS Unfolding> Unable to set initial value for fit procedure : x = " << x << endl ;
246 return kFALSE;
247 }
248 gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ;
249 index++ ;
250 if(ierflg != 0){
251 cout << "PHOS Unfolding> Unable to set initial value for fit procedure : z = " << z << endl ;
252 return kFALSE;
253 }
254 gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;
255 index++ ;
256 if(ierflg != 0){
257 cout << "PHOS Unfolding> Unable to set initial value for fit procedure : energy = " << energy << endl ;
258 return kFALSE;
259 }
260 }
261
262 Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
263 // depends on it.
264 Double_t p1 = 1.0 ;
265 Double_t p2 = 0.0 ;
266
267 gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
268 gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
269 gMinuit->SetMaxIterations(5);
270 gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
271
272 gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
273
274 if(ierflg == 4){ // Minimum not found
275 cout << "PHOS Unfolding> Fit not converged, cluster abandoned "<< endl ;
276 return kFALSE ;
277 }
278 for(index = 0; index < nPar; index++){
279 Double_t err ;
280 Double_t val ;
281 gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index
282 fitparameters[index] = val ;
283 }
284
285 delete toMinuit ;
286 return kTRUE;
c3c187e5 287
d15a28e7 288}
289
290//____________________________________________________________________________
9a1398dd 291void AliPHOSClusterizerv1::Init(){
f035f6ce 292 //Make all memory allocations which can not be done in default constructor.
9a1398dd 293 if(!fIsInitialized){
294 if(fHeaderFileName.IsNull())
295 fHeaderFileName = "galice.root" ;
296
297
298 TFile * file = (TFile*) gROOT->GetFile(fHeaderFileName.Data() ) ;
299
300 if(file == 0){
301 file = new TFile(fHeaderFileName.Data(),"update") ;
302 gAlice = (AliRun *) file->Get("gAlice") ;
303 }
304
305 AliPHOS * phos = (AliPHOS *) gAlice->GetDetector("PHOS") ;
306 fGeom = AliPHOSGeometry::GetInstance(phos->GetGeometry()->GetName(),phos->GetGeometry()->GetTitle() );
307
308 fDigits = new TClonesArray("AliPHOSDigit",10) ;
309 fDigitizer = new AliPHOSDigitizer() ;
310 fEmcRecPoints = new TObjArray(200) ;
311 fCpvRecPoints = new TObjArray(200) ;
7f5d510c 312 fEmcRecPoints->SetOwner(); // This lets Clear() really detete rec.points in array
313 fCpvRecPoints->SetOwner();
9a1398dd 314
315 if(!gMinuit) gMinuit = new TMinuit(100) ;
316
317 // add Task to //root/Tasks folder
318 TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ;
319 roottasks->Add(this) ;
320
321 fIsInitialized = kTRUE ;
322 }
323}
324//____________________________________________________________________________
325Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
d15a28e7 326{
b2a60966 327 // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
328 // = 1 are neighbour
329 // = 2 are not neighbour but do not continue searching
9a1398dd 330 // neighbours are defined as digits having at least common vertex
b2a60966 331 // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
332 // which is compared to a digit (d2) not yet in a cluster
333
d15a28e7 334 Int_t rv = 0 ;
335
d15a28e7 336 Int_t relid1[4] ;
6b9b768e 337 fGeom->AbsToRelNumbering(d1->GetId(), relid1) ;
d15a28e7 338
339 Int_t relid2[4] ;
6b9b768e 340 fGeom->AbsToRelNumbering(d2->GetId(), relid2) ;
d15a28e7 341
342 if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module and the same PPSD Module
92862013 343 Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
344 Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
d15a28e7 345
92862013 346 if (( coldiff <= 1 ) && ( rowdiff <= 1 )){
d15a28e7 347 rv = 1 ;
348 }
349 else {
350 if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
351 rv = 2; // Difference in row numbers is too large to look further
352 }
353
354 }
355 else {
356
357 if( (relid1[0] < relid2[0]) || (relid1[1] < relid2[1]) )
358 rv=2 ;
359
360 }
d72dfbc3 361
362 //Do NOT clusterize upper PPSD
ed4205d8 363 if( IsInPpsd(d1) && IsInPpsd(d2) &&
364 relid1[1] > 0 &&
d72dfbc3 365 relid1[1] < fGeom->GetNumberOfPadsPhi()*fGeom->GetNumberOfPadsPhi() ) rv = 2 ;
366
d15a28e7 367 return rv ;
368}
369
d15a28e7 370
371//____________________________________________________________________________
9a1398dd 372Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
d15a28e7 373{
b2a60966 374 // Tells if (true) or not (false) the digit is in a PHOS-EMC module
375
9f616d61 376 Bool_t rv = kFALSE ;
d15a28e7 377
d15a28e7 378 Int_t relid[4] ;
6b9b768e 379 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
d15a28e7 380
ed4205d8 381 if ( relid[1] == 0 ) rv = kTRUE;
382
383 return rv ;
384}
385
386//____________________________________________________________________________
9a1398dd 387Bool_t AliPHOSClusterizerv1::IsInPpsd(AliPHOSDigit * digit) const
ed4205d8 388{
fad3e5b9 389 // Tells if (true) or not (false) the digit is in a PHOS-PPSD module
ed4205d8 390
391 Bool_t rv = kFALSE ;
392
393 Int_t relid[4] ;
394 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
395
396 if ( relid[1] > 0 && relid[0] > fGeom->GetNCPVModules() ) rv = kTRUE;
397
398 return rv ;
399}
400
401//____________________________________________________________________________
9a1398dd 402Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
ed4205d8 403{
fad3e5b9 404 // Tells if (true) or not (false) the digit is in a PHOS-CPV module
ed4205d8 405
406 Bool_t rv = kFALSE ;
407
408 Int_t relid[4] ;
409 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
410
411 if ( relid[1] > 0 && relid[0] <= fGeom->GetNCPVModules() ) rv = kTRUE;
d15a28e7 412
413 return rv ;
414}
d15a28e7 415//____________________________________________________________________________
9a1398dd 416Bool_t AliPHOSClusterizerv1::ReadDigits(){
f035f6ce 417 //reads digitis with specified title from TreeD
c161df70 418
b73f246d 419 fNumberOfEmcClusters = 0 ;
b73f246d 420 fNumberOfCpvClusters = 0 ;
c161df70 421
9a1398dd 422 // Get Digits Tree header from file
9a1398dd 423 gAlice->GetEvent(fEvent) ;
424 gAlice->SetEvent(fEvent) ;
425
f035f6ce 426 TTree * treeD = gAlice->TreeD() ;
9a1398dd 427
428 if(treeD==0){
f035f6ce 429 char treeName[20];
430 sprintf(treeName,"TreeD%d",fEvent);
9a1398dd 431 cout << "Error in AliPHOSClusterizerv1 : no "<<treeName << endl ;
432 cout << " Do nothing " << endl ;
433 return kFALSE ;
434 }
435
436 TBranch * digitsBranch = 0;
437 TBranch * digitizerBranch = 0;
438
439 TObjArray * branches = treeD->GetListOfBranches() ;
440 Int_t ibranch;
441 Bool_t phosNotFound = kTRUE ;
442 Bool_t digitizerNotFound = kTRUE ;
443
444 for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
445
446 if(phosNotFound){
447 digitsBranch=(TBranch *) branches->At(ibranch) ;
448 if( fDigitsBranchTitle.CompareTo(digitsBranch->GetTitle())==0 )
449 if( strcmp(digitsBranch->GetName(),"PHOS") == 0)
450 phosNotFound = kFALSE ;
451 }
452
453 if(digitizerNotFound){
454 digitizerBranch = (TBranch *) branches->At(ibranch) ;
455 if( fDigitsBranchTitle.CompareTo(digitizerBranch->GetTitle()) == 0)
456 if( strcmp(digitizerBranch->GetName(),"AliPHOSDigitizer") == 0)
457 digitizerNotFound = kFALSE ;
458 }
459
460 }
461
462 if(digitizerNotFound || phosNotFound){
463 cout << "ERROR in AliPHOSClusterizerv1: " << endl ;
464 cout << " Can't find Branch with digits or Digitizer "<< endl ; ;
465 cout << " Do nothing" <<endl ;
466 return kFALSE ;
467 }
468
469 digitsBranch->SetAddress(&fDigits) ;
470 digitizerBranch->SetAddress(&fDigitizer) ;
761e34c0 471
472 digitsBranch->GetEntry(0) ;
473 digitizerBranch->GetEntry(0) ;
9a1398dd 474
475 fPedestal = fDigitizer->GetPedestal() ;
476 fSlope = fDigitizer->GetSlope() ;
477 return kTRUE ;
478}
479
480//____________________________________________________________________________
481void AliPHOSClusterizerv1::WriteRecPoints(){
f035f6ce 482 // checks, if PHOSEmcRP etc. branches with given title already exist,
483 // exits without writing, otherwise create new branches with given title
484 // fills and wrights TreeR.
9a1398dd 485
486 Int_t index ;
487 //Evaluate poisition, dispersion and other RecPoint properties...
488 for(index = 0; index < fEmcRecPoints->GetEntries(); index++)
489 ((AliPHOSEmcRecPoint *)fEmcRecPoints->At(index))->EvalAll(fW0,fDigits) ;
490
491 fEmcRecPoints->Sort() ;
492
493 for(index = 0; index < fEmcRecPoints->GetEntries(); index++)
494 ((AliPHOSEmcRecPoint *)fEmcRecPoints->At(index))->SetIndexInList(index) ;
495
761e34c0 496 fEmcRecPoints->Expand(fEmcRecPoints->GetEntriesFast()) ;
497
9a1398dd 498 //Now the same for CPV
499 for(index = 0; index < fCpvRecPoints->GetEntries(); index++)
500 ((AliPHOSRecPoint *)fCpvRecPoints->At(index))->EvalAll(fW0CPV,fDigits) ;
501
502 fCpvRecPoints->Sort() ;
503
504 for(index = 0; index < fCpvRecPoints->GetEntries(); index++)
505 ((AliPHOSRecPoint *)fCpvRecPoints->At(index))->SetIndexInList(index) ;
506
761e34c0 507 fCpvRecPoints->Expand(fCpvRecPoints->GetEntriesFast()) ;
508
9a1398dd 509 if(gAlice->TreeR()==0)
510 gAlice->MakeTree("R") ;
511
512
513 //Check, if branches already exist
514 TBranch * emcBranch = 0;
515 TBranch * cpvBranch = 0;
516 TBranch * clusterizerBranch = 0;
517
518 TObjArray * branches = gAlice->TreeR()->GetListOfBranches() ;
519 Int_t ibranch;
520 Bool_t emcNotFound = kTRUE ;
521 Bool_t cpvNotFound = kTRUE ;
522 Bool_t clusterizerNotFound = kTRUE ;
523
524 for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
525
526 if(emcNotFound){
527 emcBranch=(TBranch *) branches->At(ibranch) ;
528 if( fRecPointsBranchTitle.CompareTo(emcBranch->GetTitle())==0 ){
529 if( strcmp(emcBranch->GetName(),"PHOSEmcRP") == 0) {
530 emcNotFound = kFALSE ;
531 }
532 }
533 }
534 if(cpvNotFound){
535 cpvBranch=(TBranch *) branches->At(ibranch) ;
536 if( fRecPointsBranchTitle.CompareTo(cpvBranch->GetTitle())==0 ){
537 if( strcmp(cpvBranch->GetName(),"PHOSCpvRP") == 0) {
538 cpvNotFound = kFALSE ;
539 }
540 }
541 }
542
543 if(clusterizerNotFound){
544 clusterizerBranch = (TBranch *) branches->At(ibranch) ;
545 if( fRecPointsBranchTitle.CompareTo(clusterizerBranch->GetTitle()) == 0){
546 if( strcmp(clusterizerBranch->GetName(),"AliPHOSClusterizer") == 0) {
547 clusterizerNotFound = kFALSE ;
548 }
549 }
550 }
551
552 }
553
554 if(!(clusterizerNotFound && emcNotFound && cpvNotFound)){
555 cout << "AliPHOSClusterizer error" << endl;
556 cout << " Branches PHOSEmcRP, PHOSCpvRP and AliPHOSClusterizer " << endl ;
557 cout << " with title '" << fRecPointsBranchTitle.Data() <<"' already exist" << endl ;
558 cout << " can not overwrite " << endl ;
559 return ;
560 }
561
562 //Make branches in TreeR for RecPoints and Clusterizer
563 char * filename = 0;
564 if(gSystem->Getenv("CONFIG_SPLIT_FILE")!=0){ //generating file name
565 filename = new char[strlen(gAlice->GetBaseFile())+20] ;
566 sprintf(filename,"%s/PHOS.Reco.root",gAlice->GetBaseFile()) ;
567 }
568
569 //Make new branches
570 TDirectory *cwd = gDirectory;
571
572 //First EMC
573 Int_t bufferSize = 32000 ;
574 Int_t splitlevel = 0 ;
575 emcBranch = gAlice->TreeR()->Branch("PHOSEmcRP","TObjArray",&fEmcRecPoints,bufferSize,splitlevel);
576 emcBranch->SetTitle(fRecPointsBranchTitle.Data());
577 if (filename) {
578 emcBranch->SetFile(filename);
579 TIter next( emcBranch->GetListOfBranches());
761e34c0 580 TBranch * sb ;
581 while ((sb=(TBranch*)next())) {
582 sb->SetFile(filename);
9a1398dd 583 }
584 cwd->cd();
585 }
586
587 //Now CPV branch
588 cpvBranch = gAlice->TreeR()->Branch("PHOSCpvRP","TObjArray",&fCpvRecPoints,bufferSize,splitlevel);
589 cpvBranch->SetTitle(fRecPointsBranchTitle.Data());
590 if (filename) {
591 cpvBranch->SetFile(filename);
592 TIter next( cpvBranch->GetListOfBranches());
761e34c0 593 TBranch * sb;
594 while ((sb=(TBranch*)next())) {
595 sb->SetFile(filename);
9a1398dd 596 }
597 cwd->cd();
598 }
599
600 //And Finally clusterizer branch
601 AliPHOSClusterizerv1 * cl = this ;
602 clusterizerBranch = gAlice->TreeR()->Branch("AliPHOSClusterizer","AliPHOSClusterizerv1",
603 &cl,bufferSize,splitlevel);
604 clusterizerBranch->SetTitle(fRecPointsBranchTitle.Data());
605 if (filename) {
606 clusterizerBranch->SetFile(filename);
607 TIter next( clusterizerBranch->GetListOfBranches());
761e34c0 608 TBranch * sb ;
609 while ((sb=(TBranch*)next())) {
610 sb->SetFile(filename);
9a1398dd 611 }
612 cwd->cd();
613 }
614
9a1398dd 615
761e34c0 616 emcBranch->Fill() ;
617 cpvBranch->Fill() ;
618 clusterizerBranch->Fill() ;
9a1398dd 619 gAlice->TreeR()->Write(0,kOverwrite) ;
620
621}
622
623//____________________________________________________________________________
624void AliPHOSClusterizerv1::MakeClusters()
625{
626 // Steering method to construct the clusters stored in a list of Reconstructed Points
627 // A cluster is defined as a list of neighbour digits
628 fEmcRecPoints->Clear() ;
629 fCpvRecPoints->Clear() ;
630
d15a28e7 631
632 // Clusterization starts
9a1398dd 633 TClonesArray * digits = (TClonesArray*)fDigits->Clone() ;
634
635 TIter nextdigit(digits) ;
d15a28e7 636 AliPHOSDigit * digit ;
92862013 637 Bool_t notremoved = kTRUE ;
6ad0bfa0 638
d15a28e7 639 while ( (digit = (AliPHOSDigit *)nextdigit()) ) { // scan over the list of digits
9a1398dd 640 AliPHOSRecPoint * clu = 0 ;
c161df70 641
9a1398dd 642 TArrayI clusterdigitslist(1000) ;
d15a28e7 643 Int_t index ;
f2bc1b87 644
b73f246d 645 if (( IsInEmc (digit) && Calibrate(digit->GetAmp()) > fEmcClusteringThreshold ) ||
646 ( IsInPpsd(digit) && Calibrate(digit->GetAmp()) > fPpsdClusteringThreshold ) ||
647 ( IsInCpv (digit) && Calibrate(digit->GetAmp()) > fCpvClusteringThreshold ) ) {
f2bc1b87 648
d15a28e7 649 Int_t iDigitInCluster = 0 ;
650
6ad0bfa0 651 if ( IsInEmc(digit) ) {
83974468 652 // start a new EMC RecPoint
9a1398dd 653 if(fNumberOfEmcClusters >= fEmcRecPoints->GetSize()) fEmcRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
7f5d510c 654 fEmcRecPoints->AddAt(new AliPHOSEmcRecPoint(), fNumberOfEmcClusters) ;
9a1398dd 655 clu = (AliPHOSEmcRecPoint *) fEmcRecPoints->At(fNumberOfEmcClusters) ;
7f5d510c 656 fNumberOfEmcClusters++ ;
9f616d61 657 clu->AddDigit(*digit, Calibrate(digit->GetAmp())) ;
9a1398dd 658 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
9f616d61 659 iDigitInCluster++ ;
9a1398dd 660 digits->Remove(digit) ;
ca77b032 661
d72dfbc3 662 } else {
83974468 663
fad3e5b9 664 // start a new PPSD/CPV cluster
9a1398dd 665 if(fNumberOfCpvClusters >= fCpvRecPoints->GetSize()) fCpvRecPoints->Expand(2*fNumberOfCpvClusters+1);
666 if(IsInPpsd(digit))
667 fCpvRecPoints->AddAt(new AliPHOSPpsdRecPoint(),fNumberOfCpvClusters) ;
fad3e5b9 668 else
9a1398dd 669 fCpvRecPoints->AddAt(new AliPHOSCpvRecPoint(), fNumberOfCpvClusters) ;
670 clu = (AliPHOSPpsdRecPoint *) fCpvRecPoints->At(fNumberOfCpvClusters) ;
671 fNumberOfCpvClusters++ ;
fad3e5b9 672
6ad0bfa0 673 clu->AddDigit(*digit, Calibrate(digit->GetAmp()) ) ;
9a1398dd 674 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
9f616d61 675 iDigitInCluster++ ;
9a1398dd 676 digits->Remove(digit) ;
d15a28e7 677 nextdigit.Reset() ;
678
9f616d61 679 // Here we remove resting EMC digits, which cannot make cluster
fad3e5b9 680
92862013 681 if( notremoved ) {
9f616d61 682 while( ( digit = (AliPHOSDigit *)nextdigit() ) ) {
9f616d61 683 if( IsInEmc(digit) )
9a1398dd 684 digits->Remove(digit) ;
d15a28e7 685 else
686 break ;
d72dfbc3 687 }
688 notremoved = kFALSE ;
689 }
d15a28e7 690
691 } // else
692
693 nextdigit.Reset() ;
fad3e5b9 694
d15a28e7 695 AliPHOSDigit * digitN ;
696 index = 0 ;
9f616d61 697 while (index < iDigitInCluster){ // scan over digits already in cluster
9a1398dd 698 digit = (AliPHOSDigit*)fDigits->At(clusterdigitslist[index]) ;
9f616d61 699 index++ ;
d15a28e7 700 while ( (digitN = (AliPHOSDigit *)nextdigit()) ) { // scan over the reduced list of digits
d72dfbc3 701 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
d15a28e7 702 switch (ineb ) {
9f616d61 703 case 0 : // not a neighbour
d72dfbc3 704 break ;
9f616d61 705 case 1 : // are neighbours
83974468 706 clu->AddDigit(*digitN, Calibrate( digitN->GetAmp() ) ) ;
9a1398dd 707 clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ;
9f616d61 708 iDigitInCluster++ ;
9a1398dd 709 digits->Remove(digitN) ;
d15a28e7 710 break ;
9f616d61 711 case 2 : // too far from each other
d15a28e7 712 goto endofloop;
713 } // switch
714
715 } // while digitN
fad3e5b9 716
d15a28e7 717 endofloop: ;
718 nextdigit.Reset() ;
719
720 } // loop over cluster
ad8cfaf4 721
ad8cfaf4 722 } // energy theshold
ad8cfaf4 723
31aa6d6c 724
d15a28e7 725 } // while digit
726
9a1398dd 727 delete digits ;
728
729}
730
731//____________________________________________________________________________
732void AliPHOSClusterizerv1::MakeUnfolding(){
733 //Unfolds clusters using the shape of ElectroMagnetic shower
734 // Performs unfolding of all EMC/CPV but NOT ppsd clusters
735
736 //Unfold first EMC clusters
737 if(fNumberOfEmcClusters > 0){
738
739 Int_t nModulesToUnfold = fGeom->GetNModules() ;
740
741 Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
742 Int_t index ;
743 for(index = 0 ; index < numberofNotUnfolded ; index++){
744
745 AliPHOSEmcRecPoint * emcRecPoint = (AliPHOSEmcRecPoint *) fEmcRecPoints->At(index) ;
746 if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
747 break ;
748
749 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
750 Int_t * maxAt = new Int_t[nMultipl] ;
751 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
752 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,fDigits) ;
753
754 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
755 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
756 fEmcRecPoints->Remove(emcRecPoint);
757 fEmcRecPoints->Compress() ;
758 index-- ;
759 fNumberOfEmcClusters -- ;
760 numberofNotUnfolded-- ;
761 }
762
763 delete[] maxAt ;
764 delete[] maxAtEnergy ;
765 }
766 }
767 //Unfolding of EMC clusters finished
768
769
770 //Unfold now CPV clusters
771 if(fNumberOfCpvClusters > 0){
772
773 Int_t nModulesToUnfold = fGeom->GetNCPVModules() ;
774
775 Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
776 Int_t index ;
777 for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
778
779 AliPHOSRecPoint * recPoint = (AliPHOSRecPoint *) fCpvRecPoints->At(index) ;
780
781 if(recPoint->GetPHOSMod()> nModulesToUnfold)
782 break ;
783
784 AliPHOSEmcRecPoint * emcRecPoint = (AliPHOSEmcRecPoint*) recPoint ;
785
786 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
787 Int_t * maxAt = new Int_t[nMultipl] ;
788 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
789 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,fDigits) ;
790
791 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
792 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
793 fCpvRecPoints->Remove(emcRecPoint);
794 fCpvRecPoints->Compress() ;
795 index-- ;
796 numberofCpvNotUnfolded-- ;
797 fNumberOfCpvClusters-- ;
798 }
799
800 delete[] maxAt ;
801 delete[] maxAtEnergy ;
802 }
803 }
804 //Unfolding of Cpv clusters finished
805
806}
807
9a1398dd 808//____________________________________________________________________________
809Double_t AliPHOSClusterizerv1::ShowerShape(Double_t r)
810{
811 // Shape of the shower (see PHOS TDR)
812 // If you change this function, change also the gradien evaluation in ChiSquare()
813
814 Double_t r4 = r*r*r*r ;
815 Double_t r295 = TMath::Power(r, 2.95) ;
816 Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
817 return shape ;
818}
819
820//____________________________________________________________________________
821void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
822 Int_t nMax,
823 int * maxAt,
824 Float_t * maxAtEnergy)
825{
826 // Performs the unfolding of a cluster with nMax overlapping showers
827
828 Int_t nPar = 3 * nMax ;
829 Float_t * fitparameters = new Float_t[nPar] ;
830
831 Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
832 if( !rv ) {
833 // Fit failed, return and remove cluster
834 delete[] fitparameters ;
835 return ;
836 }
837
838 // create ufolded rec points and fill them with new energy lists
839 // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
840 // and later correct this number in acordance with actual energy deposition
841
842 Int_t nDigits = iniEmc->GetMultiplicity() ;
843 Float_t * efit = new Float_t[nDigits] ;
844 Float_t xDigit,zDigit,distance ;
845 Float_t xpar,zpar,epar ;
846 Int_t relid[4] ;
847 AliPHOSDigit * digit ;
848 Int_t * emcDigits = iniEmc->GetDigitsList() ;
849
850 Int_t iparam ;
851 Int_t iDigit ;
852 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
853 digit = (AliPHOSDigit*) fDigits->At(emcDigits[iDigit] ) ;
854 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
855 fGeom->RelPosInModule(relid, xDigit, zDigit) ;
856 efit[iDigit] = 0;
857
858 iparam = 0 ;
859 while(iparam < nPar ){
860 xpar = fitparameters[iparam] ;
861 zpar = fitparameters[iparam+1] ;
862 epar = fitparameters[iparam+2] ;
863 iparam += 3 ;
864 distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ;
865 distance = TMath::Sqrt(distance) ;
866 efit[iDigit] += epar * ShowerShape(distance) ;
867 }
868 }
869
870
871 // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
872 // so that energy deposited in each cell is distributed betwin new clusters proportionally
873 // to its contribution to efit
874
875 Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
876 Float_t ratio ;
877
878 iparam = 0 ;
879 while(iparam < nPar ){
880 xpar = fitparameters[iparam] ;
881 zpar = fitparameters[iparam+1] ;
882 epar = fitparameters[iparam+2] ;
883 iparam += 3 ;
884
885 AliPHOSEmcRecPoint * emcRP ;
886
887 if(iniEmc->IsEmc()){ //create new entries in fEmcRecPoints...
888
889 if(fNumberOfEmcClusters >= fEmcRecPoints->GetSize())
890 fEmcRecPoints->Expand(2*fNumberOfEmcClusters) ;
891
892 (*fEmcRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint() ;
893 emcRP = (AliPHOSEmcRecPoint *) fEmcRecPoints->At(fNumberOfEmcClusters);
894 fNumberOfEmcClusters++ ;
895 }
896 else{//create new entries in fCpvRecPoints
897 if(fNumberOfCpvClusters >= fCpvRecPoints->GetSize())
898 fCpvRecPoints->Expand(2*fNumberOfCpvClusters) ;
899
900 (*fCpvRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint() ;
901 emcRP = (AliPHOSEmcRecPoint *) fCpvRecPoints->At(fNumberOfCpvClusters);
902 fNumberOfCpvClusters++ ;
903 }
904
905 Float_t eDigit ;
906 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
907 digit = (AliPHOSDigit*) fDigits->At( emcDigits[iDigit] ) ;
908 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
909 fGeom->RelPosInModule(relid, xDigit, zDigit) ;
910 distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ;
911 distance = TMath::Sqrt(distance) ;
912 ratio = epar * ShowerShape(distance) / efit[iDigit] ;
913 eDigit = emcEnergies[iDigit] * ratio ;
914 emcRP->AddDigit( *digit, eDigit ) ;
915 }
916 }
917
918 delete[] fitparameters ;
919 delete[] efit ;
920
921}
922
923//_____________________________________________________________________________
924void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
925{
926 // Calculates th Chi square for the cluster unfolding minimization
927 // Number of parameters, Gradient, Chi squared, parameters, what to do
928
929
930 TList * toMinuit = (TList*) gMinuit->GetObjectFit() ;
931
932 AliPHOSEmcRecPoint * emcRP = (AliPHOSEmcRecPoint*) toMinuit->At(0) ;
933 TClonesArray * digits = (TClonesArray*)toMinuit->At(1) ;
934
935
936
937 // AliPHOSEmcRecPoint * emcRP = (AliPHOSEmcRecPoint *) gMinuit->GetObjectFit() ; // EmcRecPoint to fit
938
939 Int_t * emcDigits = emcRP->GetDigitsList() ;
940
941 Int_t nOfDigits = emcRP->GetDigitsMultiplicity() ;
942
943 Float_t * emcEnergies = emcRP->GetEnergiesList() ;
944
945 AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
946
947 fret = 0. ;
948 Int_t iparam ;
949
950 if(iflag == 2)
951 for(iparam = 0 ; iparam < nPar ; iparam++)
952 Grad[iparam] = 0 ; // Will evaluate gradient
953
954 Double_t efit ;
955
956 AliPHOSDigit * digit ;
957 Int_t iDigit ;
958
959 for( iDigit = 0 ; iDigit < nOfDigits ; iDigit++) {
960
961 digit = (AliPHOSDigit*) digits->At( emcDigits[iDigit] ) ;
962
963 Int_t relid[4] ;
964 Float_t xDigit ;
965 Float_t zDigit ;
966
967 geom->AbsToRelNumbering(digit->GetId(), relid) ;
968
969 geom->RelPosInModule(relid, xDigit, zDigit) ;
970
971 if(iflag == 2){ // calculate gradient
972 Int_t iParam = 0 ;
973 efit = 0 ;
974 while(iParam < nPar ){
975 Double_t distance = (xDigit - x[iParam]) * (xDigit - x[iParam]) ;
976 iParam++ ;
977 distance += (zDigit - x[iParam]) * (zDigit - x[iParam]) ;
978 distance = TMath::Sqrt( distance ) ;
979 iParam++ ;
980 efit += x[iParam] * ShowerShape(distance) ;
981 iParam++ ;
982 }
983 Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
984 iParam = 0 ;
985 while(iParam < nPar ){
986 Double_t xpar = x[iParam] ;
987 Double_t zpar = x[iParam+1] ;
988 Double_t epar = x[iParam+2] ;
989 Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
990 Double_t shape = sum * ShowerShape(dr) ;
991 Double_t r4 = dr*dr*dr*dr ;
992 Double_t r295 = TMath::Power(dr,2.95) ;
993 Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
994 0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
995
996 Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x
997 iParam++ ;
998 Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z
999 iParam++ ;
1000 Grad[iParam] += shape ; // Derivative over energy
1001 iParam++ ;
1002 }
1003 }
1004 efit = 0;
1005 iparam = 0 ;
1006
1007 while(iparam < nPar ){
1008 Double_t xpar = x[iparam] ;
1009 Double_t zpar = x[iparam+1] ;
1010 Double_t epar = x[iparam+2] ;
1011 iparam += 3 ;
1012 Double_t distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ;
1013 distance = TMath::Sqrt(distance) ;
1014 efit += epar * ShowerShape(distance) ;
1015 }
1016
1017 fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
1018 // Here we assume, that sigma = sqrt(E)
1019 }
83974468 1020
d15a28e7 1021}
1022
1023//____________________________________________________________________________
9a1398dd 1024void AliPHOSClusterizerv1::Print(Option_t * option)const
d15a28e7 1025{
9a1398dd 1026 if(fIsInitialized){
1027
1028 // Print parameters
1029
1030 cout << "---------------"<< GetName() << " " << GetTitle()<< "-----------" << endl
1031 << "Clusterizing digits from the file: " << fHeaderFileName.Data() << endl
1032 << " Branch: " << fDigitsBranchTitle.Data() << endl
1033 << endl
1034 << " EMC Clustering threshold = " << fEmcClusteringThreshold << endl
1035 << " EMC Local Maximum cut = " << fEmcLocMaxCut << endl
1036 << " EMC Logarothmic weight = " << fW0 << endl
1037 << endl
1038 << " CPV Clustering threshold = " << fCpvClusteringThreshold << endl
1039 << " CPV Local Maximum cut = " << fCpvLocMaxCut << endl
1040 << " CPV Logarothmic weight = " << fW0CPV << endl
1041 << endl
1042 << " PPSD Clustering threshold = " << fPpsdClusteringThreshold << endl;
1043 if(fToUnfold)
1044 cout << " Unfolding on " << endl ;
1045 else
1046 cout << " Unfolding off " << endl ;
1047
1048 cout << "------------------------------------------------------------------" <<endl ;
1049 }
1050 else
1051 cout << " AliPHOSClusterizerv1 not initialized " << endl ;
1052}
1053//____________________________________________________________________________
1054void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option){
1055 //Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
1056
1057 cout << "AliPHOSClusterizerv1: " << endl ;
1058 cout << " Found "<< fEmcRecPoints->GetEntriesFast() << " EMC Rec Points and "
1059 << fCpvRecPoints->GetEntriesFast() << " CPV RecPoints" << endl ;
1060
1061 if(strstr(option,"all")) {
1062 cout << "EMC clusters " << endl ;
1063 cout << " Index "
1064 << " Ene(MeV) "
1065 << " Multi "
1066 << " Module "
1067 << " X "
1068 << " Y "
1069 << " Z "
1070 << " Lambda 1 "
1071 << " Lambda 2 "
1072 << " MaxEnergy "
1073 << " # of prim "
1074 << " Primaries list " << endl;
1075
1076 Int_t index ;
1077 for (index = 0 ; index < fEmcRecPoints->GetEntries() ; index++) {
1078 AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )fEmcRecPoints->At(index) ;
1079
1080 cout << setw(6) << rp->GetIndexInList() << " ";
1081 cout << setw(6) << rp->GetEnergy() << " ";
1082 cout << setw(6) << rp->GetMultiplicity()<< " ";
1083 cout << setw(6) << rp->GetPHOSMod() << " ";
1084
1085 TVector3 locpos;
1086 rp->GetLocalPosition(locpos);
1087 cout << setw(8) << locpos.X() << " ";
1088 cout << setw(8) << locpos.Y() << " ";
1089 cout << setw(8) << locpos.Z() << " ";
1090
1091 Float_t lambda[2];
1092 rp->GetElipsAxis(lambda);
1093 cout << setw(10)<< lambda[0] << " ";
1094 cout << setw(10)<< lambda[1] << " ";
1095
1096
1097 Int_t * primaries;
1098 Int_t nprimaries;
1099 primaries = rp->GetPrimaries(nprimaries);
1100 cout << setw(8) << primaries << " ";
1101
1102 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++)
1103 cout << setw(4) << primaries[iprimary] << " ";
1104 cout << endl;
1105 }
1106 cout << endl ;
1107
1108 //Now plot PCV/PPSD recPoints
1109 cout << "EMC clusters " << endl ;
1110 cout << " Index "
1111 << " Multi "
1112 << " Module "
1113 << " Layer "
1114 << " X "
1115 << " Y "
1116 << " Z "
1117 << " # of prim "
1118 << " Primaries list " << endl;
1119
1120 for (index = 0 ; index < fEmcRecPoints->GetEntries() ; index++) {
1121 AliPHOSRecPoint * rp = (AliPHOSRecPoint * )fCpvRecPoints->At(index) ;
1122 cout << setw(6) << rp->GetIndexInList() << " ";
1123 cout << setw(6) << rp->GetPHOSMod() << " ";
1124
1125 if( (strcmp(rp->ClassName() , "AliPHOSPpsdRecPoint" )) == 0){
1126 AliPHOSPpsdRecPoint * ppsd = (AliPHOSPpsdRecPoint*) rp ;
1127 if(ppsd->GetUp())
1128 cout <<" CPV ";
1129 else
1130 cout <<" PPSD ";
1131 }
1132 else
1133 cout <<" CPV ";
1134
1135 TVector3 locpos;
1136 rp->GetLocalPosition(locpos);
1137 cout << setw(8) << locpos.X() << " ";
1138 cout << setw(8) << locpos.Y() << " ";
1139 cout << setw(8) << locpos.Z() << " ";
1140
1141 Int_t * primaries;
1142 Int_t nprimaries;
1143 primaries = rp->GetPrimaries(nprimaries);
1144 cout << setw(8) << primaries << " ";
1145
1146 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++)
1147 cout << setw(4) << primaries[iprimary] << " ";
1148 cout << endl;
1149 }
1150
1151
1152 cout << "-------------------------------------------------"<<endl ;
1153 }
d15a28e7 1154}
9a1398dd 1155