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