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