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 class TROOT;
  50 #include <TH1.h>
  51 #include <TFile.h>
  52 class TFolder;
  53 #include <TMath.h>
  54 #include <TMinuit.h>
  55 #include <TTree.h>
  56 class TSystem;
  57 #include <TBenchmark.h>
  58 #include <TBrowser.h>
  59
  60 // --- Standard library ---
  61
  62
  63 // --- AliRoot header files ---
  64 #include "AliRunLoader.h"
  65 #include "AliRun.h"
  66 #include "AliESD.h"
  67 #include "AliEMCALLoader.h"
  68 #include "AliEMCALClusterizerv1.h"
  69 #include "AliEMCALRecPoint.h"
  70 #include "AliEMCALDigit.h"
  71 #include "AliEMCALDigitizer.h"
  72 #include "AliEMCAL.h"
  73 #include "AliEMCALGeometry.h"
  74 #include "AliEMCALHistoUtilities.h"
  75 #include "AliCDBManager.h"
  76
  77 class AliCDBStorage;
  78 #include "AliCDBEntry.h"
  79
  80 ClassImp(AliEMCALClusterizerv1)
  81
  82 //____________________________________________________________________________
  83   AliEMCALClusterizerv1::AliEMCALClusterizerv1() : AliEMCALClusterizer()
  84 {
  85   // default ctor (to be used mainly by Streamer)
  86
  87   InitParameters() ;
  88   fDefaultInit = kTRUE ;
  89   fGeom = AliEMCALGeometry::GetInstance();
  90   fGeom->GetTransformationForSM(); // Global <-> Local
  91 }
  92
  93 //____________________________________________________________________________
  94 AliEMCALClusterizerv1::AliEMCALClusterizerv1(const TString alirunFileName, const TString eventFolderName)
  95 :AliEMCALClusterizer(alirunFileName, eventFolderName)
  96 {
  97   // ctor with the indication of the file where header Tree and digits Tree are stored
  98
  99   InitParameters() ;
 100   Init() ;
 101   fDefaultInit = kFALSE ;
 102 }
 103
 104 //____________________________________________________________________________
 105   AliEMCALClusterizerv1::~AliEMCALClusterizerv1()
 106 {
 107   // dtor
 108 }
 109
 110 //____________________________________________________________________________
 111 const TString AliEMCALClusterizerv1::BranchName() const
 112 {
 113    return GetName();
 114 }
 115
 116 //____________________________________________________________________________
 117 Float_t  AliEMCALClusterizerv1::Calibrate(Int_t amp, Int_t AbsId)
 118 {
 119
 120   // Convert digitized amplitude into energy.
 121   // Calibration parameters are taken from calibration data base for raw data,
 122   // or from digitizer parameters for simulated data.
 123
 124   if(fCalibData){
 125
 126     //JLK 13-Mar-2006
 127     //We now get geometry at a higher level
 128     //
 129     // Loader
 130     //    AliRunLoader *rl = AliRunLoader::GetRunLoader();
 131
 132     // Load EMCAL Geomtry
 133     //    rl->LoadgAlice();
 134     //AliRun   * gAlice = rl->GetAliRun();
 135     //AliEMCAL * emcal  = (AliEMCAL*)gAlice->GetDetector("EMCAL");
 136     //AliEMCALGeometry * geom = emcal->GetGeometry();
 137
 138     if (fGeom==0)
 139       AliFatal("Did not get geometry from EMCALLoader") ;
 140
 141     Int_t iSupMod = -1;
 142     Int_t nTower  = -1;
 143     Int_t nIphi   = -1;
 144     Int_t nIeta   = -1;
 145     Int_t iphi    = -1;
 146     Int_t ieta    = -1;
 147
 148     Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nTower, nIphi, nIeta) ;
 149     if(!bCell)
 150       Error("DigitizeEnergy","Wrong cell id number") ;
 151     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nTower,nIphi, nIeta,iphi,ieta);
 152
 153     fADCchannelECA  = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
 154     fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
 155     return -fADCpedestalECA + amp * fADCchannelECA ;
 156
 157   }
 158   else //Return energy with default parameters if calibration is not available
 159     return -fADCpedestalECA + amp * fADCchannelECA ;
 160
 161 }
 162
 163 //____________________________________________________________________________
 164 void AliEMCALClusterizerv1::Exec(Option_t * option)
 165 {
 166   // Steering method to perform clusterization for events
 167   // in the range from fFirstEvent to fLastEvent.
 168   // This range is optionally set by SetEventRange().
 169   // if fLastEvent=-1 (by default), then process events until the end.
 170
 171   if(strstr(option,"tim"))
 172     gBenchmark->Start("EMCALClusterizer");
 173
 174   if(strstr(option,"print"))
 175     Print("") ;
 176
 177   AliRunLoader *rl = AliRunLoader::GetRunLoader();
 178   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
 179
 180   //Get calibration parameters from file or digitizer default values.
 181   GetCalibrationParameters() ;
 182
 183   if (fLastEvent == -1)
 184     fLastEvent = rl->GetNumberOfEvents() - 1;
 185   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
 186
 187   Int_t ievent ;
 188   rl->LoadDigits("EMCAL");
 189   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
 190     rl->GetEvent(ievent);
 191
 192     fNumberOfECAClusters = 0;
 193
 194     if(strstr(option,"pseudo"))
 195       MakeClusters("pseudo") ;  //both types
 196     else
 197       MakeClusters("") ;  //only the real clusters
 198
 199     if(fToUnfold)
 200       MakeUnfolding() ;
 201
 202     WriteRecPoints() ;
 203
 204     if(strstr(option,"deb"))
 205       PrintRecPoints(option) ;
 206
 207     //increment the total number of recpoints per run
 208     fRecPointsInRun += emcalLoader->RecPoints()->GetEntriesFast() ;
 209   }
 210
 211   Unload();
 212
 213   if(strstr(option,"tim")){
 214     gBenchmark->Stop("EMCALClusterizer");
 215     printf("Exec took %f seconds for Clusterizing %f seconds per event",
 216          gBenchmark->GetCpuTime("EMCALClusterizer"), gBenchmark->GetCpuTime("EMCALClusterizer")/nEvents );
 217   }
 218
 219 }
 220
 221 //____________________________________________________________________________
 222 Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALRecPoint * emcRP, AliEMCALDigit ** maxAt, Float_t * maxAtEnergy,
 223                                     Int_t nPar, Float_t * fitparameters) const
 224 {
 225   // Calls TMinuit to fit the energy distribution of a cluster with several maxima
 226   // The initial values for fitting procedure are set equal to the positions of local maxima.
 227   // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
 228
 229   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 230   TClonesArray *digits = emcalLoader->Digits();
 231
 232   gMinuit->mncler();                     // Reset Minuit's list of paramters
 233   gMinuit->SetPrintLevel(-1) ;           // No Printout
 234   gMinuit->SetFCN(AliEMCALClusterizerv1::UnfoldingChiSquare) ;
 235                                          // To set the address of the minimization function
 236   TList * toMinuit = new TList();
 237   toMinuit->AddAt(emcRP,0) ;
 238   toMinuit->AddAt(digits,1) ;
 239
 240   gMinuit->SetObjectFit(toMinuit) ;         // To tranfer pointer to UnfoldingChiSquare
 241
 242   // filling initial values for fit parameters
 243   AliEMCALDigit * digit ;
 244
 245   Int_t ierflg  = 0;
 246   Int_t index   = 0 ;
 247   Int_t nDigits = (Int_t) nPar / 3 ;
 248
 249   Int_t iDigit ;
 250
 251   for(iDigit = 0; iDigit < nDigits; iDigit++){
 252     digit = maxAt[iDigit];
 253
 254     Int_t relid[2] ;
 255     Float_t x = 0.;
 256     Float_t z = 0.;
 257     fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 258     fGeom->PosInAlice(relid, x, z) ;
 259
 260     Float_t energy = maxAtEnergy[iDigit] ;
 261
 262     gMinuit->mnparm(index, "x",  x, 0.1, 0, 0, ierflg) ;
 263     index++ ;
 264     if(ierflg != 0){
 265       Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : x = %f",  x ) ;
 266       return kFALSE;
 267     }
 268     gMinuit->mnparm(index, "z",  z, 0.1, 0, 0, ierflg) ;
 269     index++ ;
 270     if(ierflg != 0){
 271        Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : z = %f", z) ;
 272       return kFALSE;
 273     }
 274     gMinuit->mnparm(index, "Energy",  energy , 0.05*energy, 0., 4.*energy, ierflg) ;
 275     index++ ;
 276     if(ierflg != 0){
 277      Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : energy = %f", energy) ;
 278       return kFALSE;
 279     }
 280   }
 281
 282   Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
 283                       //  depends on it.
 284   Double_t p1 = 1.0 ;
 285   Double_t p2 = 0.0 ;
 286
 287   gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls
 288   gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient
 289   gMinuit->SetMaxIterations(5);
 290   gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings
 291
 292   gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize
 293
 294   if(ierflg == 4){  // Minimum not found
 295     Error("FindFit", "EMCAL Unfolding  Fit not converged, cluster abandoned " ) ;
 296     return kFALSE ;
 297   }
 298   for(index = 0; index < nPar; index++){
 299     Double_t err ;
 300     Double_t val ;
 301     gMinuit->GetParameter(index, val, err) ;    // Returns value and error of parameter index
 302     fitparameters[index] = val ;
 303    }
 304
 305   delete toMinuit ;
 306   return kTRUE;
 307
 308 }
 309
 310 //____________________________________________________________________________
 311 void AliEMCALClusterizerv1::GetCalibrationParameters()
 312 {
 313   // Set calibration parameters:
 314   // if calibration database exists, they are read from database,
 315   // otherwise, they are taken from digitizer.
 316   //
 317   // It is a user responsilibity to open CDB before reconstruction,
 318   // for example:
 319   // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
 320
 321   //Check if calibration is stored in data base
 322    if(AliCDBManager::Instance()->IsDefaultStorageSet()){
 323      AliCDBEntry *entry = (AliCDBEntry*)
 324      AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
 325      if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject();
 326    }
 327    if(!fCalibData)
 328      {
 329        //If calibration is not available use default parameters
 330        //Loader
 331        AliEMCALLoader *emcalLoader =
 332          dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 333        if ( !emcalLoader->Digitizer() )
 334          emcalLoader->LoadDigitizer();
 335        AliEMCALDigitizer * dig = dynamic_cast<AliEMCALDigitizer*>(emcalLoader->Digitizer());
 336
 337        fADCchannelECA   = dig->GetECAchannel() ;
 338        fADCpedestalECA  = dig->GetECApedestal();
 339      }
 340 }
 341
 342 //____________________________________________________________________________
 343 void AliEMCALClusterizerv1::Init()
 344 {
 345   // Make all memory allocations which can not be done in default constructor.
 346   // Attach the Clusterizer task to the list of EMCAL tasks
 347
 348   AliRunLoader *rl = AliRunLoader::GetRunLoader();
 349   fGeom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
 350   fGeom->GetTransformationForSM(); // Global <-> Local
 351   AliInfo(Form("geom 0x%x",fGeom));
 352
 353   if(!gMinuit)
 354     gMinuit = new TMinuit(100) ;
 355
 356   fHists = BookHists();
 357 }
 358
 359 //____________________________________________________________________________
 360 void AliEMCALClusterizerv1::InitParameters()
 361 {
 362   // Initializes the parameters for the Clusterizer
 363   fNumberOfECAClusters = 0;
 364
 365   fNTowerInGroup = 36;  //Produces maximum of 80 pseudoclusters per event
 366
 367   fECAClusteringThreshold   = 0.5;  // value obtained from Alexei
 368   fECALocMaxCut = 0.03; // ??
 369
 370   fECAW0     = 4.5 ;
 371   fTimeGate = 1.e-8 ;
 372   fToUnfold = kFALSE ;
 373   fRecPointsInRun  = 0 ;
 374   fMinECut = 0.01; // have to be tune
 375
 376   fCalibData               = 0 ;
 377 }
 378
 379 //____________________________________________________________________________
 380 Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2) const
 381 {
 382   // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
 383   //                                       = 1 are neighbour
 384   //                                       = 2 is in different SM, quit from searching
 385   // neighbours are defined as digits having at least a common vertex
 386   // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
 387   //                                      which is compared to a digit (d2)  not yet in a cluster
 388
 389   static Int_t rv;
 390   static Int_t nSupMod1=0, nTower1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
 391   static Int_t nSupMod2=0, nTower2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
 392   static Int_t rowdiff, coldiff;
 393   rv = 0 ;
 394
 395   fGeom->GetCellIndex(d1->GetId(), nSupMod1,nTower1,nIphi1,nIeta1);
 396   fGeom->GetCellIndex(d2->GetId(), nSupMod2,nTower2,nIphi2,nIeta2);
 397   if(nSupMod1 != nSupMod2) return 2; // different SM
 398
 399   fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nTower1,nIphi1,nIeta1, iphi1,ieta1);
 400   fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nTower2,nIphi2,nIeta2, iphi2,ieta2);
 401
 402   rowdiff = TMath::Abs(iphi1 - iphi2);
 403   coldiff = TMath::Abs(ieta1 - ieta2) ;
 404
 405   if (( coldiff <= 1 )  && ( rowdiff <= 1 )) rv = 1;  // neighbours with at least commom vertex
 406
 407   if (gDebug == 2 && rv==1)
 408   printf("AreNeighbours: neighbours=%d, id1=%d, relid1=%d,%d \n id2=%d, relid2=%d,%d \n",
 409          rv, d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2);
 410
 411   return rv ;
 412 }
 413
 414 //____________________________________________________________________________
 415 Int_t AliEMCALClusterizerv1::AreInGroup(AliEMCALDigit * d1, AliEMCALDigit * d2) const
 416 {
 417   // Tells whether two digits fall within the same supermodule and
 418   // tower grouping.  The number of towers in a group is controlled by
 419   // the parameter nTowersInGroup
 420   //    = 0 are not in same group but continue searching
 421   //    = 1 same group
 422   //    = 2 is in different SM, quit from searching
 423   //    = 3 different tower group, quit from searching
 424   //
 425   // The order of d1 and d2 is important: first (d1) should be a digit
 426   // already in a cluster which is compared to a digit (d2)  not yet in a cluster
 427
 428   //JLK Question: does the quit from searching assume that the digits
 429   //are ordered, so that once you are in a different SM, you'll not
 430   //find another in the list that will match?  How about my TowerGroup search?
 431
 432   static Int_t rv;
 433   static Int_t nSupMod1=0, nTower1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
 434   static Int_t nSupMod2=0, nTower2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
 435   static Int_t towerGroup1 = -1, towerGroup2 = -1;
 436   rv = 0 ;
 437
 438   fGeom->GetCellIndex(d1->GetId(), nSupMod1,nTower1,nIphi1,nIeta1);
 439   fGeom->GetCellIndex(d2->GetId(), nSupMod2,nTower2,nIphi2,nIeta2);
 440   if(nSupMod1 != nSupMod2) return 2; // different SM
 441
 442   static Int_t nTowerInSM = fGeom->GetNCellsInSupMod()/fGeom->GetNCellsInTower();
 443
 444   //figure out which tower grouping each digit belongs to
 445   for(int it = 0; it < nTowerInSM/fNTowerInGroup; it++) {
 446     if(nTower1 <= nTowerInSM - it*fNTowerInGroup) towerGroup1 = it;
 447     if(nTower2 <= nTowerInSM - it*fNTowerInGroup) towerGroup2 = it;
 448   }
 449   if(towerGroup1 != towerGroup2) return 3; //different Towergroup
 450
 451   //same SM, same towergroup, we're happy
 452   if(towerGroup1 == towerGroup2 && towerGroup2 >= 0)
 453     rv = 1;
 454
 455   if (gDebug == 2 && rv==1)
 456   printf("AreInGroup: neighbours=%d, id1=%d, relid1=%d,%d \n id2=%d, relid2=%d,%d \n",
 457          rv, d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2);
 458
 459   return rv ;
 460 }
 461
 462 //____________________________________________________________________________
 463 void AliEMCALClusterizerv1::Unload()
 464 {
 465   // Unloads the Digits and RecPoints
 466   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 467
 468   emcalLoader->UnloadDigits() ;
 469   emcalLoader->UnloadRecPoints() ;
 470 }
 471
 472 //____________________________________________________________________________
 473 void AliEMCALClusterizerv1::WriteRecPoints()
 474 {
 475
 476   // Creates new branches with given title
 477   // fills and writes into TreeR.
 478
 479   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 480
 481   TObjArray * aECARecPoints = emcalLoader->RecPoints() ;
 482
 483   TClonesArray * digits = emcalLoader->Digits() ;
 484   TTree * treeR = emcalLoader->TreeR();
 485   if ( treeR==0 ) {
 486     emcalLoader->MakeRecPointsContainer();
 487     treeR = emcalLoader->TreeR();
 488   }
 489   Int_t index ;
 490
 491   //Evaluate position, dispersion and other RecPoint properties for EC section
 492   for(index = 0; index < aECARecPoints->GetEntries(); index++)
 493     (dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(index)))->EvalAll(fECAW0,digits) ;
 494
 495   aECARecPoints->Sort() ;
 496
 497   for(index = 0; index < aECARecPoints->GetEntries(); index++) {
 498     (dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(index)))->SetIndexInList(index) ;
 499     (dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(index)))->Print();
 500   }
 501
 502   Int_t bufferSize = 32000 ;
 503   Int_t splitlevel = 0 ;
 504
 505   //EC section branch
 506
 507   TBranch * branchECA = 0;
 508   if ((branchECA = treeR->GetBranch("EMCALECARP")))
 509     branchECA->SetAddress(&aECARecPoints);
 510   else
 511     treeR->Branch("EMCALECARP","TObjArray",&aECARecPoints,bufferSize,splitlevel);
 512   treeR->Fill() ;
 513
 514   emcalLoader->WriteRecPoints("OVERWRITE");
 515
 516 }
 517
 518 //____________________________________________________________________________
 519 void AliEMCALClusterizerv1::MakeClusters(char* option)
 520 {
 521   // Steering method to construct the clusters stored in a list of Reconstructed Points
 522   // A cluster is defined as a list of neighbour digits
 523
 524   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 525
 526   TObjArray * aECARecPoints = emcalLoader->RecPoints() ;
 527
 528   if (fGeom==0)
 529      AliFatal("Did not get geometry from EMCALLoader");
 530
 531   aECARecPoints->Clear();
 532
 533   //Start with pseudoclusters, if option
 534   if(strstr(option,"pseudo")) {
 535     TClonesArray * digits  = emcalLoader->Digits() ;
 536     TClonesArray * digitsC =  dynamic_cast<TClonesArray*>(digits->Clone());
 537
 538     TIter nextdigit(digitsC) ;
 539     AliEMCALDigit * digit;
 540
 541     AliEMCALRecPoint * recPoint = 0 ;
 542     int ineb=0;
 543     nextdigit.Reset();
 544
 545     // PseudoClusterization starts
 546     while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // scan over the list of digitsC
 547       recPoint = 0 ;
 548       TArrayI clusterECAdigitslist(1000); // what is this
 549
 550       if(fGeom->CheckAbsCellId(digit->GetId()) ) { //Is this an EMCAL digit? Just maing sure...
 551         Int_t iDigitInECACluster = 0; // start a new Tower RecPoint
 552         if(fNumberOfECAClusters >= aECARecPoints->GetSize()) aECARecPoints->Expand(2*fNumberOfECAClusters+1) ;
 553         AliEMCALRecPoint * rp = new  AliEMCALRecPoint("") ;
 554         aECARecPoints->AddAt(rp, fNumberOfECAClusters) ;
 555         recPoint = dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(fNumberOfECAClusters)) ;
 556         fNumberOfECAClusters++ ;
 557
 558         recPoint->SetClusterType(AliESDCaloCluster::kPseudoCluster);
 559
 560         recPoint->AddDigit(*digit, digit->GetAmp()) ;
 561         clusterECAdigitslist[iDigitInECACluster] = digit->GetIndexInList() ;
 562         iDigitInECACluster++ ;
 563         digitsC->Remove(digit) ;
 564         AliDebug(1,Form("MakePseudoClusters: OK id = %d, adc = %f \n", digit->GetId(), digit->GetAmp()));
 565         nextdigit.Reset(); // will start from beggining
 566
 567         AliEMCALDigit * digitN = 0; // digi neighbor
 568         Int_t index = 0 ;
 569
 570         // Find the neighbours
 571         while (index < iDigitInECACluster){ // scan over digits already in cluster
 572           digit =  (AliEMCALDigit*)digits->At(clusterECAdigitslist[index]);
 573           index++ ;
 574           while ( (digitN = (AliEMCALDigit *)nextdigit())) { // scan over the reduced list of digits
 575             ineb = AreInGroup(digit, digitN);       // call (digit,digitN) in THAT oder !!!!!
 576             switch (ineb ) {
 577               case 0 :   // not a neighbour
 578                break ;
 579               case 1 :   // are neighbours
 580                recPoint->AddDigit(*digitN, digitN->GetAmp() ) ;
 581                clusterECAdigitslist[iDigitInECACluster] = digitN->GetIndexInList() ;
 582                iDigitInECACluster++ ;
 583                digitsC->Remove(digitN) ;
 584                break ;
 585               case 2 :   // different SM
 586                break ;
 587               case 3 :   // different Tower group
 588                break ;
 589             } // switch
 590           } // scan over the reduced list of digits
 591         } // scan over digits already in cluster
 592         nextdigit.Reset() ;  // will start from beggining
 593       }
 594     }
 595     if(recPoint)
 596       AliDebug(1,Form("MakeClusters: cl.e %f \n", recPoint->GetEnergy()));
 597     //if(recPoint) cout << "cl.e " << recPoint->GetEnergy() << endl;
 598     delete digitsC ;
 599   }
 600
 601   //Now do real clusters
 602   TClonesArray * digits  = emcalLoader->Digits() ;
 603   TClonesArray * digitsC =  dynamic_cast<TClonesArray*>(digits->Clone());
 604
 605   TIter nextdigit(digitsC) ;
 606   AliEMCALDigit * digit;
 607
 608   AliEMCALRecPoint * recPoint = 0 ;
 609   int ineb=0;
 610   nextdigit.Reset();
 611
 612   double e=0.0, ehs = 0.0;
 613   while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // clean up digits
 614     e = Calibrate(digit->GetAmp(), digit->GetId());
 615     AliEMCALHistoUtilities::FillH1(fHists, 10, digit->GetAmp());
 616     AliEMCALHistoUtilities::FillH1(fHists, 11, e);
 617     if(e < fMinECut ) digitsC->Remove(digit);
 618     else              ehs += e;
 619   }
 620   AliDebug(1,Form("MakeClusters: Number of digits %d  -> (e %f), ehs %d\n",
 621                   digits->GetEntries(),fMinECut,ehs));
 622   //cout << " Number of digits " << digits->GetEntries() << " -> (e>" <<fMinECut <<")";
 623   //cout << digitsC->GetEntries()<< " ehs "<<ehs<<endl;
 624
 625   // Clusterization starts
 626   //  cout << "Outer Loop" << endl;
 627   ineb=0;
 628   nextdigit.Reset();
 629   recPoint = 0 ;
 630   while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // scan over the list of digitsC
 631     recPoint = 0 ;
 632     TArrayI clusterECAdigitslist(1000); // what is this
 633
 634     if(fGeom->CheckAbsCellId(digit->GetId()) && (Calibrate(digit->GetAmp(), digit->GetId()) > fECAClusteringThreshold  ) ){
 635       Int_t iDigitInECACluster = 0; // start a new Tower RecPoint
 636       if(fNumberOfECAClusters >= aECARecPoints->GetSize()) aECARecPoints->Expand(2*fNumberOfECAClusters+1) ;
 637       AliEMCALRecPoint * rp = new  AliEMCALRecPoint("") ;
 638       aECARecPoints->AddAt(rp, fNumberOfECAClusters) ;
 639       recPoint = dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(fNumberOfECAClusters)) ;
 640       fNumberOfECAClusters++ ;
 641
 642       recPoint->SetClusterType(AliESDCaloCluster::kClusterv1);
 643
 644       recPoint->AddDigit(*digit, Calibrate(digit->GetAmp(), digit->GetId())) ;
 645       clusterECAdigitslist[iDigitInECACluster] = digit->GetIndexInList() ;
 646       iDigitInECACluster++ ;
 647       digitsC->Remove(digit) ;
 648       AliDebug(1,Form("MakeClusters: OK id = %d, ene = %f , thre = %f \n", digit->GetId(),
 649       Calibrate(digit->GetAmp(),digit->GetId()), fECAClusteringThreshold));
 650       nextdigit.Reset(); // will start from beggining
 651
 652       AliEMCALDigit * digitN = 0; // digi neighbor
 653       Int_t index = 0 ;
 654
 655       // Find the neighbours
 656       while (index < iDigitInECACluster){ // scan over digits already in cluster
 657         digit =  (AliEMCALDigit*)digits->At(clusterECAdigitslist[index]);
 658         index++ ;
 659         while ( (digitN = (AliEMCALDigit *)nextdigit())) { // scan over the reduced list of digits
 660           ineb = AreNeighbours(digit, digitN);       // call (digit,digitN) in THAT oder !!!!!
 661           switch (ineb ) {
 662             case 0 :   // not a neighbour
 663               break ;
 664             case 1 :   // are neighbours
 665               recPoint->AddDigit(*digitN, Calibrate(digitN->GetAmp(),digitN->GetId()) ) ;
 666               clusterECAdigitslist[iDigitInECACluster] = digitN->GetIndexInList() ;
 667               iDigitInECACluster++ ;
 668               digitsC->Remove(digitN) ;
 669                break ;
 670              case 2 :   // different SM
 671                break ;
 672           } // switch
 673         } // scan over the reduced list of digits
 674       } // scan over digits already in cluster
 675       nextdigit.Reset() ;  // will start from beggining
 676     }
 677   } // while digit
 678   if(recPoint)
 679     AliDebug(1,Form("MakeClusters: cl.e %f \n", recPoint->GetEnergy()));
 680   //if(recPoint) cout << "cl.e " << recPoint->GetEnergy() << endl;
 681   delete digitsC ;
 682
 683   AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,digits->GetEntriesFast()));
 684 }
 685
 686 //____________________________________________________________________________
 687 void AliEMCALClusterizerv1::MakeUnfolding() const
 688 {
 689   Fatal("AliEMCALClusterizerv1::MakeUnfolding", "--> Unfolding not implemented") ;
 690
 691 }
 692
 693 //____________________________________________________________________________
 694 Double_t  AliEMCALClusterizerv1::ShowerShape(Double_t r)
 695 {
 696   // Shape of the shower (see EMCAL TDR)
 697   // If you change this function, change also the gradient evaluation in ChiSquare()
 698
 699   Double_t r4    = r*r*r*r ;
 700   Double_t r295  = TMath::Power(r, 2.95) ;
 701   Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
 702   return shape ;
 703 }
 704
 705 //____________________________________________________________________________
 706 void  AliEMCALClusterizerv1::UnfoldCluster(AliEMCALRecPoint * /*iniTower*/,
 707                                            Int_t /*nMax*/,
 708                                            AliEMCALDigit ** /*maxAt*/,
 709                                            Float_t * /*maxAtEnergy*/) const
 710 {
 711   // Performs the unfolding of a cluster with nMax overlapping showers
 712
 713   Fatal("UnfoldCluster", "--> Unfolding not implemented") ;
 714
 715 }
 716
 717 //_____________________________________________________________________________
 718 void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & /*nPar*/, Double_t * /*Grad*/,
 719                                                Double_t & /*fret*/,
 720                                                Double_t * /*x*/, Int_t /*iflag*/)
 721 {
 722   // Calculates the Chi square for the cluster unfolding minimization
 723   // Number of parameters, Gradient, Chi squared, parameters, what to do
 724
 725   ::Fatal("UnfoldingChiSquare","Unfolding not implemented") ;
 726 }
 727 //____________________________________________________________________________
 728 void AliEMCALClusterizerv1::Print(Option_t * /*option*/)const
 729 {
 730   // Print clusterizer parameters
 731
 732   TString message("\n") ;
 733
 734   if( strcmp(GetName(), "") !=0 ){
 735
 736     // Print parameters
 737
 738     TString taskName(GetName()) ;
 739     taskName.ReplaceAll(Version(), "") ;
 740
 741     printf("--------------- ");
 742     printf(taskName.Data()) ;
 743     printf(" ");
 744     printf(GetTitle()) ;
 745     printf("-----------\n");
 746     printf("Clusterizing digits from the file: ");
 747     printf(taskName.Data());
 748     printf("\n                           Branch: ");
 749     printf(GetName());
 750     printf("\n                       ECA Local Maximum cut    = %f", fECALocMaxCut);
 751     printf("\n                       ECA Logarithmic weight   = %f", fECAW0);
 752     if(fToUnfold)
 753       printf("\nUnfolding on\n");
 754     else
 755       printf("\nUnfolding off\n");
 756
 757     printf("------------------------------------------------------------------");
 758   }
 759   else
 760     printf("AliEMCALClusterizerv1 not initialized ") ;
 761 }
 762
 763 //____________________________________________________________________________
 764 void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option)
 765 {
 766   // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
 767   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 768   TObjArray * aECARecPoints = emcalLoader->RecPoints() ;
 769
 770   printf("PrintRecPoints: Clusterization result:") ;
 771
 772   printf("event # %d\n", emcalLoader->GetRunLoader()->GetEventNumber() ) ;
 773   printf("           Found %d ECA Rec Points\n ",
 774          aECARecPoints->GetEntriesFast()) ;
 775
 776   fRecPointsInRun +=  aECARecPoints->GetEntriesFast() ;
 777
 778   if(strstr(option,"all")) {
 779     Int_t index =0;
 780     printf("\n-----------------------------------------------------------------------\n") ;
 781     printf("Clusters in ECAL section\n") ;
 782     printf("Index    Ene(GeV) Multi Module     phi     r   theta    X    Y      Z   Dispersion Lambda 1   Lambda 2  # of prim  Primaries list\n") ;
 783    Float_t maxE=0;
 784    Float_t maxL1=0;
 785    Float_t maxL2=0;
 786    Float_t maxDis=0;
 787
 788     for (index = 0 ; index < aECARecPoints->GetEntries() ; index++) {
 789       AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint * >(aECARecPoints->At(index)) ;
 790       TVector3  globalpos;
 791       //rp->GetGlobalPosition(globalpos);
 792       TVector3  localpos;
 793       rp->GetLocalPosition(localpos);
 794       Float_t lambda[2];
 795       rp->GetElipsAxis(lambda);
 796       Int_t * primaries;
 797       Int_t nprimaries;
 798       primaries = rp->GetPrimaries(nprimaries);
 799       printf("\n%6d  %8.4f  %3d     %4.1f    %4.1f %4.1f  %4.1f %4.1f %4.1f    %4.1f   %4f  %4f    %2d     : ",
 800              rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(),
 801              globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(),
 802              rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ;
 803   /////////////
 804       if(rp->GetEnergy()>maxE){
 805               maxE=rp->GetEnergy();
 806               maxL1=lambda[0];
 807               maxL2=lambda[1];
 808               maxDis=rp->GetDispersion();
 809       }
 810       fPointE->Fill(rp->GetEnergy());
 811       fPointL1->Fill(lambda[0]);
 812       fPointL2->Fill(lambda[1]);
 813       fPointDis->Fill(rp->GetDispersion());
 814       fPointMult->Fill(rp->GetMultiplicity());
 815       /////////////
 816       for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
 817         printf("%d ", primaries[iprimary] ) ;
 818       }
 819     }
 820
 821       fMaxE->Fill(maxE);
 822       fMaxL1->Fill(maxL1);
 823       fMaxL2->Fill(maxL2);
 824       fMaxDis->Fill(maxDis);
 825
 826
 827     printf("\n-----------------------------------------------------------------------\n");
 828   }
 829 }
 830 TList* AliEMCALClusterizerv1::BookHists()
 831 {
 832   //set up histograms for monitoring clusterizer performance
 833
 834   gROOT->cd();
 835
 836         fPointE = new TH1F("pointE","point energy", 2000, 0.0, 150.);
 837         fPointL1 = new TH1F("pointL1","point L1", 1000, 0.0, 3.);
 838         fPointL2 = new TH1F("pointL2","point L2", 1000, 0.0, 3.);
 839         fPointDis = new TH1F("pointDis","point Dis", 1000, 0.0, 3.);
 840         fPointMult = new TH1F("pointMult","point Mult", 100, 0.0, 100.);
 841         fDigitAmp = new TH1F("digitAmp","Digit Amplitude", 2000, 0.0, 5000.);
 842         fMaxE = new TH1F("maxE","Max point energy", 2000, 0.0, 150.);
 843         fMaxL1 = new TH1F("maxL1","Max point L1", 1000, 0.0, 3.);
 844         fMaxL2 = new TH1F("maxL2","Max point L2", 1000, 0.0, 3.);
 845         fMaxDis = new TH1F("maxDis","Max point Dis", 1000, 0.0, 3.); // 9
 846         //
 847         new TH1F("adcOfDigits","adc of digits(threshold control)", 1001, -0.5, 1000.5);   // 10
 848         new TH1F("energyOfDigits","energy of digits(threshold control)", 1000, 0.0, 1.);  // 11
 849
 850   return AliEMCALHistoUtilities::MoveHistsToList("EmcalClusterizerv1ControlHists", kFALSE);
 851 }
 852
 853 void AliEMCALClusterizerv1::SaveHists(const char *fn)
 854 {
 855   AliEMCALHistoUtilities::SaveListOfHists(fHists, fn, kTRUE);
 856 }
 857
 858 void AliEMCALClusterizerv1::Browse(TBrowser* b)
 859 {
 860   if(fHists) b->Add(fHists);
 861   TTask::Browse(b);
 862 }