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