EMCAL/AliEMCALClusterizerv1.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17 /* $Log:
  18    1 October 2000. Yuri Kharlov:
  19      AreNeighbours()
  20      PPSD upper layer is considered if number of layers>1
  21    18 October 2000. Yuri Kharlov:
  22      AliEMCALClusterizerv1()
  23      CPV clusterizing parameters added
  24      MakeClusters()
  25      After first PPSD digit remove EMC digits only once
  26 */
  27 //*-- Author: Yves Schutz (SUBATECH)  & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
  28 //  August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction
  29 //                           of new  IO (à la PHOS)
  30 //////////////////////////////////////////////////////////////////////////////
  31 //  Clusterization class. Performs clusterization (collects neighbouring active cells) and
  32 //  unfolds the clusters having several local maxima.
  33 //  Results are stored in TreeR#, branches EMCALTowerRP (EMC recPoints),
  34 //  EMCALPreShoRP (CPV RecPoints) and AliEMCALClusterizer (Clusterizer with all
  35 //  parameters including input digits branch title, thresholds etc.)
  36 //  This TTask is normally called from Reconstructioner, but can as well be used in
  37 //  standalone mode.
  38 // Use Case:
  39 //  root [0] AliEMCALClusterizerv1 * cl = new AliEMCALClusterizerv1("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")
  45 //               //sets another title for Digitis (input) branch
  46 //  root [3] cl->SetRecPointsBranch("recp2")
  47 //               //sets another title four output branches
  48 //  root [4] cl->SetTowerLocalMaxCut(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
  55
  56 // --- ROOT system ---
  57
  58 #include "TROOT.h"
  59 #include "TFile.h"
  60 #include "TFolder.h"
  61 #include "TMath.h"
  62 #include "TMinuit.h"
  63 #include "TTree.h"
  64 #include "TSystem.h"
  65 #include "TBenchmark.h"
  66
  67 // --- Standard library ---
  68
  69
  70 // --- AliRoot header files ---
  71
  72 #include "AliEMCALClusterizerv1.h"
  73 #include "AliEMCALDigit.h"
  74 #include "AliEMCALDigitizer.h"
  75 #include "AliEMCALTowerRecPoint.h"
  76 #include "AliEMCAL.h"
  77 #include "AliEMCALGetter.h"
  78 #include "AliEMCALGeometry.h"
  79 #include "AliRun.h"
  80
  81 ClassImp(AliEMCALClusterizerv1)
  82
  83 //____________________________________________________________________________
  84   AliEMCALClusterizerv1::AliEMCALClusterizerv1() : AliEMCALClusterizer()
  85 {
  86   // default ctor (to be used mainly by Streamer)
  87
  88   InitParameters() ;
  89   fDefaultInit = kTRUE ;
  90 }
  91
  92 //____________________________________________________________________________
  93 AliEMCALClusterizerv1::AliEMCALClusterizerv1(const char* headerFile, const char* name, const Bool_t toSplit)
  94 :AliEMCALClusterizer(headerFile, name, toSplit)
  95 {
  96   // ctor with the indication of the file where header Tree and digits Tree are stored
  97
  98   InitParameters() ;
  99   Init() ;
 100   fDefaultInit = kFALSE ;
 101
 102 }
 103
 104 //____________________________________________________________________________
 105   AliEMCALClusterizerv1::~AliEMCALClusterizerv1()
 106 {
 107   // dtor
 108   fSplitFile = 0 ;
 109
 110 }
 111
 112 //____________________________________________________________________________
 113 const TString AliEMCALClusterizerv1::BranchName() const
 114 {
 115   TString branchName(GetName() ) ;
 116   branchName.Remove(branchName.Index(Version())-1) ;
 117   return branchName ;
 118 }
 119
 120 //____________________________________________________________________________
 121 Float_t  AliEMCALClusterizerv1::Calibrate(Int_t amp, Int_t where) const
 122 {
 123   //To be replased later by the method, reading individual parameters from the database
 124   // where = 0 == PRE ; where = 1 == ECAL ; where = 2 == HCAL
 125   if ( where == 0 ) // calibrate as PRE section
 126     return -fADCpedestalPRE + amp * fADCchannelPRE ;
 127   else if (where == 1) //calibrate as EC section
 128     return -fADCpedestalEC + amp * fADCchannelEC ;
 129   else if (where == 2) //calibrate as HC section
 130     return -fADCpedestalHC + amp * fADCchannelHC ;
 131   else
 132     Fatal("Calibrate", "Something went wrong!") ;
 133   return -9999999. ;
 134 }
 135
 136 //____________________________________________________________________________
 137 void AliEMCALClusterizerv1::Exec(Option_t * option)
 138 {
 139   // Steering method
 140
 141   if( strcmp(GetName(), "")== 0 )
 142     Init() ;
 143
 144   if(strstr(option,"tim"))
 145     gBenchmark->Start("EMCALClusterizer");
 146
 147   if(strstr(option,"print"))
 148     Print("") ;
 149
 150   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
 151   if(gime->BranchExists("RecPoints"))
 152     return ;
 153   Int_t nevents = gime->MaxEvent() ;
 154   Int_t ievent ;
 155
 156   for(ievent = 0; ievent < nevents; ievent++){
 157
 158     gime->Event(ievent,"D") ;
 159
 160     if(ievent == 0)
 161       GetCalibrationParameters() ;
 162
 163     fNumberOfPREClusters = fNumberOfECClusters = fNumberOfHCClusters = 0 ;
 164
 165     MakeClusters() ;
 166
 167     if(fToUnfold)
 168       MakeUnfolding() ;
 169
 170     WriteRecPoints(ievent) ;
 171
 172     if(strstr(option,"deb"))
 173       PrintRecPoints(option) ;
 174
 175     //increment the total number of digits per run
 176     fRecPointsInRun += gime->PRERecPoints()->GetEntriesFast() ;
 177     fRecPointsInRun += gime->ECALRecPoints()->GetEntriesFast() ;
 178     fRecPointsInRun += gime->HCALRecPoints()->GetEntriesFast() ;
 179  }
 180
 181   if(strstr(option,"tim")){
 182     gBenchmark->Stop("EMCALClusterizer");
 183     Info("Exec", "took %f seconds for Clusterizing %f seconds per event",
 184          gBenchmark->GetCpuTime("EMCALClusterizer"), gBenchmark->GetCpuTime("EMCALClusterizer")/nevents ) ;
 185   }
 186
 187 }
 188
 189 //____________________________________________________________________________
 190 Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALTowerRecPoint * emcRP, AliEMCALDigit ** maxAt, Float_t * maxAtEnergy,
 191                                     Int_t nPar, Float_t * fitparameters) const
 192 {
 193   // Calls TMinuit to fit the energy distribution of a cluster with several maxima
 194   // The initial values for fitting procedure are set equal to the positions of local maxima.
 195   // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
 196
 197   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
 198   TClonesArray * digits = gime->Digits() ;
 199
 200
 201   gMinuit->mncler();                     // Reset Minuit's list of paramters
 202   gMinuit->SetPrintLevel(-1) ;           // No Printout
 203   gMinuit->SetFCN(AliEMCALClusterizerv1::UnfoldingChiSquare) ;
 204                                          // To set the address of the minimization function
 205   TList * toMinuit = new TList();
 206   toMinuit->AddAt(emcRP,0) ;
 207   toMinuit->AddAt(digits,1) ;
 208
 209   gMinuit->SetObjectFit(toMinuit) ;         // To tranfer pointer to UnfoldingChiSquare
 210
 211   // filling initial values for fit parameters
 212   AliEMCALDigit * digit ;
 213
 214   Int_t ierflg  = 0;
 215   Int_t index   = 0 ;
 216   Int_t nDigits = (Int_t) nPar / 3 ;
 217
 218   Int_t iDigit ;
 219
 220   AliEMCALGeometry * geom = gime->EMCALGeometry() ;
 221
 222   for(iDigit = 0; iDigit < nDigits; iDigit++){
 223     digit = maxAt[iDigit];
 224
 225     Int_t relid[4] ;
 226     Float_t x = 0.;
 227     Float_t z = 0.;
 228     geom->AbsToRelNumbering(digit->GetId(), relid) ;
 229     geom->PosInAlice(relid, x, z) ;
 230
 231     Float_t energy = maxAtEnergy[iDigit] ;
 232
 233     gMinuit->mnparm(index, "x",  x, 0.1, 0, 0, ierflg) ;
 234     index++ ;
 235     if(ierflg != 0){
 236       Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : x = %f",  x ) ;
 237       return kFALSE;
 238     }
 239     gMinuit->mnparm(index, "z",  z, 0.1, 0, 0, ierflg) ;
 240     index++ ;
 241     if(ierflg != 0){
 242        Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : z = %f", z) ;
 243       return kFALSE;
 244     }
 245     gMinuit->mnparm(index, "Energy",  energy , 0.05*energy, 0., 4.*energy, ierflg) ;
 246     index++ ;
 247     if(ierflg != 0){
 248      Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : energy = %f", energy) ;
 249       return kFALSE;
 250     }
 251   }
 252
 253   Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
 254                       //  depends on it.
 255   Double_t p1 = 1.0 ;
 256   Double_t p2 = 0.0 ;
 257
 258   gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls
 259   gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient
 260   gMinuit->SetMaxIterations(5);
 261   gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings
 262
 263   gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize
 264
 265   if(ierflg == 4){  // Minimum not found
 266     Error("FindFit", "EMCAL Unfolding  Fit not converged, cluster abandoned " ) ;
 267     return kFALSE ;
 268   }
 269   for(index = 0; index < nPar; index++){
 270     Double_t err ;
 271     Double_t val ;
 272     gMinuit->GetParameter(index, val, err) ;    // Returns value and error of parameter index
 273     fitparameters[index] = val ;
 274    }
 275
 276   delete toMinuit ;
 277   return kTRUE;
 278
 279 }
 280
 281 //____________________________________________________________________________
 282 void AliEMCALClusterizerv1::GetCalibrationParameters()
 283 {
 284   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
 285   const AliEMCALDigitizer * dig = gime->Digitizer(BranchName()) ;
 286
 287   fADCchannelPRE  = dig->GetPREchannel() ;
 288   fADCpedestalPRE = dig->GetPREpedestal() ;
 289
 290   fADCchannelEC   = dig->GetECchannel() ;
 291   fADCpedestalEC  = dig->GetECpedestal();
 292
 293   fADCchannelHC   = dig->GetHCchannel() ;
 294   fADCpedestalHC  = dig->GetHCpedestal();
 295 }
 296
 297 //____________________________________________________________________________
 298 void AliEMCALClusterizerv1::Init()
 299 {
 300   // Make all memory allocations which can not be done in default constructor.
 301   // Attach the Clusterizer task to the list of EMCAL tasks
 302
 303   if ( strcmp(GetTitle(), "") == 0 )
 304     SetTitle("galice.root") ;
 305
 306   TString branchname = GetName() ;
 307   branchname.Remove(branchname.Index(Version())-1) ;
 308
 309   AliEMCALGetter * gime = AliEMCALGetter::GetInstance(GetTitle(), branchname.Data(), fToSplit ) ;
 310   if ( gime == 0 ) {
 311     Error("Init", "Could not obtain the Getter object !" ) ;
 312     return ;
 313   }
 314
 315   fSplitFile = 0 ;
 316   if(fToSplit){
 317     // construct the name of the file as /path/EMCAL.SDigits.root
 318     //First - extract full path if necessary
 319     TString fileName(GetTitle()) ;
 320     Ssiz_t islash = fileName.Last('/') ;
 321     if(islash<fileName.Length())
 322       fileName.Remove(islash+1,fileName.Length()) ;
 323     else
 324       fileName="" ;
 325     // Next - append the file name
 326     fileName+="EMCAL.RecData." ;
 327     if((strcmp(branchname.Data(),"Default")!=0)&&(strcmp(branchname.Data(),"")!=0)){
 328       fileName+=branchname ;
 329       fileName+="." ;
 330     }
 331     fileName+="root" ;
 332     // Finally - check if the file already opened or open the file
 333     fSplitFile = static_cast<TFile*>(gROOT->GetFile(fileName.Data()));
 334     if(!fSplitFile)
 335       fSplitFile =  TFile::Open(fileName.Data(),"update") ;
 336   }
 337
 338   const AliEMCALGeometry * geom = gime->EMCALGeometry() ;
 339   fNTowers = geom->GetNZ() *  geom->GetNPhi() ;
 340   if(!gMinuit)
 341     gMinuit = new TMinuit(100) ;
 342
 343   gime->PostClusterizer(this) ;
 344   gime->PostRecPoints(branchname ) ;
 345
 346 }
 347
 348 //____________________________________________________________________________
 349 void AliEMCALClusterizerv1::InitParameters()
 350 {
 351   fNumberOfPREClusters = fNumberOfECClusters = fNumberOfHCClusters = 0 ;
 352   fPREClusteringThreshold  = 0.0001; // must be adjusted according to the noise leve set by digitizer
 353   fECClusteringThreshold   = 0.0045;  // must be adjusted according to the noise leve set by digitizer
 354   fHCClusteringThreshold   = 0.001;  // must be adjusted according to the noise leve set by digitizer
 355   fPRELocMaxCut = 0.03 ;
 356   fECLocMaxCut  = 0.03 ;
 357   fHCLocMaxCut  = 0.03 ;
 358
 359   fPREW0    = 4.0 ;
 360   fECW0     = 4.5 ;
 361   fHCW0     = 4.5 ;
 362
 363   fTimeGate = 1.e-8 ;
 364
 365   fToUnfold = kFALSE ;
 366
 367   TString clusterizerName( GetName()) ;
 368   if (clusterizerName.IsNull() )
 369     clusterizerName = "Default" ;
 370   clusterizerName.Append(":") ;
 371   clusterizerName.Append(Version()) ;
 372   SetName(clusterizerName) ;
 373   fRecPointsInRun          = 0 ;
 374
 375 }
 376
 377 //____________________________________________________________________________
 378 Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2)const
 379 {
 380   // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
 381   //                                       = 1 are neighbour
 382   //                                       = 2 are not neighbour but do not continue searching
 383   // neighbours are defined as digits having at least a common vertex
 384   // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
 385   //                                      which is compared to a digit (d2)  not yet in a cluster
 386
 387    AliEMCALGeometry * geom = AliEMCALGetter::GetInstance()->EMCALGeometry() ;
 388
 389   Int_t rv = 0 ;
 390
 391   Int_t relid1[4] ;
 392   geom->AbsToRelNumbering(d1->GetId(), relid1) ;
 393
 394   Int_t relid2[4] ;
 395   geom->AbsToRelNumbering(d2->GetId(), relid2) ;
 396
 397   if ( (relid1[0] == relid2[0]) && // inside the same EMCAL Arm
 398        (relid1[1]==relid2[1]) ) {  // and same tower section
 399     Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
 400     Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
 401
 402     if (( coldiff <= 1 )  && ( rowdiff <= 1 )){
 403       if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fTimeGate))
 404       rv = 1 ;
 405     }
 406     else {
 407       if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
 408         rv = 2; //  Difference in row numbers is too large to look further
 409     }
 410
 411   }
 412   else {
 413
 414     if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) )
 415       rv=0 ;
 416   }
 417
 418   if (gDebug == 2 )
 419     Info("AreNeighbours", "neighbours=%d, id1=%d, relid1=%d,%d,%d,%d \n id2=%d, relid2=%d,%d,%d,%d",
 420          rv, d1->GetId(), relid1[0], relid1[1], relid1[2], relid1[3],
 421          d2->GetId(), relid2[0], relid2[1], relid2[2], relid2[3]) ;
 422
 423   return rv ;
 424 }
 425
 426
 427 //____________________________________________________________________________
 428 void AliEMCALClusterizerv1::WriteRecPoints(Int_t event)
 429 {
 430
 431   // Creates new branches with given title
 432   // fills and writes into TreeR.
 433
 434   AliEMCALGetter *gime = AliEMCALGetter::GetInstance() ;
 435
 436   TObjArray * aPRERecPoints = gime->PRERecPoints() ;
 437   TObjArray * aECRecPoints = gime->ECALRecPoints() ;
 438   TObjArray * aHCRecPoints = gime->HCALRecPoints() ;
 439
 440   TClonesArray * digits = gime->Digits() ;
 441   TTree * treeR ;
 442
 443   if(fToSplit){
 444     if(!fSplitFile)
 445       return ;
 446     fSplitFile->cd() ;
 447     TString name("TreeR") ;
 448     name += event ;
 449     treeR = dynamic_cast<TTree*>(fSplitFile->Get(name));
 450   }
 451   else{
 452     treeR = gAlice->TreeR();
 453   }
 454
 455   if(!treeR){
 456     gAlice->MakeTree("R", fSplitFile);
 457     treeR = gAlice->TreeR() ;
 458   }
 459
 460   Int_t index ;
 461
 462   //Evaluate position, dispersion and other RecPoint properties for PRE section
 463   for(index = 0; index < aPRERecPoints->GetEntries(); index++)
 464     (dynamic_cast<AliEMCALRecPoint *>(aPRERecPoints->At(index)))->EvalAll(fPREW0,digits)  ;
 465   aPRERecPoints->Sort() ;
 466
 467   for(index = 0; index < aPRERecPoints->GetEntries(); index++)
 468     (dynamic_cast<AliEMCALRecPoint *>(aPRERecPoints->At(index)))->SetIndexInList(index) ;
 469
 470   aPRERecPoints->Expand(aPRERecPoints->GetEntriesFast()) ;
 471
 472   //Evaluate position, dispersion and other RecPoint properties for EC section
 473   for(index = 0; index < aECRecPoints->GetEntries(); index++)
 474     (dynamic_cast<AliEMCALTowerRecPoint *>(aECRecPoints->At(index)))->EvalAll(fECW0,digits) ;
 475
 476   aECRecPoints->Sort() ;
 477
 478   for(index = 0; index < aECRecPoints->GetEntries(); index++)
 479     (dynamic_cast<AliEMCALTowerRecPoint *>(aECRecPoints->At(index)))->SetIndexInList(index) ;
 480
 481   aECRecPoints->Expand(aECRecPoints->GetEntriesFast()) ;
 482
 483   //Evaluate position, dispersion and other RecPoint properties for HC section
 484   for(index = 0; index < aHCRecPoints->GetEntries(); index++)
 485     (dynamic_cast<AliEMCALTowerRecPoint *>(aHCRecPoints->At(index)))->EvalAll(fHCW0,digits) ;
 486
 487   aHCRecPoints->Sort() ;
 488
 489   for(index = 0; index < aHCRecPoints->GetEntries(); index++)
 490     (dynamic_cast<AliEMCALTowerRecPoint *>(aHCRecPoints->At(index)))->SetIndexInList(index) ;
 491
 492   aHCRecPoints->Expand(aHCRecPoints->GetEntriesFast()) ;
 493
 494   Int_t bufferSize = 32000 ;
 495   Int_t splitlevel = 0 ;
 496
 497   //PRE section branch
 498   TBranch * branchPRE = treeR->Branch("EMCALPRERP","TObjArray",&aPRERecPoints,bufferSize,splitlevel);
 499   branchPRE->SetTitle(BranchName());
 500
 501   //EC section branch
 502   TBranch * branchEC = treeR->Branch("EMCALECRP","TObjArray",&aECRecPoints,bufferSize,splitlevel);
 503   branchEC->SetTitle(BranchName());
 504
 505   //HC section branch
 506   TBranch * branchHC = treeR->Branch("EMCALHCRP","TObjArray",&aHCRecPoints,bufferSize,splitlevel);
 507   branchHC->SetTitle(BranchName());
 508
 509   //And Finally  clusterizer branch
 510   AliEMCALClusterizerv1 * cl = (AliEMCALClusterizerv1*)gime->Clusterizer(BranchName()) ;
 511   TBranch * clusterizerBranch = treeR->Branch("AliEMCALClusterizer","AliEMCALClusterizerv1",
 512                                               &cl,bufferSize,splitlevel);
 513   clusterizerBranch->SetTitle(BranchName());
 514
 515   branchPRE->Fill() ;
 516   branchEC->Fill() ;
 517   branchHC->Fill() ;
 518   clusterizerBranch->Fill() ;
 519
 520   treeR->AutoSave() ; //Write(0,kOverwrite) ;
 521   if(gAlice->TreeR()!=treeR)
 522     treeR->Delete();
 523 }
 524
 525 //____________________________________________________________________________
 526 void AliEMCALClusterizerv1::MakeClusters()
 527 {
 528   // Steering method to construct the clusters stored in a list of Reconstructed Points
 529   // A cluster is defined as a list of neighbour digits
 530
 531   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
 532
 533   AliEMCALGeometry * geom = gime->EMCALGeometry() ;
 534
 535
 536   TObjArray * aPRERecPoints = gime->PRERecPoints(BranchName()) ;
 537   TObjArray * aECRecPoints  = gime->ECALRecPoints(BranchName()) ;
 538   TObjArray * aHCRecPoints  = gime->HCALRecPoints(BranchName()) ;
 539
 540   aPRERecPoints->Delete() ;
 541   aECRecPoints->Delete() ;
 542   aHCRecPoints->Delete() ;
 543
 544   TClonesArray * digits = gime->Digits() ;
 545   if ( !digits ) {
 546     Fatal("MakeClusters -> Digits with name %s not found", GetName() ) ;
 547   }
 548
 549   TIter next(digits) ;
 550   AliEMCALDigit * digit ;
 551   Int_t ndigEC=0, ndigPRE=0, ndigHC=0 ;
 552
 553   // count the number of digits in EC section
 554   while ( (digit = dynamic_cast<AliEMCALDigit *>(next())) ) { // scan over the list of digits
 555     if (geom->IsInECAL(digit->GetId()))
 556       ndigEC++ ;
 557     else if (geom->IsInPRE(digit->GetId()))
 558       ndigPRE++;
 559     else if (geom->IsInHCAL(digit->GetId()))
 560       ndigHC++;
 561     else {
 562       Error("MakeClusters", "id = %d is a wrong ID!", digit->GetId()) ;
 563       abort() ;
 564     }
 565   }
 566
 567   // add amplitude of PRE and EC sections
 568   Int_t digEC ;
 569   for (digEC = 0 ; digEC < ndigEC ; digEC++) {
 570     digit = dynamic_cast<AliEMCALDigit *>(digits->At(digEC)) ;
 571     Int_t digPRE ;
 572     for (digPRE = ndigEC ; digPRE < ndigEC+ndigPRE ; digPRE++) {
 573       AliEMCALDigit *  digitPRE = dynamic_cast<AliEMCALDigit *>(digits->At(digPRE)) ;
 574       if ( geom->AreInSameTower(digit->GetId(), digitPRE->GetId()) ){
 575         Float_t  amp = static_cast<Float_t>(digit->GetAmp()) + geom->GetSummationFraction() * static_cast<Float_t>(digitPRE->GetAmp()) + 0.5 ;
 576         digit->SetAmp(static_cast<Int_t>(amp)) ;
 577         break ;
 578       }
 579     }
 580     if (gDebug)
 581       Info("MakeClusters","id = %d amp = %d", digit->GetId(), digit->GetAmp()) ;
 582   }
 583
 584   TClonesArray * digitsC =  dynamic_cast<TClonesArray*>(digits->Clone()) ;
 585
 586
 587   // Clusterization starts
 588
 589   TIter nextdigit(digitsC) ;
 590   Bool_t notremovedEC = kTRUE, notremovedPRE = kTRUE ;
 591
 592   while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // scan over the list of digitsC
 593     AliEMCALRecPoint * clu = 0 ;
 594
 595     TArrayI clusterPREdigitslist(50), clusterECdigitslist(50), clusterHCdigitslist(50);
 596
 597     Bool_t inPRE = kFALSE, inECAL = kFALSE, inHCAL = kFALSE ;
 598     if( geom->IsInPRE(digit->GetId()) ) {
 599       inPRE = kTRUE ;
 600     }
 601     else if( geom->IsInECAL(digit->GetId()) ) {
 602       inECAL = kTRUE ;
 603     }
 604     else if( geom->IsInHCAL(digit->GetId()) ) {
 605       inHCAL = kTRUE ;
 606     }
 607
 608     if (gDebug == 2) {
 609       if (inPRE)
 610         Info("MakeClusters","id = %d, ene = %f , thre = %f ",
 611              digit->GetId(),Calibrate(digit->GetAmp(), 0), fPREClusteringThreshold) ;
 612       if (inECAL)
 613         Info("MakeClusters","id = %d, ene = %f , thre = %f",
 614              digit->GetId(),Calibrate(digit->GetAmp(), 1), fECClusteringThreshold) ;
 615       if (inHCAL)
 616         Info("MakeClusters","id = %d, ene = %f , thre = %f",
 617              digit->GetId(),Calibrate(digit->GetAmp(), 2), fHCClusteringThreshold ) ;
 618     }
 619
 620     if ( (inPRE  && (Calibrate(digit->GetAmp(), 0) > fPREClusteringThreshold  )) ||
 621          (inECAL && (Calibrate(digit->GetAmp(), 1) > fECClusteringThreshold  ))  ||
 622          (inHCAL && (Calibrate(digit->GetAmp(), 2) > fHCClusteringThreshold  )) ) {
 623
 624       Int_t  iDigitInPRECluster = 0, iDigitInECCluster = 0, iDigitInHCCluster = 0;
 625       Int_t where ; // PRE = 0, ECAl = 1, HCAL = 2
 626
 627       // Find the seed in each of the section ECAL/PRE/HCAL
 628
 629       if( geom->IsInECAL(digit->GetId()) ) {
 630         where = 1 ; // to tell we are in ECAL
 631         // start a new Tower RecPoint
 632         if(fNumberOfECClusters >= aECRecPoints->GetSize())
 633           aECRecPoints->Expand(2*fNumberOfECClusters+1) ;
 634         AliEMCALTowerRecPoint * rp = new  AliEMCALTowerRecPoint("") ;
 635         rp->SetECAL() ;
 636         aECRecPoints->AddAt(rp, fNumberOfECClusters) ;
 637         clu = dynamic_cast<AliEMCALTowerRecPoint *>(aECRecPoints->At(fNumberOfECClusters)) ;
 638         fNumberOfECClusters++ ;
 639         clu->AddDigit(*digit, Calibrate(digit->GetAmp(), where)) ;
 640         clusterECdigitslist[iDigitInECCluster] = digit->GetIndexInList() ;
 641         iDigitInECCluster++ ;
 642         digitsC->Remove(digit) ;
 643         if (gDebug == 2 )
 644           Info("MakeClusters","OK id = %d, ene = %f , thre = %f ", digit->GetId(),Calibrate(digit->GetAmp(), 1), fECClusteringThreshold) ;
 645
 646       }
 647       else if( geom->IsInPRE(digit->GetId()) ) {
 648         where = 0 ; // to tell we are in PRE
 649         // start a new Pre Shower cluster
 650         if(fNumberOfPREClusters >= aPRERecPoints->GetSize())
 651           aPRERecPoints->Expand(2*fNumberOfPREClusters+1);
 652         AliEMCALTowerRecPoint * rp = new AliEMCALTowerRecPoint("") ;
 653         rp->SetPRE() ;
 654         aPRERecPoints->AddAt(rp, fNumberOfPREClusters) ;
 655         clu =  dynamic_cast<AliEMCALTowerRecPoint *>(aPRERecPoints->At(fNumberOfPREClusters))  ;
 656         fNumberOfPREClusters++ ;
 657         clu->AddDigit(*digit, Calibrate(digit->GetAmp(), where));
 658         clusterPREdigitslist[iDigitInPRECluster] = digit->GetIndexInList()  ;
 659         iDigitInPRECluster++ ;
 660         digitsC->Remove(digit) ;
 661         if (gDebug == 2 )
 662           Info("MakeClusters","OK id = %d, ene = %f , thre = %f ", digit->GetId(),Calibrate(digit->GetAmp(), 0), fPREClusteringThreshold) ;
 663
 664         nextdigit.Reset() ;
 665
 666         // Here we remove remaining EC digits, which cannot make a cluster
 667
 668         if( notremovedEC ) {
 669           while( ( digit = dynamic_cast<AliEMCALDigit *>(nextdigit()) ) ) {
 670             if( geom->IsInECAL(digit->GetId()) )
 671               digitsC->Remove(digit) ;
 672             else
 673               break ;
 674           }
 675           notremovedEC = kFALSE ;
 676         }
 677
 678       }
 679       else if( geom->IsInHCAL(digit->GetId()) ) {
 680         where = 2 ; // to tell we are in HCAL
 681         // start a new HCAL cluster
 682         if(fNumberOfHCClusters >= aHCRecPoints->GetSize())
 683           aHCRecPoints->Expand(2*fNumberOfHCClusters+1);
 684         AliEMCALTowerRecPoint * rp = new AliEMCALTowerRecPoint("") ;
 685         rp->SetHCAL() ;
 686         aHCRecPoints->AddAt(rp, fNumberOfHCClusters) ;
 687         clu =  dynamic_cast<AliEMCALTowerRecPoint *>(aHCRecPoints->At(fNumberOfHCClusters))  ;
 688         fNumberOfHCClusters++ ;
 689         clu->AddDigit(*digit, Calibrate(digit->GetAmp(), where));
 690         clusterHCdigitslist[iDigitInHCCluster] = digit->GetIndexInList()  ;
 691         iDigitInHCCluster++ ;
 692         digitsC->Remove(digit) ;
 693         if (gDebug == 2 )
 694           Info("MakeClusters","OK id = %d, ene = %f , thre = %f ", digit->GetId(),Calibrate(digit->GetAmp(), 2), fHCClusteringThreshold) ;
 695
 696         nextdigit.Reset() ;
 697
 698         // Here we remove remaining PRE digits, which cannot make a cluster
 699
 700         if( notremovedPRE ) {
 701           while( ( digit = dynamic_cast<AliEMCALDigit *>(nextdigit()) ) ) {
 702             if( geom->IsInPRE(digit->GetId()) )
 703               digitsC->Remove(digit) ;
 704             else
 705               break ;
 706           }
 707           notremovedPRE = kFALSE ;
 708         }
 709       }
 710
 711       nextdigit.Reset() ;
 712
 713       AliEMCALDigit * digitN ;
 714       Int_t index = 0 ;
 715
 716       // Do the Clustering in each of the three section ECAL/PRE/HCAL
 717
 718       while (index < iDigitInECCluster){ // scan over digits already in cluster
 719         digit =  (AliEMCALDigit*)digits->At(clusterECdigitslist[index])  ;
 720         index++ ;
 721         while ( (digitN = (AliEMCALDigit *)nextdigit()) ) { // scan over the reduced list of digits
 722           Int_t ineb = AreNeighbours(digit, digitN);       // call (digit,digitN) in THAT oder !!!!!
 723           // Info("MakeClusters","id1 = %d, id2 = %d , neighbours = %d", digit->GetId(), digitN->GetId(), ineb) ;
 724          switch (ineb ) {
 725           case 0 :   // not a neighbour
 726             break ;
 727           case 1 :   // are neighbours
 728             clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), 1) ) ;
 729             clusterECdigitslist[iDigitInECCluster] = digitN->GetIndexInList() ;
 730             iDigitInECCluster++ ;
 731             digitsC->Remove(digitN) ;
 732             break ;
 733           case 2 :   // too far from each other
 734             goto endofloop1;
 735           } // switch
 736
 737         } // while digitN
 738
 739       endofloop1: ;
 740         nextdigit.Reset() ;
 741       } // loop over EC cluster
 742
 743       index = 0 ;
 744       while (index < iDigitInPRECluster){ // scan over digits already in cluster
 745         digit =  (AliEMCALDigit*)digits->At(clusterPREdigitslist[index])  ;
 746         index++ ;
 747         while ( (digitN = (AliEMCALDigit *)nextdigit()) ) { // scan over the reduced list of digits
 748           Int_t ineb = AreNeighbours(digit, digitN);       // call (digit,digitN) in THAT oder !!!!!
 749           //      Info("MakeClusters","id1 = %d, id2 = %d , neighbours = %d", digit->GetId(), digitN->GetId(), ineb) ;
 750           switch (ineb ) {
 751           case 0 :   // not a neighbour
 752             break ;
 753           case 1 :   // are neighbours
 754             clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), 0) ) ;
 755             clusterPREdigitslist[iDigitInPRECluster] = digitN->GetIndexInList() ;
 756             iDigitInPRECluster++ ;
 757             digitsC->Remove(digitN) ;
 758             break ;
 759           case 2 :   // too far from each other
 760             goto endofloop2;
 761           } // switch
 762
 763         } // while digitN
 764
 765       endofloop2: ;
 766         nextdigit.Reset() ;
 767       } // loop over PRE cluster
 768
 769       index = 0 ;
 770       while (index < iDigitInHCCluster){ // scan over digits already in cluster
 771         digit =  (AliEMCALDigit*)digits->At(clusterHCdigitslist[index])  ;
 772         index++ ;
 773         while ( (digitN = (AliEMCALDigit *)nextdigit()) ) { // scan over the reduced list of digits
 774           Int_t ineb = AreNeighbours(digit, digitN);       // call (digit,digitN) in THAT oder !!!!!
 775           //Info("MakeClusters","id1 = %d, id2 = %d , neighbours = %d", digit->GetId(), digitN->GetId(), ineb) ;
 776           switch (ineb ) {
 777           case 0 :   // not a neighbour
 778             break ;
 779           case 1 :   // are neighbours
 780             clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), 2) ) ;
 781             clusterHCdigitslist[iDigitInHCCluster] = digitN->GetIndexInList() ;
 782             iDigitInHCCluster++ ;
 783             digitsC->Remove(digitN) ;
 784             break ;
 785           case 2 :   // too far from each other
 786             goto endofloop3;
 787           } // switch
 788         } // while digitN
 789
 790       endofloop3: ;
 791         nextdigit.Reset() ;
 792       } // loop over HC cluster
 793
 794     } // energy theshold
 795   } // while digit
 796   delete digitsC ;
 797 }
 798
 799 //____________________________________________________________________________
 800 void AliEMCALClusterizerv1::MakeUnfolding()
 801 {
 802   Fatal("AliEMCALClusterizerv1::MakeUnfolding", "--> Unfolding not implemented") ;
 803
 804 }
 805
 806 //____________________________________________________________________________
 807 Double_t  AliEMCALClusterizerv1::ShowerShape(Double_t r)
 808 {
 809   // Shape of the shower (see EMCAL TDR)
 810   // If you change this function, change also the gradient evaluation in ChiSquare()
 811
 812   Double_t r4    = r*r*r*r ;
 813   Double_t r295  = TMath::Power(r, 2.95) ;
 814   Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
 815   return shape ;
 816 }
 817
 818 //____________________________________________________________________________
 819 void  AliEMCALClusterizerv1::UnfoldCluster(AliEMCALTowerRecPoint * iniTower,
 820                                                  Int_t nMax,
 821                                                  AliEMCALDigit ** maxAt,
 822                                                  Float_t * maxAtEnergy)
 823 {
 824   // Performs the unfolding of a cluster with nMax overlapping showers
 825
 826   Fatal("UnfoldCluster", "--> Unfolding not implemented") ;
 827
 828 }
 829
 830 //_____________________________________________________________________________
 831 void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
 832 {
 833   // Calculates the Chi square for the cluster unfolding minimization
 834   // Number of parameters, Gradient, Chi squared, parameters, what to do
 835
 836   ::Fatal("UnfoldingChiSquare","Unfolding not implemented") ;
 837 }
 838 //____________________________________________________________________________
 839 void AliEMCALClusterizerv1::Print(Option_t * option)const
 840 {
 841   // Print clusterizer parameters
 842
 843   TString message("\n") ;
 844
 845   if( strcmp(GetName(), "") !=0 ){
 846
 847     // Print parameters
 848
 849     TString taskName(GetName()) ;
 850     taskName.ReplaceAll(Version(), "") ;
 851
 852     message += "--------------- " ;
 853     message += taskName.Data() ;
 854     message += " " ;
 855     message += GetTitle() ;
 856     message += "-----------\n" ;
 857     message += "Clusterizing digits from the file: " ;
 858     message += taskName.Data() ;
 859     message += "\n                           Branch: " ;
 860     message += GetName() ;
 861     message += "\n                       Pre Shower Clustering threshold = " ;
 862     message += fPREClusteringThreshold ;
 863     message += "\n                       Pre Shower  Local Maximum cut    = " ;
 864     message += fPRELocMaxCut ;
 865     message += "\n                       Pre Shower Logarothmic weight   = " ;
 866     message += fPREW0 ;
 867     message += "\n                       EC Clustering threshold = " ;
 868     message += fECClusteringThreshold ;
 869     message += "\n                       EC Local Maximum cut    = " ;
 870     message += fECLocMaxCut ;
 871     message += "\n                       EC Logarothmic weight   = " ;
 872     message += fECW0 ;
 873     message += "\n                       Pre Shower Clustering threshold = " ;
 874     message += fHCClusteringThreshold ;
 875     message += "\n                       HC Local Maximum cut    = " ;
 876     message += fHCLocMaxCut ;
 877     message += "\n                       HC Logarothmic weight   = " ;
 878     message += fHCW0 ;
 879     if(fToUnfold)
 880       message +="\nUnfolding on\n" ;
 881     else
 882       message += "\nUnfolding off\n";
 883
 884     message += "------------------------------------------------------------------" ;
 885   }
 886   else
 887     message += "AliEMCALClusterizerv1 not initialized " ;
 888
 889   Info("Print", message.Data() ) ;
 890 }
 891
 892 //____________________________________________________________________________
 893 void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option)
 894 {
 895   // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
 896
 897   TObjArray * aPRERecPoints = AliEMCALGetter::GetInstance()->PRERecPoints() ;
 898   TObjArray * aECRecPoints = AliEMCALGetter::GetInstance()->ECALRecPoints() ;
 899   TObjArray * aHCRecPoints = AliEMCALGetter::GetInstance()->HCALRecPoints() ;
 900
 901   Info("PrintRecPoints", "Clusterization result:") ;
 902
 903   printf("event # %d\n", gAlice->GetEvNumber() ) ;
 904   printf("           Found %d PRE SHOWER RecPoints, %d EC Rec Points and %d HC Rec Points\n ",
 905          aPRERecPoints->GetEntriesFast(), aECRecPoints->GetEntriesFast(), aHCRecPoints->GetEntriesFast() ) ;
 906
 907   fRecPointsInRun +=  aPRERecPoints->GetEntriesFast() ;
 908   fRecPointsInRun +=  aECRecPoints->GetEntriesFast() ;
 909   fRecPointsInRun +=  aHCRecPoints->GetEntriesFast() ;
 910
 911   if(strstr(option,"all")) {
 912
 913     //Pre shower recPoints
 914
 915     printf("-----------------------------------------------------------------------\n") ;
 916     printf("Clusters in PRE section\n") ;
 917     printf("Index    Ene(GeV) Multi Module     phi     r   theta    X    Y      Z   Dispersion Lambda 1   Lambda 2  # of prim  Primaries list\n") ;
 918
 919     Int_t index ;
 920
 921     for (index = 0 ; index < aPRERecPoints->GetEntries() ; index++) {
 922       AliEMCALTowerRecPoint * rp = dynamic_cast<AliEMCALTowerRecPoint *>(aPRERecPoints->At(index)) ;
 923       TVector3  globalpos;
 924       rp->GetGlobalPosition(globalpos);
 925       TVector3  localpos;
 926       rp->GetLocalPosition(localpos);
 927       Float_t lambda[2];
 928       rp->GetElipsAxis(lambda);
 929       Int_t * primaries;
 930       Int_t nprimaries;
 931       primaries = rp->GetPrimaries(nprimaries);
 932       printf("\n%6d  %8.4f  %3d     %2d     %4.1f    %4.1f %4.1f  %4.1f %4.1f %4.1f    %4.1f   %4f  %4f    %2d     : ",
 933              rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetEMCALArm(),
 934              globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(),
 935              rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ;
 936       for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
 937         printf("%d ", primaries[iprimary] ) ;
 938       }
 939     }
 940
 941     printf("\n-----------------------------------------------------------------------\n") ;
 942     printf("Clusters in ECAL section\n") ;
 943     printf("Index    Ene(GeV) Multi Module     phi     r   theta    X    Y      Z   Dispersion Lambda 1   Lambda 2  # of prim  Primaries list\n") ;
 944
 945     for (index = 0 ; index < aECRecPoints->GetEntries() ; index++) {
 946       AliEMCALTowerRecPoint * rp = dynamic_cast<AliEMCALTowerRecPoint * >(aECRecPoints->At(index)) ;
 947       TVector3  globalpos;
 948       rp->GetGlobalPosition(globalpos);
 949       TVector3  localpos;
 950       rp->GetLocalPosition(localpos);
 951       Float_t lambda[2];
 952       rp->GetElipsAxis(lambda);
 953       Int_t * primaries;
 954       Int_t nprimaries;
 955       primaries = rp->GetPrimaries(nprimaries);
 956       printf("\n%6d  %8.4f  %3d     %2d     %4.1f    %4.1f %4.1f  %4.1f %4.1f %4.1f    %4.1f   %4f  %4f    %2d     : ",
 957              rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetEMCALArm(),
 958              globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(),
 959              rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ;
 960       for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
 961         printf("%d ", primaries[iprimary] ) ;
 962       }
 963     }
 964
 965     printf("\n-----------------------------------------------------------------------\n") ;
 966     printf("Clusters in HCAL section\n") ;
 967     printf("Index    Ene(GeV) Multi Module     phi     r   theta    X    Y      Z   Dispersion Lambda 1   Lambda 2  # of prim  Primaries list\n") ;
 968
 969     for (index = 0 ; index < aHCRecPoints->GetEntries() ; index++) {
 970       AliEMCALTowerRecPoint * rp = dynamic_cast<AliEMCALTowerRecPoint * >(aHCRecPoints->At(index)) ;
 971       TVector3  globalpos;
 972       rp->GetGlobalPosition(globalpos);
 973       TVector3  localpos;
 974       rp->GetLocalPosition(localpos);
 975       Float_t lambda[2];
 976       rp->GetElipsAxis(lambda);
 977       Int_t * primaries;
 978       Int_t nprimaries;
 979       primaries = rp->GetPrimaries(nprimaries);
 980       printf("\n%6d  %8.4f  %3d     %2d     %4.1f    %4.1f %4.1f  %4.1f %4.1f %4.1f    %4.1f   %4f  %4f    %2d     : ",
 981              rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetEMCALArm(),
 982              globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(),
 983              rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ;
 984       for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
 985         printf("%d ", primaries[iprimary] ) ;
 986       }
 987     }
 988
 989     printf("\n-----------------------------------------------------------------------\n");
 990   }
 991 }