fCalcMass not copied in ctor. with TParticle parameter
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALClusterizerv1.cxx
CommitLineData
483b0559 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 **************************************************************************/
173558f2 15
483b0559 16/* $Id$ */
17
173558f2 18
483b0559 19/* $Log:
173558f2 20
483b0559 21 1 October 2000. Yuri Kharlov:
173558f2 22
483b0559 23 AreNeighbours()
173558f2 24
483b0559 25 PPSD upper layer is considered if number of layers>1
26
173558f2 27
28
483b0559 29 18 October 2000. Yuri Kharlov:
173558f2 30
483b0559 31 AliEMCALClusterizerv1()
173558f2 32
483b0559 33 CPV clusterizing parameters added
34
173558f2 35
36
483b0559 37 MakeClusters()
173558f2 38
483b0559 39 After first PPSD digit remove EMC digits only once
173558f2 40
483b0559 41*/
173558f2 42
483b0559 43//*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
173558f2 44
05a92d59 45// August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction
173558f2 46
05a92d59 47// of new IO (à la PHOS)
173558f2 48
483b0559 49//////////////////////////////////////////////////////////////////////////////
173558f2 50
483b0559 51// Clusterization class. Performs clusterization (collects neighbouring active cells) and
173558f2 52
483b0559 53// unfolds the clusters having several local maxima.
173558f2 54
483b0559 55// Results are stored in TreeR#, branches EMCALTowerRP (EMC recPoints),
173558f2 56
483b0559 57// EMCALPreShoRP (CPV RecPoints) and AliEMCALClusterizer (Clusterizer with all
173558f2 58
483b0559 59// parameters including input digits branch title, thresholds etc.)
173558f2 60
483b0559 61// This TTask is normally called from Reconstructioner, but can as well be used in
173558f2 62
483b0559 63// standalone mode.
173558f2 64
483b0559 65// Use Case:
173558f2 66
483b0559 67// root [0] AliEMCALClusterizerv1 * cl = new AliEMCALClusterizerv1("galice.root")
173558f2 68
483b0559 69// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
173558f2 70
483b0559 71// //reads gAlice from header file "..."
173558f2 72
483b0559 73// root [1] cl->ExecuteTask()
173558f2 74
483b0559 75// //finds RecPoints in all events stored in galice.root
173558f2 76
483b0559 77// root [2] cl->SetDigitsBranch("digits2")
173558f2 78
483b0559 79// //sets another title for Digitis (input) branch
173558f2 80
483b0559 81// root [3] cl->SetRecPointsBranch("recp2")
173558f2 82
483b0559 83// //sets another title four output branches
173558f2 84
483b0559 85// root [4] cl->SetTowerLocalMaxCut(0.03)
173558f2 86
483b0559 87// //set clusterization parameters
173558f2 88
483b0559 89// root [5] cl->ExecuteTask("deb all time")
173558f2 90
483b0559 91// //once more finds RecPoints options are
173558f2 92
483b0559 93// // deb - print number of found rec points
173558f2 94
483b0559 95// // deb all - print number of found RecPoints and some their characteristics
173558f2 96
483b0559 97// // time - print benchmarking results
98
173558f2 99
100
483b0559 101// --- ROOT system ---
102
173558f2 103
104
483b0559 105#include "TROOT.h"
173558f2 106
483b0559 107#include "TFile.h"
173558f2 108
483b0559 109#include "TFolder.h"
173558f2 110
483b0559 111#include "TMath.h"
173558f2 112
483b0559 113#include "TMinuit.h"
173558f2 114
483b0559 115#include "TTree.h"
173558f2 116
483b0559 117#include "TSystem.h"
173558f2 118
483b0559 119#include "TBenchmark.h"
120
173558f2 121
122
483b0559 123// --- Standard library ---
124
173558f2 125
126
70479d0e 127#include <Riostream.h>
483b0559 128
173558f2 129
130
483b0559 131// --- AliRoot header files ---
132
173558f2 133
134
483b0559 135#include "AliEMCALClusterizerv1.h"
173558f2 136
483b0559 137#include "AliEMCALDigit.h"
173558f2 138
483b0559 139#include "AliEMCALDigitizer.h"
173558f2 140
483b0559 141#include "AliEMCALTowerRecPoint.h"
173558f2 142
483b0559 143#include "AliEMCAL.h"
173558f2 144
483b0559 145#include "AliEMCALGetter.h"
173558f2 146
05a92d59 147#include "AliEMCALGeometry.h"
173558f2 148
483b0559 149#include "AliRun.h"
150
173558f2 151
152
483b0559 153ClassImp(AliEMCALClusterizerv1)
173558f2 154
483b0559 155
173558f2 156
483b0559 157//____________________________________________________________________________
173558f2 158
483b0559 159 AliEMCALClusterizerv1::AliEMCALClusterizerv1() : AliEMCALClusterizer()
173558f2 160
483b0559 161{
173558f2 162
483b0559 163 // default ctor (to be used mainly by Streamer)
173558f2 164
483b0559 165
173558f2 166
839828a6 167 InitParameters() ;
173558f2 168
ef305168 169 fDefaultInit = kTRUE ;
173558f2 170
483b0559 171}
172
173558f2 173
174
483b0559 175//____________________________________________________________________________
173558f2 176
05a92d59 177AliEMCALClusterizerv1::AliEMCALClusterizerv1(const char* headerFile, const char* name, const Bool_t toSplit)
173558f2 178
05a92d59 179:AliEMCALClusterizer(headerFile, name, toSplit)
173558f2 180
483b0559 181{
173558f2 182
483b0559 183 // ctor with the indication of the file where header Tree and digits Tree are stored
173558f2 184
483b0559 185
173558f2 186
839828a6 187 InitParameters() ;
173558f2 188
483b0559 189 Init() ;
173558f2 190
05a92d59 191 fDefaultInit = kFALSE ;
483b0559 192
173558f2 193
194
483b0559 195}
05a92d59 196
173558f2 197
198
483b0559 199//____________________________________________________________________________
173558f2 200
483b0559 201 AliEMCALClusterizerv1::~AliEMCALClusterizerv1()
173558f2 202
483b0559 203{
173558f2 204
ef305168 205 // dtor
173558f2 206
05a92d59 207 fSplitFile = 0 ;
173558f2 208
839828a6 209
173558f2 210
ef305168 211}
212
173558f2 213
214
ef305168 215//____________________________________________________________________________
173558f2 216
ef305168 217const TString AliEMCALClusterizerv1::BranchName() const
173558f2 218
ef305168 219{
173558f2 220
ef305168 221 TString branchName(GetName() ) ;
173558f2 222
ef305168 223 branchName.Remove(branchName.Index(Version())-1) ;
173558f2 224
ef305168 225 return branchName ;
173558f2 226
483b0559 227}
839828a6 228
173558f2 229
230
483b0559 231//____________________________________________________________________________
173558f2 232
483b0559 233Float_t AliEMCALClusterizerv1::Calibrate(Int_t amp, Bool_t inpresho) const
173558f2 234
05a92d59 235{//To be replased later by the method, reading individual parameters from the database
236
173558f2 237
238
483b0559 239 if ( inpresho ) // calibrate as pre shower
173558f2 240
483b0559 241 return -fADCpedestalPreSho + amp * fADCchannelPreSho ;
242
173558f2 243
244
483b0559 245 else //calibrate as tower
173558f2 246
483b0559 247 return -fADCpedestalTower + amp * fADCchannelTower ;
173558f2 248
483b0559 249}
05a92d59 250
173558f2 251
252
483b0559 253//____________________________________________________________________________
173558f2 254
483b0559 255void AliEMCALClusterizerv1::Exec(Option_t * option)
173558f2 256
483b0559 257{
173558f2 258
483b0559 259 // Steering method
260
173558f2 261
262
483b0559 263 if( strcmp(GetName(), "")== 0 )
173558f2 264
483b0559 265 Init() ;
266
173558f2 267
268
483b0559 269 if(strstr(option,"tim"))
173558f2 270
483b0559 271 gBenchmark->Start("EMCALClusterizer");
173558f2 272
483b0559 273
173558f2 274
483b0559 275 if(strstr(option,"print"))
173558f2 276
483b0559 277 Print("") ;
278
173558f2 279
280
483b0559 281 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
173558f2 282
05a92d59 283 if(gime->BranchExists("RecPoints"))
173558f2 284
05a92d59 285 return ;
173558f2 286
05a92d59 287 Int_t nevents = gime->MaxEvent() ;
173558f2 288
483b0559 289 Int_t ievent ;
290
173558f2 291
292
483b0559 293 for(ievent = 0; ievent < nevents; ievent++){
294
173558f2 295
296
05a92d59 297 gime->Event(ievent,"D") ;
298
173558f2 299
300
483b0559 301 if(ievent == 0)
173558f2 302
483b0559 303 GetCalibrationParameters() ;
304
173558f2 305
306
483b0559 307 fNumberOfTowerClusters = fNumberOfPreShoClusters = 0 ;
173558f2 308
05a92d59 309
173558f2 310
483b0559 311 MakeClusters() ;
173558f2 312
483b0559 313
173558f2 314
483b0559 315 if(fToUnfold)
173558f2 316
483b0559 317 MakeUnfolding() ;
318
173558f2 319
320
483b0559 321 WriteRecPoints(ievent) ;
322
173558f2 323
324
483b0559 325 if(strstr(option,"deb"))
173558f2 326
483b0559 327 PrintRecPoints(option) ;
328
173558f2 329
330
483b0559 331 //increment the total number of digits per run
173558f2 332
483b0559 333 fRecPointsInRun += gime->TowerRecPoints()->GetEntriesFast() ;
173558f2 334
483b0559 335 fRecPointsInRun += gime->PreShowerRecPoints()->GetEntriesFast() ;
173558f2 336
483b0559 337 }
173558f2 338
483b0559 339
173558f2 340
483b0559 341 if(strstr(option,"tim")){
173558f2 342
483b0559 343 gBenchmark->Stop("EMCALClusterizer");
173558f2 344
483b0559 345 cout << "AliEMCALClusterizer:" << endl ;
173558f2 346
483b0559 347 cout << " took " << gBenchmark->GetCpuTime("EMCALClusterizer") << " seconds for Clusterizing "
173558f2 348
483b0559 349 << gBenchmark->GetCpuTime("EMCALClusterizer")/nevents << " seconds per event " << endl ;
173558f2 350
483b0559 351 cout << endl ;
173558f2 352
483b0559 353 }
173558f2 354
483b0559 355
173558f2 356
483b0559 357}
358
173558f2 359
360
483b0559 361//____________________________________________________________________________
173558f2 362
a0636361 363Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALTowerRecPoint * emcRP, AliEMCALDigit ** maxAt, Float_t * maxAtEnergy,
173558f2 364
483b0559 365 Int_t nPar, Float_t * fitparameters) const
173558f2 366
483b0559 367{
173558f2 368
483b0559 369 // Calls TMinuit to fit the energy distribution of a cluster with several maxima
173558f2 370
483b0559 371 // The initial values for fitting procedure are set equal to the positions of local maxima.
173558f2 372
483b0559 373 // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
374
173558f2 375
376
483b0559 377 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
173558f2 378
483b0559 379 TClonesArray * digits = gime->Digits() ;
173558f2 380
483b0559 381
382
173558f2 383
384
483b0559 385 gMinuit->mncler(); // Reset Minuit's list of paramters
173558f2 386
483b0559 387 gMinuit->SetPrintLevel(-1) ; // No Printout
173558f2 388
483b0559 389 gMinuit->SetFCN(AliEMCALClusterizerv1::UnfoldingChiSquare) ;
173558f2 390
483b0559 391 // To set the address of the minimization function
392
173558f2 393
394
483b0559 395 TList * toMinuit = new TList();
173558f2 396
483b0559 397 toMinuit->AddAt(emcRP,0) ;
173558f2 398
483b0559 399 toMinuit->AddAt(digits,1) ;
173558f2 400
483b0559 401
173558f2 402
483b0559 403 gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare
404
173558f2 405
406
483b0559 407 // filling initial values for fit parameters
173558f2 408
483b0559 409 AliEMCALDigit * digit ;
410
173558f2 411
412
483b0559 413 Int_t ierflg = 0;
173558f2 414
483b0559 415 Int_t index = 0 ;
173558f2 416
483b0559 417 Int_t nDigits = (Int_t) nPar / 3 ;
418
173558f2 419
420
483b0559 421 Int_t iDigit ;
422
173558f2 423
424
483b0559 425 AliEMCALGeometry * geom = gime->EMCALGeometry() ;
426
173558f2 427
428
483b0559 429 for(iDigit = 0; iDigit < nDigits; iDigit++){
173558f2 430
a0636361 431 digit = maxAt[iDigit];
483b0559 432
173558f2 433
434
483b0559 435 Int_t relid[4] ;
173558f2 436
483b0559 437 Float_t x = 0.;
173558f2 438
483b0559 439 Float_t z = 0.;
173558f2 440
483b0559 441 geom->AbsToRelNumbering(digit->GetId(), relid) ;
173558f2 442
483b0559 443 geom->PosInAlice(relid, x, z) ;
444
173558f2 445
446
483b0559 447 Float_t energy = maxAtEnergy[iDigit] ;
448
173558f2 449
450
483b0559 451 gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ;
173558f2 452
483b0559 453 index++ ;
173558f2 454
483b0559 455 if(ierflg != 0){
173558f2 456
483b0559 457 cout << "EMCAL Unfolding> Unable to set initial value for fit procedure : x = " << x << endl ;
173558f2 458
483b0559 459 return kFALSE;
173558f2 460
483b0559 461 }
173558f2 462
483b0559 463 gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ;
173558f2 464
483b0559 465 index++ ;
173558f2 466
483b0559 467 if(ierflg != 0){
173558f2 468
483b0559 469 cout << "EMCAL Unfolding> Unable to set initial value for fit procedure : z = " << z << endl ;
173558f2 470
483b0559 471 return kFALSE;
173558f2 472
483b0559 473 }
173558f2 474
483b0559 475 gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;
173558f2 476
483b0559 477 index++ ;
173558f2 478
483b0559 479 if(ierflg != 0){
173558f2 480
483b0559 481 cout << "EMCAL Unfolding> Unable to set initial value for fit procedure : energy = " << energy << endl ;
173558f2 482
483b0559 483 return kFALSE;
173558f2 484
483b0559 485 }
173558f2 486
483b0559 487 }
488
173558f2 489
490
483b0559 491 Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
173558f2 492
483b0559 493 // depends on it.
173558f2 494
483b0559 495 Double_t p1 = 1.0 ;
173558f2 496
483b0559 497 Double_t p2 = 0.0 ;
498
173558f2 499
500
483b0559 501 gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
173558f2 502
483b0559 503 gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
173558f2 504
483b0559 505 gMinuit->SetMaxIterations(5);
173558f2 506
483b0559 507 gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
508
173558f2 509
510
483b0559 511 gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
512
173558f2 513
514
483b0559 515 if(ierflg == 4){ // Minimum not found
173558f2 516
483b0559 517 cout << "EMCAL Unfolding> Fit not converged, cluster abandoned "<< endl ;
173558f2 518
483b0559 519 return kFALSE ;
173558f2 520
483b0559 521 }
173558f2 522
483b0559 523 for(index = 0; index < nPar; index++){
173558f2 524
483b0559 525 Double_t err ;
173558f2 526
483b0559 527 Double_t val ;
173558f2 528
483b0559 529 gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index
173558f2 530
483b0559 531 fitparameters[index] = val ;
173558f2 532
483b0559 533 }
534
173558f2 535
536
483b0559 537 delete toMinuit ;
173558f2 538
483b0559 539 return kTRUE;
540
173558f2 541
542
483b0559 543}
544
173558f2 545
546
483b0559 547//____________________________________________________________________________
173558f2 548
483b0559 549void AliEMCALClusterizerv1::GetCalibrationParameters()
173558f2 550
483b0559 551{
173558f2 552
483b0559 553 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
173558f2 554
05a92d59 555 const AliEMCALDigitizer * dig = gime->Digitizer(BranchName()) ;
483b0559 556
173558f2 557
558
483b0559 559 fADCchannelTower = dig->GetTowerchannel() ;
173558f2 560
483b0559 561 fADCpedestalTower = dig->GetTowerpedestal();
562
173558f2 563
564
483b0559 565 fADCchannelPreSho = dig->GetPreShochannel() ;
173558f2 566
483b0559 567 fADCpedestalPreSho = dig->GetPreShopedestal() ;
568
173558f2 569
570
483b0559 571}
05a92d59 572
173558f2 573
574
483b0559 575//____________________________________________________________________________
173558f2 576
483b0559 577void AliEMCALClusterizerv1::Init()
173558f2 578
483b0559 579{
173558f2 580
483b0559 581 // Make all memory allocations which can not be done in default constructor.
173558f2 582
483b0559 583 // Attach the Clusterizer task to the list of EMCAL tasks
173558f2 584
483b0559 585
173558f2 586
483b0559 587 if ( strcmp(GetTitle(), "") == 0 )
173558f2 588
483b0559 589 SetTitle("galice.root") ;
590
173558f2 591
592
483b0559 593 TString branchname = GetName() ;
173558f2 594
483b0559 595 branchname.Remove(branchname.Index(Version())-1) ;
596
173558f2 597
598
05a92d59 599 AliEMCALGetter * gime = AliEMCALGetter::GetInstance(GetTitle(), branchname.Data(), fToSplit ) ;
173558f2 600
483b0559 601 if ( gime == 0 ) {
173558f2 602
483b0559 603 cerr << "ERROR: AliEMCALClusterizerv1::Init -> Could not obtain the Getter object !" << endl ;
173558f2 604
483b0559 605 return ;
173558f2 606
483b0559 607 }
05a92d59 608
173558f2 609
610
05a92d59 611 fSplitFile = 0 ;
173558f2 612
05a92d59 613 if(fToSplit){
173558f2 614
05a92d59 615 // construct the name of the file as /path/EMCAL.SDigits.root
173558f2 616
05a92d59 617 //First - extract full path if necessary
173558f2 618
05a92d59 619 TString fileName(GetTitle()) ;
173558f2 620
05a92d59 621 Ssiz_t islash = fileName.Last('/') ;
173558f2 622
05a92d59 623 if(islash<fileName.Length())
173558f2 624
05a92d59 625 fileName.Remove(islash+1,fileName.Length()) ;
173558f2 626
05a92d59 627 else
173558f2 628
05a92d59 629 fileName="" ;
173558f2 630
05a92d59 631 // Next - append the file name
173558f2 632
05a92d59 633 fileName+="EMCAL.RecData." ;
173558f2 634
05a92d59 635 if((strcmp(branchname.Data(),"Default")!=0)&&(strcmp(branchname.Data(),"")!=0)){
173558f2 636
05a92d59 637 fileName+=branchname ;
173558f2 638
05a92d59 639 fileName+="." ;
173558f2 640
05a92d59 641 }
173558f2 642
05a92d59 643 fileName+="root" ;
173558f2 644
05a92d59 645 // Finally - check if the file already opened or open the file
173558f2 646
05a92d59 647 fSplitFile = static_cast<TFile*>(gROOT->GetFile(fileName.Data()));
173558f2 648
05a92d59 649 if(!fSplitFile)
173558f2 650
05a92d59 651 fSplitFile = TFile::Open(fileName.Data(),"update") ;
173558f2 652
05a92d59 653 }
654
173558f2 655
656
483b0559 657
173558f2 658
483b0559 659 const AliEMCALGeometry * geom = gime->EMCALGeometry() ;
173558f2 660
483b0559 661 fNTowers = geom->GetNZ() * geom->GetNPhi() ;
662
173558f2 663
664
483b0559 665 if(!gMinuit)
173558f2 666
483b0559 667 gMinuit = new TMinuit(100) ;
668
173558f2 669
670
483b0559 671 gime->PostClusterizer(this) ;
173558f2 672
483b0559 673 gime->PostRecPoints(branchname ) ;
173558f2 674
05a92d59 675
173558f2 676
483b0559 677}
678
173558f2 679
680
483b0559 681//____________________________________________________________________________
173558f2 682
839828a6 683void AliEMCALClusterizerv1::InitParameters()
173558f2 684
839828a6 685{
173558f2 686
839828a6 687 fNumberOfPreShoClusters = fNumberOfTowerClusters = 0 ;
688
173558f2 689
690
839828a6 691
173558f2 692
839828a6 693
173558f2 694
839828a6 695 fPreShoClusteringThreshold = 0.0001;
173558f2 696
839828a6 697 fTowerClusteringThreshold = 0.2;
173558f2 698
839828a6 699
173558f2 700
839828a6 701 fTowerLocMaxCut = 0.03 ;
173558f2 702
839828a6 703 fPreShoLocMaxCut = 0.03 ;
173558f2 704
839828a6 705
173558f2 706
839828a6 707 fW0 = 4.5 ;
173558f2 708
839828a6 709 fW0CPV = 4.0 ;
710
173558f2 711
712
839828a6 713 fTimeGate = 1.e-8 ;
173558f2 714
839828a6 715
173558f2 716
839828a6 717 fToUnfold = kFALSE ;
173558f2 718
05a92d59 719
173558f2 720
839828a6 721 TString clusterizerName( GetName()) ;
173558f2 722
839828a6 723 if (clusterizerName.IsNull() )
173558f2 724
839828a6 725 clusterizerName = "Default" ;
173558f2 726
839828a6 727 clusterizerName.Append(":") ;
173558f2 728
839828a6 729 clusterizerName.Append(Version()) ;
173558f2 730
839828a6 731 SetName(clusterizerName) ;
173558f2 732
839828a6 733 fRecPointsInRun = 0 ;
05a92d59 734
173558f2 735
736
839828a6 737}
738
173558f2 739
740
839828a6 741//____________________________________________________________________________
173558f2 742
483b0559 743Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2)const
173558f2 744
483b0559 745{
173558f2 746
483b0559 747 // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
173558f2 748
483b0559 749 // = 1 are neighbour
173558f2 750
483b0559 751 // = 2 are not neighbour but do not continue searching
173558f2 752
483b0559 753 // neighbours are defined as digits having at least a common vertex
173558f2 754
483b0559 755 // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
173558f2 756
483b0559 757 // which is compared to a digit (d2) not yet in a cluster
758
173558f2 759
760
483b0559 761 AliEMCALGeometry * geom = AliEMCALGetter::GetInstance()->EMCALGeometry() ;
762
173558f2 763
764
483b0559 765 Int_t rv = 0 ;
766
173558f2 767
768
483b0559 769 Int_t relid1[4] ;
173558f2 770
483b0559 771 geom->AbsToRelNumbering(d1->GetId(), relid1) ;
772
173558f2 773
774
483b0559 775 Int_t relid2[4] ;
173558f2 776
483b0559 777 geom->AbsToRelNumbering(d2->GetId(), relid2) ;
173558f2 778
483b0559 779
173558f2 780
483b0559 781 if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same EMCAL Arm
173558f2 782
483b0559 783 Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
173558f2 784
483b0559 785 Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
173558f2 786
483b0559 787
173558f2 788
483b0559 789 if (( coldiff <= 1 ) && ( rowdiff <= 1 )){
173558f2 790
483b0559 791 if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fTimeGate))
173558f2 792
483b0559 793 rv = 1 ;
173558f2 794
483b0559 795 }
173558f2 796
483b0559 797 else {
173558f2 798
483b0559 799 if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
173558f2 800
483b0559 801 rv = 2; // Difference in row numbers is too large to look further
173558f2 802
483b0559 803 }
804
173558f2 805
806
483b0559 807 }
173558f2 808
483b0559 809 else {
173558f2 810
483b0559 811
173558f2 812
483b0559 813 if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) )
173558f2 814
483b0559 815 rv=2 ;
816
173558f2 817
818
483b0559 819 }
820
173558f2 821
822
483b0559 823 return rv ;
173558f2 824
483b0559 825}
826
827
173558f2 828
829
830
483b0559 831//____________________________________________________________________________
173558f2 832
483b0559 833Bool_t AliEMCALClusterizerv1::IsInTower(AliEMCALDigit * digit) const
173558f2 834
483b0559 835{
173558f2 836
483b0559 837 // Tells if (true) or not (false) the digit is in a EMCAL-Tower
173558f2 838
483b0559 839
173558f2 840
483b0559 841 Bool_t rv = kFALSE ;
173558f2 842
483b0559 843 if (!digit->IsInPreShower())
173558f2 844
483b0559 845 rv = kTRUE;
173558f2 846
483b0559 847 return rv ;
173558f2 848
483b0559 849}
850
173558f2 851
852
483b0559 853//____________________________________________________________________________
173558f2 854
483b0559 855Bool_t AliEMCALClusterizerv1::IsInPreShower(AliEMCALDigit * digit) const
173558f2 856
483b0559 857{
173558f2 858
483b0559 859 // Tells if (true) or not (false) the digit is in a EMCAL-PreShower
173558f2 860
483b0559 861
173558f2 862
483b0559 863 Bool_t rv = kFALSE ;
173558f2 864
483b0559 865 if (digit->IsInPreShower())
173558f2 866
483b0559 867 rv = kTRUE;
173558f2 868
483b0559 869 return rv ;
173558f2 870
483b0559 871}
872
173558f2 873
874
483b0559 875//____________________________________________________________________________
173558f2 876
483b0559 877void AliEMCALClusterizerv1::WriteRecPoints(Int_t event)
173558f2 878
483b0559 879{
880
173558f2 881
882
483b0559 883 // Creates new branches with given title
173558f2 884
483b0559 885 // fills and writes into TreeR.
886
173558f2 887
888
483b0559 889 AliEMCALGetter *gime = AliEMCALGetter::GetInstance() ;
173558f2 890
05a92d59 891 TObjArray * towerRecPoints = gime->TowerRecPoints() ;
173558f2 892
05a92d59 893 TObjArray * preshoRecPoints = gime->PreShowerRecPoints() ;
173558f2 894
05a92d59 895 TClonesArray * digits = gime->Digits() ;
173558f2 896
839828a6 897 TTree * treeR ;
173558f2 898
05a92d59 899
173558f2 900
05a92d59 901 if(fToSplit){
173558f2 902
05a92d59 903 if(!fSplitFile)
173558f2 904
05a92d59 905 return ;
173558f2 906
05a92d59 907 fSplitFile->cd() ;
173558f2 908
05a92d59 909 TString name("TreeR") ;
173558f2 910
05a92d59 911 name += event ;
173558f2 912
05a92d59 913 treeR = dynamic_cast<TTree*>(fSplitFile->Get(name));
173558f2 914
05a92d59 915 }
173558f2 916
05a92d59 917 else{
173558f2 918
05a92d59 919 treeR = gAlice->TreeR();
173558f2 920
05a92d59 921 }
483b0559 922
173558f2 923
924
05a92d59 925 if(!treeR){
173558f2 926
839828a6 927 gAlice->MakeTree("R", fSplitFile);
173558f2 928
05a92d59 929 treeR = gAlice->TreeR() ;
173558f2 930
05a92d59 931 }
173558f2 932
05a92d59 933
173558f2 934
05a92d59 935 Int_t index ;
173558f2 936
483b0559 937 //Evaluate position, dispersion and other RecPoint properties...
173558f2 938
483b0559 939 for(index = 0; index < towerRecPoints->GetEntries(); index++)
173558f2 940
483b0559 941 (dynamic_cast<AliEMCALTowerRecPoint *>(towerRecPoints->At(index)))->EvalAll(fW0,digits) ;
173558f2 942
05a92d59 943
173558f2 944
483b0559 945 towerRecPoints->Sort() ;
946
173558f2 947
948
483b0559 949 for(index = 0; index < towerRecPoints->GetEntries(); index++)
173558f2 950
483b0559 951 (dynamic_cast<AliEMCALTowerRecPoint *>(towerRecPoints->At(index)))->SetIndexInList(index) ;
952
173558f2 953
954
483b0559 955 towerRecPoints->Expand(towerRecPoints->GetEntriesFast()) ;
956
173558f2 957
958
ef305168 959 //Now the same for pre shower
173558f2 960
483b0559 961 for(index = 0; index < preshoRecPoints->GetEntries(); index++)
173558f2 962
483b0559 963 (dynamic_cast<AliEMCALRecPoint *>(preshoRecPoints->At(index)))->EvalAll(fW0CPV,digits) ;
964
173558f2 965
966
483b0559 967 preshoRecPoints->Sort() ;
968
173558f2 969
970
483b0559 971 for(index = 0; index < preshoRecPoints->GetEntries(); index++)
173558f2 972
483b0559 973 (dynamic_cast<AliEMCALRecPoint *>(preshoRecPoints->At(index)))->SetIndexInList(index) ;
974
173558f2 975
976
483b0559 977 preshoRecPoints->Expand(preshoRecPoints->GetEntriesFast()) ;
173558f2 978
483b0559 979
173558f2 980
483b0559 981 Int_t bufferSize = 32000 ;
173558f2 982
483b0559 983 Int_t splitlevel = 0 ;
984
173558f2 985
986
ef305168 987 //First Tower branch
173558f2 988
05a92d59 989 TBranch * towerBranch = treeR->Branch("EMCALTowerRP","TObjArray",&towerRecPoints,bufferSize,splitlevel);
173558f2 990
05a92d59 991 towerBranch->SetTitle(BranchName());
173558f2 992
05a92d59 993
173558f2 994
ef305168 995 //Now Pre Shower branch
173558f2 996
05a92d59 997 TBranch * preshoBranch = treeR->Branch("EMCALPreShoRP","TObjArray",&preshoRecPoints,bufferSize,splitlevel);
173558f2 998
05a92d59 999 preshoBranch->SetTitle(BranchName());
173558f2 1000
483b0559 1001
173558f2 1002
483b0559 1003 //And Finally clusterizer branch
173558f2 1004
ef305168 1005 AliEMCALClusterizerv1 * cl = (AliEMCALClusterizerv1*)gime->Clusterizer(BranchName()) ;
173558f2 1006
839828a6 1007 TBranch * clusterizerBranch = treeR->Branch("AliEMCALClusterizer","AliEMCALClusterizerv1",
173558f2 1008
483b0559 1009 &cl,bufferSize,splitlevel);
173558f2 1010
ef305168 1011 clusterizerBranch->SetTitle(BranchName());
839828a6 1012
173558f2 1013
1014
05a92d59 1015 towerBranch ->Fill() ;
173558f2 1016
05a92d59 1017 preshoBranch ->Fill() ;
173558f2 1018
483b0559 1019 clusterizerBranch->Fill() ;
1020
173558f2 1021
1022
839828a6 1023 treeR->AutoSave() ; //Write(0,kOverwrite) ;
173558f2 1024
05a92d59 1025 if(gAlice->TreeR()!=treeR)
173558f2 1026
05a92d59 1027 treeR->Delete();
173558f2 1028
483b0559 1029}
1030
173558f2 1031
1032
483b0559 1033//____________________________________________________________________________
173558f2 1034
483b0559 1035void AliEMCALClusterizerv1::MakeClusters()
173558f2 1036
483b0559 1037{
173558f2 1038
483b0559 1039 // Steering method to construct the clusters stored in a list of Reconstructed Points
173558f2 1040
483b0559 1041 // A cluster is defined as a list of neighbour digits
173558f2 1042
05a92d59 1043
173558f2 1044
483b0559 1045 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
173558f2 1046
483b0559 1047
173558f2 1048
05a92d59 1049 TObjArray * towerRecPoints = gime->TowerRecPoints(BranchName()) ;
173558f2 1050
05a92d59 1051 TObjArray * preshoRecPoints = gime->PreShowerRecPoints(BranchName()) ;
173558f2 1052
483b0559 1053 towerRecPoints->Delete() ;
173558f2 1054
483b0559 1055 preshoRecPoints->Delete() ;
173558f2 1056
483b0559 1057
173558f2 1058
05a92d59 1059 TClonesArray * digits = gime->Digits() ;
173558f2 1060
05a92d59 1061 if ( !digits ) {
173558f2 1062
05a92d59 1063 cerr << "ERROR: AliEMCALClusterizerv1::MakeClusters -> Digits with name "
173558f2 1064
05a92d59 1065 << GetName() << " not found ! " << endl ;
173558f2 1066
05a92d59 1067 abort() ;
173558f2 1068
05a92d59 1069 }
173558f2 1070
483b0559 1071 TClonesArray * digitsC = dynamic_cast<TClonesArray*>(digits->Clone()) ;
173558f2 1072
483b0559 1073
173558f2 1074
483b0559 1075
173558f2 1076
483b0559 1077 // Clusterization starts
173558f2 1078
483b0559 1079
173558f2 1080
483b0559 1081 TIter nextdigit(digitsC) ;
173558f2 1082
483b0559 1083 AliEMCALDigit * digit ;
173558f2 1084
483b0559 1085 Bool_t notremoved = kTRUE ;
173558f2 1086
483b0559 1087
173558f2 1088
483b0559 1089 while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // scan over the list of digitsC
173558f2 1090
483b0559 1091 AliEMCALRecPoint * clu = 0 ;
173558f2 1092
483b0559 1093
173558f2 1094
483b0559 1095 TArrayI clusterdigitslist(1500) ;
173558f2 1096
483b0559 1097 Int_t index ;
173558f2 1098
483b0559 1099
173558f2 1100
483b0559 1101 if (( IsInTower (digit) && Calibrate(digit->GetAmp(),digit->IsInPreShower()) > fTowerClusteringThreshold ) ||
173558f2 1102
483b0559 1103 ( IsInPreShower (digit) && Calibrate(digit->GetAmp(),digit->IsInPreShower()) > fPreShoClusteringThreshold ) ) {
173558f2 1104
483b0559 1105
173558f2 1106
483b0559 1107 Int_t iDigitInCluster = 0 ;
173558f2 1108
483b0559 1109
173558f2 1110
483b0559 1111 if ( IsInTower(digit) ) {
173558f2 1112
483b0559 1113 // start a new Tower RecPoint
173558f2 1114
483b0559 1115 if(fNumberOfTowerClusters >= towerRecPoints->GetSize())
173558f2 1116
483b0559 1117 towerRecPoints->Expand(2*fNumberOfTowerClusters+1) ;
173558f2 1118
483b0559 1119
173558f2 1120
483b0559 1121 towerRecPoints->AddAt(new AliEMCALTowerRecPoint(""), fNumberOfTowerClusters) ;
173558f2 1122
483b0559 1123 clu = dynamic_cast<AliEMCALTowerRecPoint *>(towerRecPoints->At(fNumberOfTowerClusters)) ;
173558f2 1124
483b0559 1125 fNumberOfTowerClusters++ ;
173558f2 1126
483b0559 1127 clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->IsInPreShower())) ;
173558f2 1128
483b0559 1129 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
173558f2 1130
483b0559 1131 iDigitInCluster++ ;
173558f2 1132
483b0559 1133 digitsC->Remove(digit) ;
173558f2 1134
483b0559 1135
173558f2 1136
483b0559 1137 } else {
173558f2 1138
483b0559 1139
173558f2 1140
483b0559 1141 // start a new Pre Shower cluster
173558f2 1142
483b0559 1143 if(fNumberOfPreShoClusters >= preshoRecPoints->GetSize())
173558f2 1144
483b0559 1145 preshoRecPoints->Expand(2*fNumberOfPreShoClusters+1);
173558f2 1146
483b0559 1147
173558f2 1148
483b0559 1149 preshoRecPoints->AddAt(new AliEMCALTowerRecPoint(""), fNumberOfPreShoClusters) ;
173558f2 1150
483b0559 1151
173558f2 1152
483b0559 1153 clu = dynamic_cast<AliEMCALTowerRecPoint *>(preshoRecPoints->At(fNumberOfPreShoClusters)) ;
173558f2 1154
483b0559 1155 fNumberOfPreShoClusters++ ;
173558f2 1156
483b0559 1157 clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->IsInPreShower() ) );
173558f2 1158
483b0559 1159 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
173558f2 1160
483b0559 1161 iDigitInCluster++ ;
173558f2 1162
483b0559 1163 digitsC->Remove(digit) ;
173558f2 1164
483b0559 1165 nextdigit.Reset() ;
173558f2 1166
483b0559 1167
173558f2 1168
483b0559 1169 // Here we remove remaining Tower digits, which cannot make a cluster
173558f2 1170
483b0559 1171
173558f2 1172
483b0559 1173 if( notremoved ) {
173558f2 1174
483b0559 1175 while( ( digit = dynamic_cast<AliEMCALDigit *>(nextdigit()) ) ) {
173558f2 1176
483b0559 1177 if( IsInTower(digit) )
173558f2 1178
483b0559 1179 digitsC->Remove(digit) ;
173558f2 1180
483b0559 1181 else
173558f2 1182
483b0559 1183 break ;
173558f2 1184
483b0559 1185 }
173558f2 1186
483b0559 1187 notremoved = kFALSE ;
173558f2 1188
483b0559 1189 }
173558f2 1190
483b0559 1191
173558f2 1192
483b0559 1193 } // else
173558f2 1194
483b0559 1195
173558f2 1196
483b0559 1197 nextdigit.Reset() ;
173558f2 1198
483b0559 1199
173558f2 1200
483b0559 1201 AliEMCALDigit * digitN ;
173558f2 1202
483b0559 1203 index = 0 ;
173558f2 1204
483b0559 1205 while (index < iDigitInCluster){ // scan over digits already in cluster
173558f2 1206
483b0559 1207 digit = (AliEMCALDigit*)digits->At(clusterdigitslist[index]) ;
173558f2 1208
483b0559 1209 index++ ;
173558f2 1210
483b0559 1211 while ( (digitN = (AliEMCALDigit *)nextdigit()) ) { // scan over the reduced list of digits
173558f2 1212
483b0559 1213 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
173558f2 1214
483b0559 1215 switch (ineb ) {
173558f2 1216
483b0559 1217 case 0 : // not a neighbour
173558f2 1218
483b0559 1219 break ;
173558f2 1220
483b0559 1221 case 1 : // are neighbours
173558f2 1222
483b0559 1223 clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), digitN->IsInPreShower() ) ) ;
173558f2 1224
483b0559 1225 clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ;
173558f2 1226
483b0559 1227 iDigitInCluster++ ;
173558f2 1228
483b0559 1229 digitsC->Remove(digitN) ;
173558f2 1230
483b0559 1231 break ;
173558f2 1232
483b0559 1233 case 2 : // too far from each other
173558f2 1234
483b0559 1235 goto endofloop;
173558f2 1236
483b0559 1237 } // switch
173558f2 1238
483b0559 1239
173558f2 1240
483b0559 1241 } // while digitN
173558f2 1242
483b0559 1243
173558f2 1244
483b0559 1245 endofloop: ;
173558f2 1246
483b0559 1247 nextdigit.Reset() ;
173558f2 1248
483b0559 1249
173558f2 1250
483b0559 1251 } // loop over cluster
173558f2 1252
483b0559 1253
173558f2 1254
483b0559 1255 } // energy theshold
173558f2 1256
483b0559 1257
173558f2 1258
483b0559 1259
173558f2 1260
483b0559 1261 } // while digit
173558f2 1262
483b0559 1263
173558f2 1264
483b0559 1265 delete digitsC ;
173558f2 1266
483b0559 1267
173558f2 1268
483b0559 1269}
1270
173558f2 1271
1272
483b0559 1273//____________________________________________________________________________
173558f2 1274
483b0559 1275void AliEMCALClusterizerv1::MakeUnfolding()
173558f2 1276
483b0559 1277{
173558f2 1278
483b0559 1279 Fatal("AliEMCALClusterizerv1::MakeUnfolding", "--> Unfolding not implemented") ;
173558f2 1280
483b0559 1281
173558f2 1282
483b0559 1283// // Unfolds clusters using the shape of an ElectroMagnetic shower
173558f2 1284
483b0559 1285// // Performs unfolding of all EMC/CPV clusters
1286
173558f2 1287
1288
483b0559 1289// AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
173558f2 1290
483b0559 1291
173558f2 1292
483b0559 1293// const AliEMCALGeometry * geom = gime->EMCALGeometry() ;
173558f2 1294
483b0559 1295// TObjArray * emcRecPoints = gime->TowerRecPoints() ;
173558f2 1296
483b0559 1297// TObjArray * cpvRecPoints = gime->PreShoRecPoints() ;
173558f2 1298
483b0559 1299// TClonesArray * digits = gime->Digits() ;
173558f2 1300
483b0559 1301
173558f2 1302
483b0559 1303// // Unfold first EMC clusters
173558f2 1304
483b0559 1305// if(fNumberOfTowerClusters > 0){
1306
173558f2 1307
1308
483b0559 1309// Int_t nModulesToUnfold = geom->GetNModules() ;
1310
173558f2 1311
1312
483b0559 1313// Int_t numberofNotUnfolded = fNumberOfTowerClusters ;
173558f2 1314
483b0559 1315// Int_t index ;
173558f2 1316
483b0559 1317// for(index = 0 ; index < numberofNotUnfolded ; index++){
173558f2 1318
483b0559 1319
173558f2 1320
483b0559 1321// AliEMCALTowerRecPoint * emcRecPoint = (AliEMCALTowerRecPoint *) emcRecPoints->At(index) ;
173558f2 1322
483b0559 1323// if(emcRecPoint->GetEMCALMod()> nModulesToUnfold)
173558f2 1324
483b0559 1325// break ;
173558f2 1326
483b0559 1327
173558f2 1328
483b0559 1329// Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
173558f2 1330
483b0559 1331// Int_t * maxAt = new Int_t[nMultipl] ;
173558f2 1332
483b0559 1333// Float_t * maxAtEnergy = new Float_t[nMultipl] ;
173558f2 1334
483b0559 1335// Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fTowerLocMaxCut,digits) ;
173558f2 1336
483b0559 1337
173558f2 1338
483b0559 1339// if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
173558f2 1340
483b0559 1341// UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
173558f2 1342
483b0559 1343// emcRecPoints->Remove(emcRecPoint);
173558f2 1344
483b0559 1345// emcRecPoints->Compress() ;
173558f2 1346
483b0559 1347// index-- ;
173558f2 1348
483b0559 1349// fNumberOfTowerClusters -- ;
173558f2 1350
483b0559 1351// numberofNotUnfolded-- ;
173558f2 1352
483b0559 1353// }
173558f2 1354
483b0559 1355
173558f2 1356
483b0559 1357// delete[] maxAt ;
173558f2 1358
483b0559 1359// delete[] maxAtEnergy ;
173558f2 1360
483b0559 1361// }
173558f2 1362
483b0559 1363// }
173558f2 1364
483b0559 1365// // Unfolding of EMC clusters finished
1366
1367
173558f2 1368
1369
1370
483b0559 1371// // Unfold now CPV clusters
173558f2 1372
483b0559 1373// if(fNumberOfPreShoClusters > 0){
173558f2 1374
483b0559 1375
173558f2 1376
483b0559 1377// Int_t nModulesToUnfold = geom->GetNModules() ;
1378
173558f2 1379
1380
483b0559 1381// Int_t numberofPreShoNotUnfolded = fNumberOfPreShoClusters ;
173558f2 1382
483b0559 1383// Int_t index ;
173558f2 1384
483b0559 1385// for(index = 0 ; index < numberofPreShoNotUnfolded ; index++){
173558f2 1386
483b0559 1387
173558f2 1388
483b0559 1389// AliEMCALRecPoint * recPoint = (AliEMCALRecPoint *) cpvRecPoints->At(index) ;
1390
173558f2 1391
1392
483b0559 1393// if(recPoint->GetEMCALMod()> nModulesToUnfold)
173558f2 1394
483b0559 1395// break ;
173558f2 1396
483b0559 1397
173558f2 1398
483b0559 1399// AliEMCALTowerRecPoint * emcRecPoint = (AliEMCALTowerRecPoint*) recPoint ;
173558f2 1400
483b0559 1401
173558f2 1402
483b0559 1403// Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
173558f2 1404
483b0559 1405// Int_t * maxAt = new Int_t[nMultipl] ;
173558f2 1406
483b0559 1407// Float_t * maxAtEnergy = new Float_t[nMultipl] ;
173558f2 1408
483b0559 1409// Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fPreShoLocMaxCut,digits) ;
173558f2 1410
483b0559 1411
173558f2 1412
483b0559 1413// if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
173558f2 1414
483b0559 1415// UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
173558f2 1416
483b0559 1417// cpvRecPoints->Remove(emcRecPoint);
173558f2 1418
483b0559 1419// cpvRecPoints->Compress() ;
173558f2 1420
483b0559 1421// index-- ;
173558f2 1422
483b0559 1423// numberofPreShoNotUnfolded-- ;
173558f2 1424
483b0559 1425// fNumberOfPreShoClusters-- ;
173558f2 1426
483b0559 1427// }
173558f2 1428
483b0559 1429
173558f2 1430
483b0559 1431// delete[] maxAt ;
173558f2 1432
483b0559 1433// delete[] maxAtEnergy ;
173558f2 1434
483b0559 1435// }
173558f2 1436
483b0559 1437// }
173558f2 1438
483b0559 1439// //Unfolding of PreSho clusters finished
173558f2 1440
483b0559 1441
173558f2 1442
483b0559 1443}
1444
173558f2 1445
1446
483b0559 1447//____________________________________________________________________________
173558f2 1448
483b0559 1449Double_t AliEMCALClusterizerv1::ShowerShape(Double_t r)
173558f2 1450
483b0559 1451{
173558f2 1452
483b0559 1453 // Shape of the shower (see EMCAL TDR)
173558f2 1454
483b0559 1455 // If you change this function, change also the gradient evaluation in ChiSquare()
1456
173558f2 1457
1458
483b0559 1459 Double_t r4 = r*r*r*r ;
173558f2 1460
483b0559 1461 Double_t r295 = TMath::Power(r, 2.95) ;
173558f2 1462
483b0559 1463 Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
173558f2 1464
483b0559 1465 return shape ;
173558f2 1466
483b0559 1467}
1468
173558f2 1469
1470
483b0559 1471//____________________________________________________________________________
173558f2 1472
483b0559 1473void AliEMCALClusterizerv1::UnfoldCluster(AliEMCALTowerRecPoint * iniTower,
173558f2 1474
483b0559 1475 Int_t nMax,
173558f2 1476
a0636361 1477 AliEMCALDigit ** maxAt,
173558f2 1478
483b0559 1479 Float_t * maxAtEnergy)
173558f2 1480
483b0559 1481{
173558f2 1482
483b0559 1483 // Performs the unfolding of a cluster with nMax overlapping showers
173558f2 1484
483b0559 1485
173558f2 1486
483b0559 1487 Fatal("AliEMCALClusterizerv1::UnfoldCluster", "--> Unfolding not implemented") ;
1488
173558f2 1489
1490
483b0559 1491 // AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
173558f2 1492
483b0559 1493// const AliEMCALGeometry * geom = gime->EMCALGeometry() ;
173558f2 1494
483b0559 1495// const TClonesArray * digits = gime->Digits() ;
173558f2 1496
483b0559 1497// TObjArray * emcRecPoints = gime->TowerRecPoints() ;
173558f2 1498
483b0559 1499// TObjArray * cpvRecPoints = gime->PreShoRecPoints() ;
1500
173558f2 1501
1502
483b0559 1503// Int_t nPar = 3 * nMax ;
173558f2 1504
483b0559 1505// Float_t * fitparameters = new Float_t[nPar] ;
1506
173558f2 1507
1508
483b0559 1509// Bool_t rv = FindFit(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ;
173558f2 1510
483b0559 1511// if( !rv ) {
173558f2 1512
483b0559 1513// // Fit failed, return and remove cluster
173558f2 1514
483b0559 1515// delete[] fitparameters ;
173558f2 1516
483b0559 1517// return ;
173558f2 1518
483b0559 1519// }
1520
173558f2 1521
1522
483b0559 1523// // create ufolded rec points and fill them with new energy lists
173558f2 1524
483b0559 1525// // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
173558f2 1526
483b0559 1527// // and later correct this number in acordance with actual energy deposition
1528
173558f2 1529
1530
483b0559 1531// Int_t nDigits = iniTower->GetMultiplicity() ;
173558f2 1532
483b0559 1533// Float_t * efit = new Float_t[nDigits] ;
173558f2 1534
483b0559 1535// Float_t xDigit=0.,zDigit=0.,distance=0. ;
173558f2 1536
483b0559 1537// Float_t xpar=0.,zpar=0.,epar=0. ;
173558f2 1538
483b0559 1539// Int_t relid[4] ;
173558f2 1540
483b0559 1541// AliEMCALDigit * digit = 0 ;
173558f2 1542
483b0559 1543// Int_t * emcDigits = iniTower->GetDigitsList() ;
1544
173558f2 1545
1546
483b0559 1547// Int_t iparam ;
173558f2 1548
483b0559 1549// Int_t iDigit ;
173558f2 1550
483b0559 1551// for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
173558f2 1552
483b0559 1553// digit = (AliEMCALDigit*) digits->At(emcDigits[iDigit] ) ;
173558f2 1554
483b0559 1555// geom->AbsToRelNumbering(digit->GetId(), relid) ;
173558f2 1556
483b0559 1557// geom->RelPosInModule(relid, xDigit, zDigit) ;
173558f2 1558
483b0559 1559// efit[iDigit] = 0;
1560
173558f2 1561
1562
483b0559 1563// iparam = 0 ;
173558f2 1564
483b0559 1565// while(iparam < nPar ){
173558f2 1566
483b0559 1567// xpar = fitparameters[iparam] ;
173558f2 1568
483b0559 1569// zpar = fitparameters[iparam+1] ;
173558f2 1570
483b0559 1571// epar = fitparameters[iparam+2] ;
173558f2 1572
483b0559 1573// iparam += 3 ;
173558f2 1574
483b0559 1575// distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ;
173558f2 1576
483b0559 1577// distance = TMath::Sqrt(distance) ;
173558f2 1578
483b0559 1579// efit[iDigit] += epar * ShowerShape(distance) ;
173558f2 1580
483b0559 1581// }
173558f2 1582
483b0559 1583// }
173558f2 1584
483b0559 1585
1586
173558f2 1587
1588
483b0559 1589// // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
173558f2 1590
483b0559 1591// // so that energy deposited in each cell is distributed betwin new clusters proportionally
173558f2 1592
483b0559 1593// // to its contribution to efit
1594
173558f2 1595
1596
483b0559 1597// Float_t * emcEnergies = iniTower->GetEnergiesList() ;
173558f2 1598
483b0559 1599// Float_t ratio ;
1600
173558f2 1601
1602
483b0559 1603// iparam = 0 ;
173558f2 1604
483b0559 1605// while(iparam < nPar ){
173558f2 1606
483b0559 1607// xpar = fitparameters[iparam] ;
173558f2 1608
483b0559 1609// zpar = fitparameters[iparam+1] ;
173558f2 1610
483b0559 1611// epar = fitparameters[iparam+2] ;
173558f2 1612
483b0559 1613// iparam += 3 ;
173558f2 1614
483b0559 1615
173558f2 1616
483b0559 1617// AliEMCALTowerRecPoint * emcRP = 0 ;
1618
173558f2 1619
1620
483b0559 1621// if(iniTower->IsTower()){ //create new entries in fTowerRecPoints...
173558f2 1622
483b0559 1623
173558f2 1624
483b0559 1625// if(fNumberOfTowerClusters >= emcRecPoints->GetSize())
173558f2 1626
483b0559 1627// emcRecPoints->Expand(2*fNumberOfTowerClusters) ;
173558f2 1628
483b0559 1629
173558f2 1630
483b0559 1631// (*emcRecPoints)[fNumberOfTowerClusters] = new AliEMCALTowerRecPoint("") ;
173558f2 1632
483b0559 1633// emcRP = (AliEMCALTowerRecPoint *) emcRecPoints->At(fNumberOfTowerClusters);
173558f2 1634
483b0559 1635// fNumberOfTowerClusters++ ;
173558f2 1636
483b0559 1637// }
173558f2 1638
483b0559 1639// else{//create new entries in fPreShoRecPoints
173558f2 1640
483b0559 1641// if(fNumberOfPreShoClusters >= cpvRecPoints->GetSize())
173558f2 1642
483b0559 1643// cpvRecPoints->Expand(2*fNumberOfPreShoClusters) ;
173558f2 1644
483b0559 1645
173558f2 1646
483b0559 1647// (*cpvRecPoints)[fNumberOfPreShoClusters] = new AliEMCALPreShoRecPoint("") ;
173558f2 1648
483b0559 1649// emcRP = (AliEMCALTowerRecPoint *) cpvRecPoints->At(fNumberOfPreShoClusters);
173558f2 1650
483b0559 1651// fNumberOfPreShoClusters++ ;
173558f2 1652
483b0559 1653// }
173558f2 1654
483b0559 1655
173558f2 1656
483b0559 1657// Float_t eDigit ;
173558f2 1658
483b0559 1659// for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
173558f2 1660
483b0559 1661// digit = (AliEMCALDigit*) digits->At( emcDigits[iDigit] ) ;
173558f2 1662
483b0559 1663// geom->AbsToRelNumbering(digit->GetId(), relid) ;
173558f2 1664
483b0559 1665// geom->RelPosInModule(relid, xDigit, zDigit) ;
173558f2 1666
483b0559 1667// distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ;
173558f2 1668
483b0559 1669// distance = TMath::Sqrt(distance) ;
173558f2 1670
483b0559 1671// ratio = epar * ShowerShape(distance) / efit[iDigit] ;
173558f2 1672
483b0559 1673// eDigit = emcEnergies[iDigit] * ratio ;
173558f2 1674
483b0559 1675// emcRP->AddDigit( *digit, eDigit ) ;
173558f2 1676
483b0559 1677// }
173558f2 1678
483b0559 1679// }
173558f2 1680
483b0559 1681
173558f2 1682
483b0559 1683// delete[] fitparameters ;
173558f2 1684
483b0559 1685// delete[] efit ;
173558f2 1686
483b0559 1687
173558f2 1688
483b0559 1689}
1690
173558f2 1691
1692
483b0559 1693//_____________________________________________________________________________
173558f2 1694
483b0559 1695void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
173558f2 1696
483b0559 1697{
173558f2 1698
483b0559 1699 // Calculates the Chi square for the cluster unfolding minimization
173558f2 1700
483b0559 1701 // Number of parameters, Gradient, Chi squared, parameters, what to do
1702
173558f2 1703
1704
483b0559 1705 abort() ;
173558f2 1706
483b0559 1707 // Fatal("AliEMCALClusterizerv1::UnfoldingChiSquare","-->Unfolding not implemented") ;
1708
173558f2 1709
1710
483b0559 1711// TList * toMinuit = (TList*) gMinuit->GetObjectFit() ;
1712
173558f2 1713
1714
483b0559 1715// AliEMCALTowerRecPoint * emcRP = (AliEMCALTowerRecPoint*) toMinuit->At(0) ;
173558f2 1716
483b0559 1717// TClonesArray * digits = (TClonesArray*)toMinuit->At(1) ;
1718
1719
173558f2 1720
1721
1722
483b0559 1723
173558f2 1724
483b0559 1725// // AliEMCALTowerRecPoint * emcRP = (AliEMCALTowerRecPoint *) gMinuit->GetObjectFit() ; // TowerRecPoint to fit
1726
173558f2 1727
1728
483b0559 1729// Int_t * emcDigits = emcRP->GetDigitsList() ;
1730
173558f2 1731
1732
483b0559 1733// Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
1734
173558f2 1735
1736
483b0559 1737// Float_t * emcEnergies = emcRP->GetEnergiesList() ;
1738
173558f2 1739
1740
483b0559 1741// const AliEMCALGeometry * geom = AliEMCALGetter::GetInstance()->EMCALGeometry() ;
173558f2 1742
483b0559 1743// fret = 0. ;
173558f2 1744
483b0559 1745// Int_t iparam ;
1746
173558f2 1747
1748
483b0559 1749// if(iflag == 2)
173558f2 1750
483b0559 1751// for(iparam = 0 ; iparam < nPar ; iparam++)
173558f2 1752
483b0559 1753// Grad[iparam] = 0 ; // Will evaluate gradient
173558f2 1754
483b0559 1755
173558f2 1756
483b0559 1757// Double_t efit ;
1758
173558f2 1759
1760
483b0559 1761// AliEMCALDigit * digit ;
173558f2 1762
483b0559 1763// Int_t iDigit ;
1764
173558f2 1765
1766
483b0559 1767// for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
1768
173558f2 1769
1770
483b0559 1771// digit = (AliEMCALDigit*) digits->At( emcDigits[iDigit] ) ;
1772
173558f2 1773
1774
483b0559 1775// Int_t relid[4] ;
173558f2 1776
483b0559 1777// Float_t xDigit ;
173558f2 1778
483b0559 1779// Float_t zDigit ;
1780
173558f2 1781
1782
483b0559 1783// geom->AbsToRelNumbering(digit->GetId(), relid) ;
1784
173558f2 1785
1786
483b0559 1787// geom->RelPosInModule(relid, xDigit, zDigit) ;
1788
173558f2 1789
1790
483b0559 1791// if(iflag == 2){ // calculate gradient
173558f2 1792
483b0559 1793// Int_t iParam = 0 ;
173558f2 1794
483b0559 1795// efit = 0 ;
173558f2 1796
483b0559 1797// while(iParam < nPar ){
173558f2 1798
483b0559 1799// Double_t distance = (xDigit - x[iParam]) * (xDigit - x[iParam]) ;
173558f2 1800
483b0559 1801// iParam++ ;
173558f2 1802
483b0559 1803// distance += (zDigit - x[iParam]) * (zDigit - x[iParam]) ;
173558f2 1804
483b0559 1805// distance = TMath::Sqrt( distance ) ;
173558f2 1806
483b0559 1807// iParam++ ;
173558f2 1808
483b0559 1809// efit += x[iParam] * ShowerShape(distance) ;
173558f2 1810
483b0559 1811// iParam++ ;
173558f2 1812
483b0559 1813// }
173558f2 1814
483b0559 1815// Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
173558f2 1816
483b0559 1817// iParam = 0 ;
173558f2 1818
483b0559 1819// while(iParam < nPar ){
173558f2 1820
483b0559 1821// Double_t xpar = x[iParam] ;
173558f2 1822
483b0559 1823// Double_t zpar = x[iParam+1] ;
173558f2 1824
483b0559 1825// Double_t epar = x[iParam+2] ;
173558f2 1826
483b0559 1827// Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
173558f2 1828
483b0559 1829// Double_t shape = sum * ShowerShape(dr) ;
173558f2 1830
483b0559 1831// Double_t r4 = dr*dr*dr*dr ;
173558f2 1832
483b0559 1833// Double_t r295 = TMath::Power(dr,2.95) ;
173558f2 1834
483b0559 1835// Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
173558f2 1836
483b0559 1837// 0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
173558f2 1838
483b0559 1839
173558f2 1840
483b0559 1841// Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x
173558f2 1842
483b0559 1843// iParam++ ;
173558f2 1844
483b0559 1845// Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z
173558f2 1846
483b0559 1847// iParam++ ;
173558f2 1848
483b0559 1849// Grad[iParam] += shape ; // Derivative over energy
173558f2 1850
483b0559 1851// iParam++ ;
173558f2 1852
483b0559 1853// }
173558f2 1854
483b0559 1855// }
173558f2 1856
483b0559 1857// efit = 0;
173558f2 1858
483b0559 1859// iparam = 0 ;
1860
173558f2 1861
1862
483b0559 1863// while(iparam < nPar ){
173558f2 1864
483b0559 1865// Double_t xpar = x[iparam] ;
173558f2 1866
483b0559 1867// Double_t zpar = x[iparam+1] ;
173558f2 1868
483b0559 1869// Double_t epar = x[iparam+2] ;
173558f2 1870
483b0559 1871// iparam += 3 ;
173558f2 1872
483b0559 1873// Double_t distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ;
173558f2 1874
483b0559 1875// distance = TMath::Sqrt(distance) ;
173558f2 1876
483b0559 1877// efit += epar * ShowerShape(distance) ;
173558f2 1878
483b0559 1879// }
1880
173558f2 1881
1882
483b0559 1883// fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
173558f2 1884
483b0559 1885// // Here we assume, that sigma = sqrt(E)
173558f2 1886
483b0559 1887// }
1888
173558f2 1889
1890
483b0559 1891}
1892
173558f2 1893
1894
483b0559 1895//____________________________________________________________________________
173558f2 1896
483b0559 1897void AliEMCALClusterizerv1::Print(Option_t * option)const
173558f2 1898
483b0559 1899{
173558f2 1900
483b0559 1901 // Print clusterizer parameters
1902
173558f2 1903
1904
483b0559 1905 if( strcmp(GetName(), "") !=0 ){
173558f2 1906
483b0559 1907
173558f2 1908
483b0559 1909 // Print parameters
173558f2 1910
483b0559 1911
173558f2 1912
483b0559 1913 TString taskName(GetName()) ;
173558f2 1914
483b0559 1915 taskName.ReplaceAll(Version(), "") ;
1916
173558f2 1917
1918
483b0559 1919 cout << "---------------"<< taskName.Data() << " " << GetTitle()<< "-----------" << endl
173558f2 1920
05a92d59 1921 << "Clusterizing digits from the file: " << taskName.Data() << endl
173558f2 1922
05a92d59 1923 << " Branch: " << GetName() << endl
173558f2 1924
483b0559 1925 << endl
173558f2 1926
483b0559 1927 << " EMC Clustering threshold = " << fTowerClusteringThreshold << endl
173558f2 1928
483b0559 1929 << " EMC Local Maximum cut = " << fTowerLocMaxCut << endl
173558f2 1930
483b0559 1931 << " EMC Logarothmic weight = " << fW0 << endl
173558f2 1932
483b0559 1933 << endl
173558f2 1934
483b0559 1935 << " CPV Clustering threshold = " << fPreShoClusteringThreshold << endl
173558f2 1936
483b0559 1937 << " CPV Local Maximum cut = " << fPreShoLocMaxCut << endl
173558f2 1938
483b0559 1939 << " CPV Logarothmic weight = " << fW0CPV << endl
173558f2 1940
483b0559 1941 << endl ;
173558f2 1942
483b0559 1943 if(fToUnfold)
173558f2 1944
483b0559 1945 cout << " Unfolding on " << endl ;
173558f2 1946
483b0559 1947 else
173558f2 1948
483b0559 1949 cout << " Unfolding off " << endl ;
173558f2 1950
483b0559 1951
173558f2 1952
483b0559 1953 cout << "------------------------------------------------------------------" <<endl ;
173558f2 1954
483b0559 1955 }
173558f2 1956
483b0559 1957 else
173558f2 1958
483b0559 1959 cout << " AliEMCALClusterizerv1 not initialized " << endl ;
173558f2 1960
483b0559 1961}
173558f2 1962
483b0559 1963//____________________________________________________________________________
173558f2 1964
483b0559 1965void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option)
173558f2 1966
483b0559 1967{
173558f2 1968
483b0559 1969 // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
1970
173558f2 1971
1972
483b0559 1973 TObjArray * towerRecPoints = AliEMCALGetter::GetInstance()->TowerRecPoints() ;
173558f2 1974
483b0559 1975 TObjArray * preshoRecPoints = AliEMCALGetter::GetInstance()->PreShowerRecPoints() ;
1976
173558f2 1977
1978
483b0559 1979 cout << "AliEMCALClusterizerv1: : event "<<gAlice->GetEvNumber() << endl ;
173558f2 1980
483b0559 1981 cout << " Found "<< towerRecPoints->GetEntriesFast() << " TOWER Rec Points and "
173558f2 1982
483b0559 1983 << preshoRecPoints->GetEntriesFast() << " PRE SHOWER RecPoints" << endl ;
1984
173558f2 1985
1986
483b0559 1987 fRecPointsInRun += towerRecPoints->GetEntriesFast() ;
173558f2 1988
483b0559 1989 fRecPointsInRun += preshoRecPoints->GetEntriesFast() ;
1990
173558f2 1991
1992
483b0559 1993 if(strstr(option,"all")) {
1994
173558f2 1995
1996
483b0559 1997 cout << "Tower clusters " << endl ;
173558f2 1998
483b0559 1999 cout << " Index Ene(MeV) Multi Module phi r theta Lambda 1 Lambda 2 # of prim Primaries list " << endl;
173558f2 2000
483b0559 2001
173558f2 2002
483b0559 2003 Int_t index ;
173558f2 2004
483b0559 2005 for (index = 0 ; index < towerRecPoints->GetEntries() ; index++) {
173558f2 2006
483b0559 2007 AliEMCALTowerRecPoint * rp = dynamic_cast<AliEMCALTowerRecPoint * >(towerRecPoints->At(index)) ;
173558f2 2008
483b0559 2009 TVector3 globalpos;
173558f2 2010
483b0559 2011 rp->GetGlobalPosition(globalpos);
173558f2 2012
483b0559 2013 Float_t lambda[2];
173558f2 2014
483b0559 2015 rp->GetElipsAxis(lambda);
173558f2 2016
483b0559 2017 Int_t * primaries;
173558f2 2018
483b0559 2019 Int_t nprimaries;
173558f2 2020
483b0559 2021 primaries = rp->GetPrimaries(nprimaries);
2022
173558f2 2023
2024
483b0559 2025 cout << setw(4) << rp->GetIndexInList() << " "
173558f2 2026
483b0559 2027 << setw(7) << setprecision(3) << rp->GetEnergy() << " "
173558f2 2028
483b0559 2029 << setw(3) << rp->GetMultiplicity() << " "
173558f2 2030
483b0559 2031 << setw(1) << rp->GetEMCALArm() << " "
173558f2 2032
483b0559 2033 << setw(5) << setprecision(4) << globalpos.X() << " "
173558f2 2034
483b0559 2035 << setw(5) << setprecision(4) << globalpos.Y() << " "
173558f2 2036
483b0559 2037 << setw(5) << setprecision(4) << globalpos.Z() << " "
173558f2 2038
483b0559 2039 << setw(4) << setprecision(2) << lambda[0] << " "
173558f2 2040
483b0559 2041 << setw(4) << setprecision(2) << lambda[1] << " "
173558f2 2042
483b0559 2043 << setw(2) << nprimaries << " " ;
173558f2 2044
483b0559 2045
173558f2 2046
483b0559 2047 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++)
173558f2 2048
483b0559 2049 cout << setw(4) << primaries[iprimary] << " " ;
173558f2 2050
483b0559 2051 cout << endl ;
173558f2 2052
483b0559 2053 }
2054
173558f2 2055
2056
483b0559 2057 //Now plot Pre shower recPoints
2058
173558f2 2059
2060
483b0559 2061 cout << "-----------------------------------------------------------------------"<<endl ;
2062
173558f2 2063
2064
483b0559 2065 cout << "PreShower clusters " << endl ;
173558f2 2066
483b0559 2067 cout << " Index Ene(MeV) Multi Module phi r theta Lambda 1 Lambda 2 # of prim Primaries list " << endl;
173558f2 2068
483b0559 2069
173558f2 2070
483b0559 2071 for (index = 0 ; index < preshoRecPoints->GetEntries() ; index++) {
173558f2 2072
483b0559 2073 AliEMCALTowerRecPoint * rp = dynamic_cast<AliEMCALTowerRecPoint *>(preshoRecPoints->At(index)) ;
173558f2 2074
483b0559 2075 TVector3 globalpos;
173558f2 2076
483b0559 2077 rp->GetGlobalPosition(globalpos);
173558f2 2078
483b0559 2079 Float_t lambda[2];
173558f2 2080
483b0559 2081 rp->GetElipsAxis(lambda);
173558f2 2082
483b0559 2083 Int_t * primaries;
173558f2 2084
483b0559 2085 Int_t nprimaries;
173558f2 2086
483b0559 2087 primaries = rp->GetPrimaries(nprimaries);
2088
173558f2 2089
2090
483b0559 2091 cout << setw(4) << rp->GetIndexInList() << " "
173558f2 2092
483b0559 2093 << setw(7) << setprecision(3) << rp->GetEnergy() << " "
173558f2 2094
483b0559 2095 << setw(3) << rp->GetMultiplicity() << " "
173558f2 2096
483b0559 2097 << setw(1) << rp->GetEMCALArm() << " "
173558f2 2098
483b0559 2099 << setw(5) << setprecision(4) << globalpos.X() << " "
173558f2 2100
483b0559 2101 << setw(5) << setprecision(4) << globalpos.Y() << " "
173558f2 2102
483b0559 2103 << setw(5) << setprecision(4) << globalpos.Z() << " "
173558f2 2104
483b0559 2105 << setw(4) << setprecision(2) << lambda[0] << " "
173558f2 2106
483b0559 2107 << setw(4) << setprecision(2) << lambda[1] << " "
173558f2 2108
483b0559 2109 << setw(2) << nprimaries << " " ;
173558f2 2110
483b0559 2111
173558f2 2112
483b0559 2113 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++)
173558f2 2114
483b0559 2115 cout << setw(4) << primaries[iprimary] << " " ;
173558f2 2116
483b0559 2117 cout << endl ;
173558f2 2118
483b0559 2119 }
2120
173558f2 2121
2122
483b0559 2123 cout << "-----------------------------------------------------------------------"<<endl ;
173558f2 2124
483b0559 2125 }
173558f2 2126
483b0559 2127}
2128
173558f2 2129
2130